////////////////////////////////////////////////////////////////////////////////
// File: testlnfactorial.c                                                    //
//                                                                            //
// Test the functions Ln_Factorial() and xLn_Factorial() in the file          //
// ln_factorial.c.                                                            //
//                                                                            //
// A table of ln factorials can be found in "Handbook of Mathematical         //
// Functions with Formulas, Graphs, and Mathematical Tables" by Abramowitz    //
// and Stegun on page 274 Table 6.4 (note the table is the log Gamma Function,//
// where the log is the common log i.e. log10).                               //
////////////////////////////////////////////////////////////////////////////////

#include <math.h>
#include <float.h>
#include <stdio.h>
                                    
                        // Externally defined routines

extern      double Ln_Factorial( int n );
extern long double xLn_Factorial( int n );
extern      double Factorial( int n );
extern long double xFactorial( int n );
extern      int    Factorial_Max_Arg( void );

int main()
{
   double z;
   double x;
   long double xx, xz;
   int i,j,k;
  
   printf("Program: testlnfactorial.c\n");
   printf("Test Ln_Factorial() \n\n");

   printf("Exceptional value: Negative argument return -DBL_MAX,\n");
   printf("ln(-1)!  %22.16e\n",Ln_Factorial(-1));
   printf("\n\nLn_Factorials\n");
   printf("   n!      Ln_Factorial(n)       ln(Factorial(n))        Relative Error /\n");
   printf("                                                         DBL_EPSILON\n");
   for (i = 0; i <= Factorial_Max_Arg(); i++) {
      z = log(Factorial(i));
      x = Ln_Factorial(i);
      if ( x != 0.0 )
         printf("%4d!  %22.16e %22.16e %22.16e\n",i,x,z,(x-z)/(x*DBL_EPSILON));
      else
         printf("%4d!  %22.16e %22.16e\n",i,x,z);
   }
   printf("\n\nLarge Factorials\n");
   for (i = Factorial_Max_Arg() + 1; i < Factorial_Max_Arg() + 3; i++) {
      z = log(Factorial(Factorial_Max_Arg()));
      k = 1;
      for (j = Factorial_Max_Arg() + 1; j <= i; j++) k *= j;
      z += log((double)k);
      x = Ln_Factorial(i);
      printf("%4d!  %22.16e %22.16e %22.16e\n",i,x,z,(x-z)/(x*DBL_EPSILON));
   }
      
   printf("\n\nTest xLn_Factorial() \n\n");

   printf("Exceptional value: Negative argument return -LDBL_MAX,\n");
   printf("ln(-1)!  %22.16Le\n",xLn_Factorial(-1));
   printf("\n\nxLn_Factorials\n");
   printf("   n!      xLn_Factorial(n)      ln(xFactorial(n))      Relative Error /\n");
   printf("                                                        LDBL_EPSILON\n");
   for (i = 0; i <= Factorial_Max_Arg(); i++) {
      xz = logl(xFactorial(i));
      xx = xLn_Factorial(i);
      if ( xx != 0.0L )
         printf("%4d!  %22.16Le %22.16Le %22.16Le\n",i,xx,xz,(xx-xz)/(xx*LDBL_EPSILON));
      else
         printf("%4d!  %22.16Le %22.16Le\n",i,xx,xz);
   }
   printf("\n\nLarge Factorials\n");
   for (i = Factorial_Max_Arg() + 1; i < Factorial_Max_Arg() + 3; i++) {
      xz = logl(xFactorial(Factorial_Max_Arg()));
      k = 1;
      for (j = Factorial_Max_Arg() + 1; j <= i; j++) k *= j;
      xz += logl((long double)k);
      xx = xLn_Factorial(i);
      printf("%4d!  %22.16Le %22.16Le %22.16Le\n",i,xx,xz,(xx-xz)/(xx*LDBL_EPSILON));
   }
}