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