/* 12.010 Homework #3 Queston 1:
  Gamma function calculation

Compile and link with
cc -lm HW03_01.c -o HW03_01

Methods are based on:

M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions with 
Formulas, Graphs and Mathematical Tables, Wiley, New York, 1970.
(Gamma finctons: page 255)
Web access:
http://www.convertit.com/Go/Convertit/Reference/AMS55.ASP
http://www.convertit.com/Go/Convertit/Reference/AMS55.ASP?Res=150&Page=255

Basic methosd here is explicit formualas of Gamma(n) = (n-1)!, 
gamma(n+1/3, 1/2 and 2/3).  For the -1 to +1 ranges, Euler's series is used
*/

/* Prototype definitions */
double factorial( int ) ;
double gammai( int ) ; 
double gammathird( int );
double gammatwothird( int );
double gammahalf( int );
double gammaeul( double , double );
double gammainf( double , double );
double gammaser( double , double );

#include <math.h>
#include <stdio.h>

int main(argc, argv) 
int argc ; // Number of arguments
char  *argv[] ; // Pointer to arguments

{
double eps = 1.e-6 ;  // Accuracy of calculation (can be passed in runstring)
double z ; // non-integer argument for gamma function
double gamma_1, gamma_2, gamma_3, gamma_4; // Methods for computing or values of gamma
int i, j ; // Loop variable

char line[80];  // Output line which will include strings and number

/* See if tolerance is passed in runsting */
 if ( argc > 1 ) {
    sscanf(argv[1],"%lf", &eps);
 }

/* Output header information */
 printf("\nTable of gamma functions of positive n integers\n");
 printf(  "-----------------------------------------------\n");
 printf(  "   n      gamma(n)    gamma(n+1/3)     gamma(n+1/2)     gamma(n+2/3)\n");

/* Loop over integer arguments */
 for ( i = 1 ; i <= 10 ; i++ ) 
 {
      gamma_1 = gammai(i);
      gamma_2 = gammathird(i);
      gamma_3 = gammahalf(i);
      gamma_4 = gammatwothird(i);
      printf("%4d  %12.0f %15.6f  %15.6f  %15.6f\n",i,gamma_1, gamma_2, gamma_3, gamma_4);
 }

/* Now do second fractional part of table */

 printf("\nTable of gamma functions for non-integer arguments\n"); 
 printf("--------------------------------------------------\n");
 printf("Three methods are used here: \n");
 printf("Eulers formula              GammaEul\n"); 
 printf("Eulers infinite product     GammaInf\n"); 
 printf("Asymptotic series expansion GammaSer \n");
 printf("Tolerance on calculation is %12.3e\n",eps); 
 printf("      z       GammaEul(z)         GammaInf(z)          GammaSer(z)\n");

/* Now generate results */

 for ( i = -10 ; i <= 10 ; i++ ){
     z = i/ (double) 10;
     gamma_1 = gammaeul(z,eps);
     gamma_2 = gammainf(z,eps);
     gamma_3 = gammaser(z,eps);

    /* Format the output */
    /* Fill line with blanks */
     for ( j = 0; j < 80 ; j++ ) { line[j] = 32; }

     if( gamma_1 == 0 ) {
        sprintf(line,"      %c%s",(char) 177,"Infinity");
     } else {
        sprintf(line," %14.6f",gamma_1);
     }

     if( gamma_2 == 0 ) {
        sprintf(&line[21],"      %c%s",(char) 177,"Infinity");
     } else {
        sprintf(&line[21]," %14.6f",gamma_2);
     }

     if( gamma_3 == 0 ) {
        sprintf(&line[41],"      %c%s",(char) 177,"Infinity");
     } else {
        sprintf(&line[41]," %14.6f",gamma_3);
     }
     for ( j = 0; j < 79 ; j++ ) { 
        if( line[j] == 0) { line[j] = 32 ; }
     }
     line[61] = 0;

/*   Now output line */
     printf("%8.2f  %s \n",z,line);
  }
}

   
double gammai( int n ) {
/* Function to return (n-1)! = gamma(n) */
double gammac ; /* Computed value of gamma */
int i;

 gammac = 1;
 if( n < 1 ) {
    printf("Invalid argument for gammai.  Passed value was %d\n",n);
    gammac = 0.0; 
    return(gammac);
 }

 for ( i = 2 ; i <= n-1; i++ ) {
    gammac = gammac * i ;
 }
 return(gammac);
}

double gammathird( int n ) {
/* Function to return gamma(n+1/3) */
double gammac = 2.6789385347077479 ; // Initial value
int i ;

  for ( i = 4; i <= (3*n-2); i = i + 3 ) {
      gammac = gammac * i ;
  }
  gammac = gammac / pow(3.0,(double) n);
  return( gammac );
}

double gammatwothird( int n ) {
/* Function to return gamma(n+2/3) */
double gammac = 1.3541179394264005 ; // Initial value
int i ;

  for ( i = 2; i <= (3*n-1); i = i + 3 ) {
      gammac = gammac * i ;
  }
  gammac = gammac / pow(3.0,(double) n);
  return( gammac );
}

double gammahalf( int n ) {
/* Function to return gamma(n+1/2) */
double gammac = 1.7724538509055161 ; // Initial value
int i ;

  for ( i = 3; i <= (2*n-1); i = i + 2 ) {
      gammac = gammac * i ;
  }
  gammac = gammac / pow(2.0,(double) n);
  return( gammac );
}

double gammaeul( double z, double eps ) {
/*  Function to return gamma(z) where z is any non-negative integer values
 The expansion here is based on Euler's formulas:
 gamma(z) = limit(n -> infinity) n! n^z/(z*(z+1)*(z_2)*...(z+n) 
 for all z except 0 and negative integers.  The restriction can be clearly
 seen because for these values of z we would divide by zero and gamma at
 these valuse is infinite

 To determine the size of n needed the above is re-arranged so that each
 term can be added i.e.,
 gamma_n = gamma_(n-1) *n/(z+n) * n**z/(n-1)**z
 by looking at second part we can see if enough terms have been included
*/
double gamma_n, gamma_prev ; // Iterated values
double err ;  // Error in current iteration
double npowz, nm1powz ;  // n and n-1 to power z 
int n ;   // Counter
int max_n =  10000000 ; // Maximum n allowed

/* see if negative integer argument 
  This code z + (int) (abs(z)+0.5) is like z-nint(z) in F77 for z < 0 */
  if( z <= 0 && fabs(z+ (int) (fabs(z)+0.5)) < eps ) { 
     /* printf("z close to -1, -2 , .. %f\n",z); */
     return ( 0.0 ) ;
  }

  n = 1 ;
  gamma_prev = 1.0/(z*(z+1));
  err = 1.0;
  nm1powz = 1.0 ;   // (n-1)^z for n = 1
  
  while ( n < max_n &&  err > eps/1.e6) { 
      n++ ; 
      npowz = pow((double) n,z) ; 
      gamma_n = gamma_prev * n/(z+n) * (npowz/nm1powz) ; 
      nm1powz = npowz  ;  // Save value for next iteration
      err = fabs(gamma_n - gamma_prev);
      gamma_prev = gamma_n ;
   }
   return(gamma_n);
}

double gammainf( double z, double eps ) {
/* Function to return gamma(z) where z is any non-negative integer values
   The expansion here is based on Euler's Infinit product formula of 1/gamma:
   1/gamma(z) = z*exp(gammacon*z)*product_n[(1+z/n)exp(-z/n)]
   for all z. When 1/gamma is zero the gamma is infinite. Gammma_con is the
   constant gamma = 0.5772 15664 ..

   To determine the size of n needed is determined by looking at the change
   at step.  The series is very slow to converge.
*/
double gammacon, gamma_prev ;  // Iterated values of gamma, gamma_n and gamma_(n-1)
double gammainv; //  Inverse of gamma
double factor ;  
double err ; // Error in current iteration
int n ; // Counter
int max_n = 10000000 ;   // Max n allowed (to stop iteration)

/* see if negative integer argument 
  This code z + (int) (abs(z)+0.5) is like z-nint(z) in F77 for z < 0 */
  if( z <= 0 && fabs(z+ (int) (fabs(z)+0.5)) < eps ) { 
     /* printf("z close to -1, -2 , .. %f\n",z); */
     return ( 0.0 ) ;
  }

  n = 0 ; 
  gammacon = (double) 0.57721566390153218606; 
  gamma_prev = z*exp(gammacon*z) ; 
  err = 1 ; 
  while ( n < max_n && err > eps/1.e6 ) {
     n++ ; 
     gammainv = gamma_prev*(1+z/n)*exp(-z/n); 
     err = fabs(gammainv-gamma_prev) ;
     gamma_prev = gammainv ; 
  }
/* Take innverse and return value */
  gammainv = 1/gammainv ;
  return ( gammainv ) ;
}


double gammaser( double z, double eps) {
/*
  Function to return gamma(z) where z is any non-negative integer values.
  This is a series expansion developed in the 1930s for 1/gamma(z)
  The coefficients c k are from H. T. Davis, Tables of
  higher mathematical functions, 2 vols., Principia Press,
  Bloomington, Ind., 1933, 1935 (with permission) ; with
  corrections due to H. E. Salaer.
  Values from
  http://www.convertit.com/Go/Convertit/Reference/AMS55.ASP?Res=150&Page=256
*/
int n ; // Counter 
double gammainv ;  // Inverse of gamma computed from series
double coeff[26] = { 
        1.0000000000000000,
        0.5772156649015329,
       -0.6558780715202538,
       -0.0420026350340952,
        0.1665386113822915,
       -0.0421977345555443,
       -0.0096219715278770,
        0.0072189432466630,
       -0.0011651675918591,
       -0.0002152416741149,
        0.0001280502823882,
       -0.0000201348547807,
       -0.0000012504934821,
        0.0000011330272320,
       -0.0000002056338417,
        0.0000000061160950,
        0.0000000050020075,
       -0.0000000011812746,
        0.0000000001043427,
        0.0000000000077823,
       -0.0000000000036868,
        0.0000000000005100,
       -0.0000000000000206,
       -0.0000000000000054,
        0.0000000000000014,
        0.0000000000000001 } ;  /* 26-coefficients in series */

/* see if negative integer argument 
  This code z + (int) (abs(z)+0.5) is like z-nint(z) in F77 for z < 0 */
  if( z <= 0 && fabs(z+ (int) (fabs(z)+0.5)) < eps ) { 
     /* printf("z close to -1, -2 , .. %f\n",z); */
     return ( 0.0 ) ;
  }

/*  Sum up the 26 terms in the series.  */
  gammainv = 0.0;
  for ( n = 1 ; n <= 26 ; n++ ) { 
     gammainv = gammainv + coeff[n-1]*pow(z, (double) n) ;
  }

  gammainv = 1/gammainv;
  return(gammainv);
}