////////////////////////////////////////////////////////////////////////////////
// File: testgammatable.c                                                     //
//                                                                            //
// Test the functions Gamma_Function() and xGamma_Function() in the file      //
// gamma_function.c                                                           //
//                                                                            //
// The file GammaTable.dat contains a list of arguments follwed by the value  //
// of the gamma function evaluated at the corresponding argument.             //
////////////////////////////////////////////////////////////////////////////////
#include <math.h>
#include <float.h>
#include <stdio.h>

extern long double xGamma_Function(long double x);
extern double Gamma_Function(double x);

FILE *infile;

int main() {

   double est;
   long double x,g;
   long double estL;
   int i;
   
   infile = fopen("GammaTable.dat","r");

   printf("Table of Positive Integral Values for the Gamma Function (double)\n");
   printf("    x           gamma                 est             Relative Err\n");
   for (i = 1; i <= 1755; i++) {
      fscanf(infile,"%6Lf %54Le",&x,&g);
      if (i > 171) continue;
      est = Gamma_Function((double)x);
      printf("%6.1Lf  %+18.8Le   %+18.8e   %+18.8e\n",x,g,est,(double)((est-g)/g));
   }
   printf("\n\n");
   printf("Table of Positive Half-Integral Values for the Gamma Function (double)\n");
   printf("    x           gamma                 est             Relative Err\n");
   for (i = 0; i <= 1755; i++) {
      fscanf(infile,"%6Le %54Le",&x,&g);
      if (i > 171) continue;
      est = Gamma_Function((double)x);
      printf("%6.1Lf  %+18.8Le   %+18.8e   %+18.8e\n",x,g,est,(double)((est-g)/g));
   }
   fclose(infile);

   infile = fopen("GammaTable.dat","r");

   printf("\n\n");
   printf("Table of Positive Integral Values for the Gamma Function (long double)\n");
   printf("    x           gamma                 est             Relative Err\n");
   for (i = 1; i <= 1755; i++) {
      fscanf(infile,"%6Lf %54Le",&x,&g);
      estL = xGamma_Function(x);
      printf("%6.1Lf  %+18.8Le   %+18.8Le   %+18.8Le\n",x,g,estL,(estL-g)/g);
   }
   printf("\n\n");
   printf("Table of Positive Half-Integral Values for the Gamma Function (long double)\n");
   printf("    x           gamma                 est             Relative Err\n");
   for (i = 0; i <= 1755; i++) {
      fscanf(infile,"%6Le %54Le",&x,&g);
      estL = xGamma_Function(x);
      printf("%6.1Lf  %+18.8Le   %+18.8Le   %+18.8Le\n",x,g,estL,(estL-g)/g);
   }
   fclose(infile);
}