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