////////////////////////////////////////////////////////////////////////////////
// File: testdoublefactorial.c                                                //
//                                                                            //
// Test the functions Double_Factorial(), Double_Factorial_Max_Arg(),  and    //
// xDouble_Factorial() in the file double_factorial.c.                        //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

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

extern double Double_Factorial( int n );
extern long double xDouble_Factorial( int n );
extern int    Double_Factorial_Max_Arg( void );

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

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

   printf("\nTest xDouble_Factorial() \n\n");
   printf("Exceptional values: Argument <= -2, return 0,\n");
   printf("                    Argument > Double_Factorial_Max_Arg() return DBL_MAX\n\n");
   printf("(-2)!!  %22.16Le\n",xDouble_Factorial(-2));
   printf("(%3d)!! %22.16Le\n",Double_Factorial_Max_Arg()+1,xDouble_Factorial(Double_Factorial_Max_Arg() + 1));  
   printf("(%3d)!! %22.16Le\n",Double_Factorial_Max_Arg()+2,xDouble_Factorial(Double_Factorial_Max_Arg() + 2));  
   printf("\n\nxDouble Factorials\n");
   printf("   n!!   xDouble_Factorial(n) n*xDouble_Factorial(n-2)   Relative Error /\n");
   printf("                                                          DBL_EPSILON\n");
   for (i = -1; i <= Double_Factorial_Max_Arg(); i++) {
      if (i <= 0) xz = 1.0L;
      else if (i <= 2) xz = (long double) i;
      else xz = xDouble_Factorial(i-2) * (long double) i;
      printf("%4d!!  %22.16Le %22.16Le %22.16Le\n",i,xDouble_Factorial(i),xz,(xDouble_Factorial(i)-xz)/(xDouble_Factorial(i)*(long double)DBL_EPSILON));
   }
}