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