Task

Implement one algorithm (or more) to compute the Gamma (



Γ


{\displaystyle \Gamma }

) function (in the real field only).

If your language has the function as built-in or you know a library which has it, compare your implementation's results with the results of the built-in/library function.

The Gamma function can be defined as:

  



Γ
(
x
)
=


∫

0


∞



t

x
−
1



e

−
t


d
t



{\displaystyle \Gamma (x)=\displaystyle \int _{0}^{\infty }t^{x-1}e^{-t}dt}


This suggests a straightforward (but inefficient) way of computing the 



Γ


{\displaystyle \Gamma }

 through numerical integration.



Better suggested methods:

 Lanczos approximation
 Stirling's approximation


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_sf_gamma.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
 
/* very simple approximation */
double st_gamma(double x)
{
  return sqrt(2.0*M_PI/x)*pow(x/M_E, x);
}
 
#define A 12
double sp_gamma(double z)
{
  const int a = A;
  static double c_space[A];
  static double *c = NULL;
  int k;
  double accm;
 
  if ( c == NULL ) {
    double k1_factrl = 1.0; /* (k - 1)!*(-1)^k with 0!==1*/
    c = c_space;
    c[0] = sqrt(2.0*M_PI);
    for(k=1; k < a; k++) {
      c[k] = exp(a-k) * pow(a-k, k-0.5) / k1_factrl;
	  k1_factrl *= -k;
    }
  }
  accm = c[0];
  for(k=1; k < a; k++) {
    accm += c[k] / ( z + k );
  }
  accm *= exp(-(z+a)) * pow(z+a, z+0.5); /* Gamma(z+1) */
  return accm/z;
}
 
int main()
{
  double x;
 
 
  printf("%15s%15s%15s%15s\n", "Stirling", "Spouge", "GSL", "libm");
  for(x=1.0; x <= 10.0; x+=1.0) {
    printf("%15.8lf%15.8lf%15.8lf%15.8lf\n", st_gamma(x/3.0), sp_gamma(x/3.0), 
	   gsl_sf_gamma(x/3.0), tgamma(x/3.0));
  }
  return 0;
}