1*3e12c5d1SDavid du Colombier #include <math.h>
2*3e12c5d1SDavid du Colombier #include <errno.h>
3*3e12c5d1SDavid du Colombier
4*3e12c5d1SDavid du Colombier /*
5*3e12c5d1SDavid du Colombier C program for floating point log gamma function
6*3e12c5d1SDavid du Colombier
7*3e12c5d1SDavid du Colombier gamma(x) computes the log of the absolute
8*3e12c5d1SDavid du Colombier value of the gamma function.
9*3e12c5d1SDavid du Colombier The sign of the gamma function is returned in the
10*3e12c5d1SDavid du Colombier external quantity signgam.
11*3e12c5d1SDavid du Colombier
12*3e12c5d1SDavid du Colombier The coefficients for expansion around zero
13*3e12c5d1SDavid du Colombier are #5243 from Hart & Cheney; for expansion
14*3e12c5d1SDavid du Colombier around infinity they are #5404.
15*3e12c5d1SDavid du Colombier
16*3e12c5d1SDavid du Colombier Calls log and sin.
17*3e12c5d1SDavid du Colombier */
18*3e12c5d1SDavid du Colombier
19*3e12c5d1SDavid du Colombier int signgam;
20*3e12c5d1SDavid du Colombier static double goobie = 0.9189385332046727417803297;
21*3e12c5d1SDavid du Colombier static double pi = 3.1415926535897932384626434;
22*3e12c5d1SDavid du Colombier
23*3e12c5d1SDavid du Colombier #define M 6
24*3e12c5d1SDavid du Colombier #define N 8
25*3e12c5d1SDavid du Colombier static double p1[] = {
26*3e12c5d1SDavid du Colombier 0.83333333333333101837e-1,
27*3e12c5d1SDavid du Colombier -.277777777735865004e-2,
28*3e12c5d1SDavid du Colombier 0.793650576493454e-3,
29*3e12c5d1SDavid du Colombier -.5951896861197e-3,
30*3e12c5d1SDavid du Colombier 0.83645878922e-3,
31*3e12c5d1SDavid du Colombier -.1633436431e-2,
32*3e12c5d1SDavid du Colombier };
33*3e12c5d1SDavid du Colombier static double p2[] = {
34*3e12c5d1SDavid du Colombier -.42353689509744089647e5,
35*3e12c5d1SDavid du Colombier -.20886861789269887364e5,
36*3e12c5d1SDavid du Colombier -.87627102978521489560e4,
37*3e12c5d1SDavid du Colombier -.20085274013072791214e4,
38*3e12c5d1SDavid du Colombier -.43933044406002567613e3,
39*3e12c5d1SDavid du Colombier -.50108693752970953015e2,
40*3e12c5d1SDavid du Colombier -.67449507245925289918e1,
41*3e12c5d1SDavid du Colombier 0.0,
42*3e12c5d1SDavid du Colombier };
43*3e12c5d1SDavid du Colombier static double q2[] = {
44*3e12c5d1SDavid du Colombier -.42353689509744090010e5,
45*3e12c5d1SDavid du Colombier -.29803853309256649932e4,
46*3e12c5d1SDavid du Colombier 0.99403074150827709015e4,
47*3e12c5d1SDavid du Colombier -.15286072737795220248e4,
48*3e12c5d1SDavid du Colombier -.49902852662143904834e3,
49*3e12c5d1SDavid du Colombier 0.18949823415702801641e3,
50*3e12c5d1SDavid du Colombier -.23081551524580124562e2,
51*3e12c5d1SDavid du Colombier 0.10000000000000000000e1,
52*3e12c5d1SDavid du Colombier };
53*3e12c5d1SDavid du Colombier
54*3e12c5d1SDavid du Colombier static double
asym(double arg)55*3e12c5d1SDavid du Colombier asym(double arg)
56*3e12c5d1SDavid du Colombier {
57*3e12c5d1SDavid du Colombier double n, argsq;
58*3e12c5d1SDavid du Colombier int i;
59*3e12c5d1SDavid du Colombier
60*3e12c5d1SDavid du Colombier argsq = 1 / (arg*arg);
61*3e12c5d1SDavid du Colombier n = 0;
62*3e12c5d1SDavid du Colombier for(i=M-1; i>=0; i--)
63*3e12c5d1SDavid du Colombier n = n*argsq + p1[i];
64*3e12c5d1SDavid du Colombier return (arg-.5)*log(arg) - arg + goobie + n/arg;
65*3e12c5d1SDavid du Colombier }
66*3e12c5d1SDavid du Colombier
67*3e12c5d1SDavid du Colombier static double
pos(double arg)68*3e12c5d1SDavid du Colombier pos(double arg)
69*3e12c5d1SDavid du Colombier {
70*3e12c5d1SDavid du Colombier double n, d, s;
71*3e12c5d1SDavid du Colombier int i;
72*3e12c5d1SDavid du Colombier
73*3e12c5d1SDavid du Colombier if(arg < 2)
74*3e12c5d1SDavid du Colombier return pos(arg+1)/arg;
75*3e12c5d1SDavid du Colombier if(arg > 3)
76*3e12c5d1SDavid du Colombier return (arg-1)*pos(arg-1);
77*3e12c5d1SDavid du Colombier
78*3e12c5d1SDavid du Colombier s = arg - 2;
79*3e12c5d1SDavid du Colombier n = 0;
80*3e12c5d1SDavid du Colombier d = 0;
81*3e12c5d1SDavid du Colombier for(i=N-1; i>=0; i--){
82*3e12c5d1SDavid du Colombier n = n*s + p2[i];
83*3e12c5d1SDavid du Colombier d = d*s + q2[i];
84*3e12c5d1SDavid du Colombier }
85*3e12c5d1SDavid du Colombier return n/d;
86*3e12c5d1SDavid du Colombier }
87*3e12c5d1SDavid du Colombier
88*3e12c5d1SDavid du Colombier static double
neg(double arg)89*3e12c5d1SDavid du Colombier neg(double arg)
90*3e12c5d1SDavid du Colombier {
91*3e12c5d1SDavid du Colombier double temp;
92*3e12c5d1SDavid du Colombier
93*3e12c5d1SDavid du Colombier arg = -arg;
94*3e12c5d1SDavid du Colombier temp = sin(pi*arg);
95*3e12c5d1SDavid du Colombier if(temp == 0) {
96*3e12c5d1SDavid du Colombier errno = EDOM;
97*3e12c5d1SDavid du Colombier return HUGE_VAL;
98*3e12c5d1SDavid du Colombier }
99*3e12c5d1SDavid du Colombier if(temp < 0)
100*3e12c5d1SDavid du Colombier temp = -temp;
101*3e12c5d1SDavid du Colombier else
102*3e12c5d1SDavid du Colombier signgam = -1;
103*3e12c5d1SDavid du Colombier return -log(arg*pos(arg)*temp/pi);
104*3e12c5d1SDavid du Colombier }
105*3e12c5d1SDavid du Colombier
106*3e12c5d1SDavid du Colombier double
gamma(double arg)107*3e12c5d1SDavid du Colombier gamma(double arg)
108*3e12c5d1SDavid du Colombier {
109*3e12c5d1SDavid du Colombier
110*3e12c5d1SDavid du Colombier signgam = 1;
111*3e12c5d1SDavid du Colombier if(arg <= 0)
112*3e12c5d1SDavid du Colombier return neg(arg);
113*3e12c5d1SDavid du Colombier if(arg > 8)
114*3e12c5d1SDavid du Colombier return asym(arg);
115*3e12c5d1SDavid du Colombier return log(pos(arg));
116*3e12c5d1SDavid du Colombier }
117