1 #ifndef lint
2 static char sccsid[] = "@(#)lgamma.c 4.4 (Berkeley) 09/11/85";
3 #endif not lint
4
5 /*
6 C program for floating point log Gamma function
7
8 lgamma(x) computes the log of the absolute
9 value of the Gamma function.
10 The sign of the Gamma function is returned in the
11 external quantity signgam.
12
13 The coefficients for expansion around zero
14 are #5243 from Hart & Cheney; for expansion
15 around infinity they are #5404.
16
17 Calls log, floor and sin.
18 */
19
20 #include <math.h>
21 #ifdef VAX
22 #include <errno.h>
23 #endif
24 int signgam = 0;
25 static double goobie = 0.9189385332046727417803297; /* log(2*pi)/2 */
26 static double pi = 3.1415926535897932384626434;
27
28 #define M 6
29 #define N 8
30 static double p1[] = {
31 0.83333333333333101837e-1,
32 -.277777777735865004e-2,
33 0.793650576493454e-3,
34 -.5951896861197e-3,
35 0.83645878922e-3,
36 -.1633436431e-2,
37 };
38 static double p2[] = {
39 -.42353689509744089647e5,
40 -.20886861789269887364e5,
41 -.87627102978521489560e4,
42 -.20085274013072791214e4,
43 -.43933044406002567613e3,
44 -.50108693752970953015e2,
45 -.67449507245925289918e1,
46 0.0,
47 };
48 static double q2[] = {
49 -.42353689509744090010e5,
50 -.29803853309256649932e4,
51 0.99403074150827709015e4,
52 -.15286072737795220248e4,
53 -.49902852662143904834e3,
54 0.18949823415702801641e3,
55 -.23081551524580124562e2,
56 0.10000000000000000000e1,
57 };
58
59 double
lgamma(arg)60 lgamma(arg)
61 double arg;
62 {
63 double log(), pos(), neg(), asym();
64
65 signgam = 1.;
66 if(arg <= 0.) return(neg(arg));
67 if(arg > 8.) return(asym(arg));
68 return(log(pos(arg)));
69 }
70
71 static double
asym(arg)72 asym(arg)
73 double arg;
74 {
75 double log();
76 double n, argsq;
77 int i;
78
79 argsq = 1./(arg*arg);
80 for(n=0,i=M-1; i>=0; i--){
81 n = n*argsq + p1[i];
82 }
83 return((arg-.5)*log(arg) - arg + goobie + n/arg);
84 }
85
86 static double
neg(arg)87 neg(arg)
88 double arg;
89 {
90 double t;
91 double log(), sin(), floor(), pos();
92
93 arg = -arg;
94 /*
95 * to see if arg were a true integer, the old code used the
96 * mathematically correct observation:
97 * sin(n*pi) = 0 <=> n is an integer.
98 * but in finite precision arithmetic, sin(n*PI) will NEVER
99 * be zero simply because n*PI is a rational number. hence
100 * it failed to work with our newer, more accurate sin()
101 * which uses true pi to do the argument reduction...
102 * temp = sin(pi*arg);
103 */
104 t = floor(arg);
105 if (arg - t > 0.5e0)
106 t += 1.e0; /* t := integer nearest arg */
107 #ifdef VAX
108 if (arg == t) {
109 extern double infnan();
110 return(infnan(ERANGE)); /* +INF */
111 }
112 #endif
113 signgam = (int) (t - 2*floor(t/2)); /* signgam = 1 if t was odd, */
114 /* 0 if t was even */
115 signgam = signgam - 1 + signgam; /* signgam = 1 if t was odd, */
116 /* -1 if t was even */
117 t = arg - t; /* -0.5 <= t <= 0.5 */
118 if (t < 0.e0) {
119 t = -t;
120 signgam = -signgam;
121 }
122 return(-log(arg*pos(arg)*sin(pi*t)/pi));
123 }
124
125 static double
pos(arg)126 pos(arg)
127 double arg;
128 {
129 double n, d, s;
130 register i;
131
132 if(arg < 2.) return(pos(arg+1.)/arg);
133 if(arg > 3.) return((arg-1.)*pos(arg-1.));
134
135 s = arg - 2.;
136 for(n=0,d=0,i=N-1; i>=0; i--){
137 n = n*s + p2[i];
138 d = d*s + q2[i];
139 }
140 return(n/d);
141 }
142