blob: c466f43407c53cbd9e6083ee08c5c42a3f86997c (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
|
#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;
}
|
