////////////////////////////////////////////////////////////////////////////////
// File: testrisingfactorial.c                                                //
//                                                                            //
// Test the functions Rising_Factorial() and xRising_Factorial() in the file  //
// rising_factorial.c                                                         //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

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

extern double Rising_Factorial( int n, int m );
extern long double xRising_Factorial( int n, int m );
extern int    Factorial_Max_Arg( void );

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

   printf("Exceptional values\n");
   printf("(-1)0  %22.16e\n",Rising_Factorial(-1,0));
   printf("\n\nRising_Factorials\n");
   printf("  (n)m    Rising_Factorial(n,m)    (n)x...x(n+m-1)        Relative Error /\n");
   printf("                                                           DBL_EPSILON\n");
   for (i = 0; i <= 10; i++) {
      for (j = 0; j < 4; j++) {
         x = Rising_Factorial(i,j);
         z = (i + j) ? 1.0 : 0.0;
         for (k = i; k < i+j; k++) z *= k;
         if (x > 0.0)  
            printf("(%4d)%1d  %22.16e %22.16e %22.16e\n",i,j,x,z,(x-z)/(x*DBL_EPSILON));
         else 
            printf("(%4d)%1d  %22.16e %22.16e\n",i,j,x,z);
      }
   }
   printf("\n\nLarge Factorials\n");
   printf("  (n)m    Rising_Factorial(n,m)    (n)x...x(n+m-1)        Relative Error /\n");
   printf("                                                           DBL_EPSILON\n");
   for (i = Factorial_Max_Arg() + 1; i < Factorial_Max_Arg() + 3; i++) {
      for (j = 0; j < 4; j++) {
         x = Rising_Factorial(i,j);
         z = 1.0;
         for (k = i; k < i+j; k++) z *= k; 
         printf("(%4d)%1d  %22.16e %22.16e %22.16e\n",i,j,x,z,(x-z)/(x*DBL_EPSILON));
      }
   }

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

   printf("Exceptional values\n");
   printf("(-1)0  %22.16Le\n",xRising_Factorial(-1,0));
   printf("\n\nxRising_Factorials\n");
   printf("  (n)m   xRising_Factorial(n,m)    (n)x...x(n+m-1)        Relative Error /\n");
   printf("                                                           DBL_EPSILON\n");
   for (i = 0; i <= 10; i++) {
      for (j = 0; j < 4; j++) {
         xx = xRising_Factorial(i,j);
         xz = (i + j) ? 1.0L : 0.0L;
         for (k = i; k < i+j; k++) xz *= (long double) k;
         if (xx > 0.0L) 
            printf("(%4d)%1d  %22.16Le %22.16Le %22.16Le\n",i,j,xx,xz,(xx-xz)/(xx*(long double)DBL_EPSILON));
         else
            printf("(%4d)%1d  %22.16Le %22.16Le\n",i,j,xx,xz);
      }
   }
   printf("\n\nLarge Factorials\n");
   printf("  (n)m   xRising_Factorial(n,m)    (n)x...x(n+m-1)        Relative Error /\n");
   printf("                                                           DBL_EPSILON\n");
   for (i = Factorial_Max_Arg() + 1; i < Factorial_Max_Arg() + 3; i++) {
      for (j = 0; j < 4; j++) {
         xx = xRising_Factorial(i,j);
         xz = 1.0L;
         for (k = i; k < i+j; k++) xz *= (long double)k; 
         printf("(%4d)%1d  %22.16Le %22.16Le %22.16Le\n",i,j,xx,xz,(xx-xz)/(xx*(long double)DBL_EPSILON));
      }
   }
}