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