////////////////////////////////////////////////////////////////////////////////
// File: testfactorial.c                                                      //
//                                                                            //
// Test the functions Factorial(), Factorial_Max_Arg(), and xFactorial()      //
// in the file factorial.c.                                                   //
//                                                                            //
// A table of factorials can be found in "Handbook of Mathematical Functions  //
// with Formulas, Graphs, and Mathematical Tables" by Abramowitz and Stegun   //
// on page 272 Table 6.3 (note the table is the Gamma Function).              //
////////////////////////////////////////////////////////////////////////////////

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

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

int main()
{
   long double xz;
   double z;
   int i;
  
   printf("Program: testfactorial.c\n");
   printf("Test Factorial() \n\n");

   printf("Exceptional values: Negative argument return 0,\n");
   printf("                    Argument > Factorial_Max_Arg() return DBL_MAX\n\n");
   printf("(-1)!  %22.16e\n",Factorial(-1));
   printf("(%3d)! %22.16e\n",Factorial_Max_Arg()+1,Factorial(Factorial_Max_Arg() + 1));  
   printf("(%3d)! %22.16e\n",Factorial_Max_Arg()+2,Factorial(Factorial_Max_Arg() + 2));  
   printf("\n\nFactorials\n");
   printf("   n!       Factorial(n)         n*Factorial(n-1)      Relative Error /\n");
   printf("                                                         DBL_EPSILON\n");
   for (i = 0; i <= Factorial_Max_Arg(); i++) {
      if (i == 0) z = 1.0;
      else z = Factorial(i-1) * i;
      printf("%4d!  %22.16e %22.16e %22.16e\n",i,Factorial(i),z,(Factorial(i)-z)/(Factorial(i)*DBL_EPSILON));
   }

   printf("\n\nTest xFactorial() \n\n");

   printf("Exceptional values: Negative argument return 0,\n");
   printf("                    Argument > Factorial_Max_Arg() return DBL_MAX\n\n");
   printf("(-1)!  %22.16e\n",(double) xFactorial(-1));
   printf("(%3d)! %22.16e\n",Factorial_Max_Arg()+1,(double)xFactorial(Factorial_Max_Arg() + 1));  
   printf("(%3d)! %22.16e\n",Factorial_Max_Arg()+2,(double)xFactorial(Factorial_Max_Arg() + 2));  
   printf("\n\nFactorials\n");
   printf("   n!      xFactorial(n)        n*xFactorial(n-1)      Relative Error /\n");
   printf("                                                         DBL_EPSILON\n");
   for (i = 0; i <= Factorial_Max_Arg(); i++) {
      if (i == 0) xz = 1.0L;
      else xz = xFactorial(i-1) * i;
      printf("%4d!  %22.16e %22.16e %22.16e\n",i,(double)xFactorial(i),(double)xz,(double)((xFactorial(i)-xz)/(xFactorial(i)*DBL_EPSILON)));
   }
}