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