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.3 (Berkeley) 09/22/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 "mathimpl.h"
27 #if defined(vax)||defined(tahoe)
28 #include <errno.h>
29 #endif	/* defined(vax)||defined(tahoe) */
30 
31 int	signgam = 0;
32 static const double goobie = 0.9189385332046727417803297;  /* log(2*pi)/2 */
33 static const double pi	   = 3.1415926535897932384626434;
34 
35 #define M 6
36 #define N 8
37 static const double p1[] = {
38 	0.83333333333333101837e-1,
39 	-.277777777735865004e-2,
40 	0.793650576493454e-3,
41 	-.5951896861197e-3,
42 	0.83645878922e-3,
43 	-.1633436431e-2,
44 };
45 static const double p2[] = {
46 	-.42353689509744089647e5,
47 	-.20886861789269887364e5,
48 	-.87627102978521489560e4,
49 	-.20085274013072791214e4,
50 	-.43933044406002567613e3,
51 	-.50108693752970953015e2,
52 	-.67449507245925289918e1,
53 	0.0,
54 };
55 static const double q2[] = {
56 	-.42353689509744090010e5,
57 	-.29803853309256649932e4,
58 	0.99403074150827709015e4,
59 	-.15286072737795220248e4,
60 	-.49902852662143904834e3,
61 	0.18949823415702801641e3,
62 	-.23081551524580124562e2,
63 	0.10000000000000000000e1,
64 };
65 
66 static double pos(), neg(), asym();
67 
68 double
69 lgamma(arg)
70 double arg;
71 {
72 
73 	signgam = 1.;
74 	if(arg <= 0.) return(neg(arg));
75 	if(arg > 8.) return(asym(arg));
76 	return(log(pos(arg)));
77 }
78 
79 static double
80 asym(arg)
81 double arg;
82 {
83 	double n, argsq;
84 	int i;
85 
86 	argsq = 1./(arg*arg);
87 	for(n=0,i=M-1; i>=0; i--){
88 		n = n*argsq + p1[i];
89 	}
90 	return((arg-.5)*log(arg) - arg + goobie + n/arg);
91 }
92 
93 static double
94 neg(arg)
95 double arg;
96 {
97 	double t;
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 	    return(infnan(ERANGE));		/* +INF */
116 	}
117 #endif	/* defined(vax)||defined(tahoe) */
118 	signgam = (int) (t - 2*floor(t/2));	/* signgam =  1 if t was odd, */
119 						/*            0 if t was even */
120 	signgam = signgam - 1 + signgam;	/* signgam =  1 if t was odd, */
121 						/*           -1 if t was even */
122 	t = arg - t;				/*  -0.5 <= t <= 0.5 */
123 	if (t < 0.e0) {
124 	    t = -t;
125 	    signgam = -signgam;
126 	}
127 	return(-log(arg*pos(arg)*sin(pi*t)/pi));
128 }
129 
130 static double
131 pos(arg)
132 double arg;
133 {
134 	double n, d, s;
135 	register i;
136 
137 	if(arg < 2.) return(pos(arg+1.)/arg);
138 	if(arg > 3.) return((arg-1.)*pos(arg-1.));
139 
140 	s = arg - 2.;
141 	for(n=0,d=0,i=N-1; i>=0; i--){
142 		n = n*s + p2[i];
143 		d = d*s + q2[i];
144 	}
145 	return(n/d);
146 }
147