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