148402Sbostic /*-
248402Sbostic  * Copyright (c) 1985 The Regents of the University of California.
348402Sbostic  * All rights reserved.
448402Sbostic  *
548402Sbostic  * %sccs.include.proprietary.c%
634120Sbostic  */
734120Sbostic 
824600Szliu #ifndef lint
9*56952Sbostic static char sccsid[] = "@(#)lgamma.c	5.5 (Berkeley) 12/02/92";
1034120Sbostic #endif /* not lint */
1124600Szliu 
12*56952Sbostic #include <math.h>
13*56952Sbostic #include <errno.h>
1424600Szliu 
15*56952Sbostic #include "mathimpl.h"
1624600Szliu 
17*56952Sbostic /* TRUNC sets trailing bits in a floating-point number to zero.
18*56952Sbostic  * ptrx points to the second half of the floating-point number.
19*56952Sbostic  * x1 is a temporary variable.
2024600Szliu */
2124600Szliu 
22*56952Sbostic #if defined(vax) || defined(tahoe)
23*56952Sbostic #define LSHFT_PLUS_1	(134217728.+1.)
24*56952Sbostic #define TRUNC(x, dummy, x1) x1 = x*LSHFT_PLUS_1, x -= x1, x += x1
25*56952Sbostic #else
26*56952Sbostic #define MASK 0xf8000000
27*56952Sbostic #define TRUNC(x, ptrx, dummy) *ptrx &= MASK
28*56952Sbostic #endif
2935679Sbostic 
30*56952Sbostic #define x0	0.461632144968362356785
31*56952Sbostic #define x0_lo	-.000000000000000015522348162858676890521
32*56952Sbostic #define a0_hi	-0.12148629053584961146
33*56952Sbostic #define a0_lo	3.07435645275902737e-18
34*56952Sbostic #define LEFT	(1.0 - (x0 + .25))
35*56952Sbostic #define RIGHT	(x0 - .218)
36*56952Sbostic #define lns2pi_hi 0.418945312500000
37*56952Sbostic #define lns2pi_lo -.000006779295327258219670263595
3824600Szliu 
39*56952Sbostic #define UNDERFL (1e-1020 * 1e-1020)
40*56952Sbostic int signgam;
41*56952Sbostic 
42*56952Sbostic #define r0	-0.02771227512955130520
43*56952Sbostic #define r1	-0.2980729795228150847
44*56952Sbostic #define r2	-0.3257411333183093394
45*56952Sbostic #define r3	-0.1126814387531706041
46*56952Sbostic #define r4	-0.01129130057170225562
47*56952Sbostic #define r5	-2.259650588213369095e-05
48*56952Sbostic #define s0	1.7144571600017144419
49*56952Sbostic #define s1	2.7864695046181946481
50*56952Sbostic #define s2	1.5645463655191798047
51*56952Sbostic #define s3	0.34858463899811098496
52*56952Sbostic #define s4	0.024677593453636563481
53*56952Sbostic 
54*56952Sbostic #define p0     -7.721566490153286087127140e-02
55*56952Sbostic #define p1	2.077324848654884653970785e-01
56*56952Sbostic #define p2	3.474331160945523535787194e-01
57*56952Sbostic #define p3	1.724375677840324429295524e-01
58*56952Sbostic #define p4	3.546181984297784658205969e-02
59*56952Sbostic #define p5	2.866163630124151557532506e-03
60*56952Sbostic #define p6	6.143168512963655570532770e-05
61*56952Sbostic #define q0	1.000000000000000000000000e+00
62*56952Sbostic #define q1	1.485897307300750567469226e+00
63*56952Sbostic #define q2	8.336064915387758261556045e-01
64*56952Sbostic #define q3	2.185469782512070977585556e-01
65*56952Sbostic #define q4	2.671060746990048983840354e-02
66*56952Sbostic #define q5	1.296631961305647947057711e-03
67*56952Sbostic #define q6	1.484566091079246905938151e-05
68*56952Sbostic 
69*56952Sbostic #define NP2 8
70*56952Sbostic double P2[] = {
71*56952Sbostic .0833333333333333148296162562474,
72*56952Sbostic -.00277777777774548123579378966497,
73*56952Sbostic .000793650778754435631476282786423,
74*56952Sbostic -.000595235082566672847950717262222,
75*56952Sbostic .000841428560346653702135821806252,
76*56952Sbostic -.00189773526463879200348872089421,
77*56952Sbostic .00569394463439411649408050664078,
78*56952Sbostic -.0144705562421428915453880392761
7924600Szliu };
8024600Szliu 
81*56952Sbostic static double neg_lgam __P((double));
82*56952Sbostic static double small_lgam __P((double));
83*56952Sbostic static double large_lgam __P((double));
8435679Sbostic 
8524600Szliu double
86*56952Sbostic lgamma(x)
87*56952Sbostic 	double x;
8824600Szliu {
89*56952Sbostic 	double zero = 0.0, one = 1.0;
90*56952Sbostic 	double r;
91*56952Sbostic 	signgam = 1;
92*56952Sbostic 	if (!finite(x)) {
93*56952Sbostic 		errno = EDOM;
94*56952Sbostic 		if (x < 0)
95*56952Sbostic 			x = -x;
96*56952Sbostic 		else if (x > 0)
97*56952Sbostic 			errno = ERANGE;
98*56952Sbostic 		return (x);
99*56952Sbostic 	}
100*56952Sbostic 	if (x > 6 + RIGHT) {
101*56952Sbostic 		if (x > 1.0e20)
102*56952Sbostic 			return (x*(log(x)-one));
103*56952Sbostic 		r = large_lgam(x);
104*56952Sbostic 		return (r);
105*56952Sbostic 	} else if (x > 1e-17)
106*56952Sbostic 		return (small_lgam(x));
107*56952Sbostic 	else if(x > -1e-17) {
108*56952Sbostic 		if (x < 0)
109*56952Sbostic 			signgam = -1, x = -x;
110*56952Sbostic 		return (-log(x));
111*56952Sbostic 	} else
112*56952Sbostic 		return (neg_lgam(x));
11324600Szliu }
11424600Szliu 
115*56952Sbostic /* Accurate to max(ulp(1/128) absolute, 2^-75 relative) error. */
11624600Szliu static double
117*56952Sbostic large_lgam(x)
118*56952Sbostic 	double x;
11924600Szliu {
120*56952Sbostic 	double z, p, x1;
12124600Szliu 	int i;
122*56952Sbostic 	long *pva, *pua;
123*56952Sbostic 	struct Double t, u, v;
124*56952Sbostic 	pua = (long *) &u.a, pua++;
125*56952Sbostic 	pva = (long *) &v.a, pva++;
126*56952Sbostic 	z = 1.0/(x*x);
127*56952Sbostic 	for (p = P2[i = NP2-1]; --i >= 0;)
128*56952Sbostic 		p = P2[i] + p*z;	/* error in approximation = 2.8e-18 */
12924600Szliu 
130*56952Sbostic 	p = p/x;			/* ulp = 1.7e-18; error < 1.6ulp */
131*56952Sbostic 					/* 0 < frac < 1/64 (at x = 5.5) */
132*56952Sbostic 	t = log__D(x);
133*56952Sbostic 	t.a -= 1.0;
134*56952Sbostic 	u.a = t.a + t.b;
135*56952Sbostic 	TRUNC (u.a, pua, x1);		/* truncate u.a */
136*56952Sbostic 	u.b = (t.a - u.a);
137*56952Sbostic 	u.b += t.b;
138*56952Sbostic 	x -= .5;
139*56952Sbostic 	v.a = x;
140*56952Sbostic 	TRUNC(v.a, pva, x1);		/* truncate v.a */
141*56952Sbostic 	v.b = x - v.a;
142*56952Sbostic 	t.a = v.a*u.a;			/* t = (x-.5)*(log(x)-1) */
143*56952Sbostic 	t.b = v.b*u.a + x*u.b;
144*56952Sbostic 	z = t.b + lns2pi_lo;		/* return t + lns2pi + p */
145*56952Sbostic 	z += p; z += lns2pi_hi;
146*56952Sbostic 	z += t.a;
147*56952Sbostic 	return (z);
14824600Szliu }
149*56952Sbostic /* Good to < 1 ulp.  (provably .90 ulp; .87 ulp on 1,000,000 runs.)
150*56952Sbostic    It also has correct monotonicity.
151*56952Sbostic  */
15224600Szliu static double
153*56952Sbostic small_lgam(x)
154*56952Sbostic 	double x;
15524600Szliu {
156*56952Sbostic 	int xi;
157*56952Sbostic 	double y, z, t, r = 0, p, q;
158*56952Sbostic 	struct Double rr;
15924600Szliu 
160*56952Sbostic 	/* Do nasty area near the minimum.  No good for 1.25 < x < 2.5 */
161*56952Sbostic 	if (x < 2.0 - LEFT && x > 1.0 + RIGHT) {
162*56952Sbostic 		t = x - 1.0; t -= x0;
163*56952Sbostic 		z = t - x0_lo;
164*56952Sbostic 		p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5))));
165*56952Sbostic 		q = s0+z*(s1+z*(s2+z*(s3+z*s4)));
166*56952Sbostic 		r = t*(z*(p/q) - x0_lo) + a0_lo;
167*56952Sbostic 		r += .5*t*t; r += a0_hi;
168*56952Sbostic 		return (r);
16924600Szliu 	}
170*56952Sbostic 	xi = (x + .5);
171*56952Sbostic 	y = x - xi;
172*56952Sbostic 	rr.a = rr.b = 0;
173*56952Sbostic 	if (y < -LEFT) {	/* necessary for 2.5 < x < 2.72.. */
174*56952Sbostic 		t = y + (1.0 - x0);
175*56952Sbostic 		z = t - x0_lo;
176*56952Sbostic 		p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5))));
177*56952Sbostic 		q = s0+z*(s1+z*(s2+z*(s3+z*s4)));
178*56952Sbostic 		r = t*(z*(p/q) - x0_lo) + a0_lo;
179*56952Sbostic 		r += .5*t*t;
180*56952Sbostic 		q = a0_hi;
181*56952Sbostic 		printf("(0)q = %.18lg r = %.18lg\n", q, r);
182*56952Sbostic 	} else {
183*56952Sbostic 		p = y*(p0+y*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*p6))))));
184*56952Sbostic 		q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*q6)))));
185*56952Sbostic 		r = p/q;
186*56952Sbostic 		q = .5*y;
187*56952Sbostic 		printf("(1)q = %.18lg r = %.18lg\n", q, r);
18824600Szliu 	}
189*56952Sbostic 	printf("y = %lg, r = %.18lg\n", y, r);
190*56952Sbostic 	z = 1.0;
191*56952Sbostic 	switch (xi) {
192*56952Sbostic 	case 6:	z  = (y + 5);
193*56952Sbostic 	case 5:	z *= (y + 4);
194*56952Sbostic 	case 4:	z *= (y + 3);
195*56952Sbostic 	case 3:	z *= (y + 2);
196*56952Sbostic 		rr = log__D(z);
197*56952Sbostic 		printf("%.21lg, %lg\n", rr.a, rr.b);
198*56952Sbostic 		r += rr.b; q += r;
199*56952Sbostic 		return(q + rr.a);
200*56952Sbostic 	case 2:	return (q + r);
201*56952Sbostic 
202*56952Sbostic 	case 0: printf("r = %lg\n", r);
203*56952Sbostic 		r -= log1p(x); printf("\t%lg\n", r);
204*56952Sbostic 	default: rr = log__D(x);
205*56952Sbostic 		r -= rr.b; printf("\t%lg\n", r);
206*56952Sbostic 	}
207*56952Sbostic 	if (q > .5 *rr.a) {
208*56952Sbostic 		printf("q = %lg, rr = %lg, r = %lg\n", q, rr.a, r);
209*56952Sbostic 		q -= rr.a;
210*56952Sbostic 		return(r + q);
211*56952Sbostic 	} else
212*56952Sbostic 		printf("q = %lg, rr = %lg, r = %lg\n", q, rr.a, r);
213*56952Sbostic 		return((r + q) - rr.a);
21424600Szliu }
21524600Szliu 
216*56952Sbostic #define lpi_hi 1.1447298858494001638
217*56952Sbostic #define lpi_lo .0000000000000000102659511627078262
218*56952Sbostic /* Error: within 3.5 ulp for x < 171.  For large x, see lgamma. */
21924600Szliu static double
220*56952Sbostic neg_lgam(x)
221*56952Sbostic 	double x;
22224600Szliu {
223*56952Sbostic 	struct Double lg, lsine;
224*56952Sbostic 	double y, z, one = 1.0, zero = 0.0;
22524600Szliu 
226*56952Sbostic 	z = floor(x + .5);
227*56952Sbostic 	if (z == x) {
228*56952Sbostic 		errno = EDOM;
229*56952Sbostic 		return (one/zero);	/* convention: G(-(integer)) -> +oo */
230*56952Sbostic 	}
231*56952Sbostic 	y = ceil(x);
232*56952Sbostic 	if (y*.5 == ceil(.5*y))
233*56952Sbostic 		signgam = -1;
23424600Szliu 
235*56952Sbostic 	x = -x;
236*56952Sbostic 	z = fabs(x + z);	/* 0 < z <= .5 */
237*56952Sbostic 	if (z < .25)
238*56952Sbostic 		z = sin(M_PI*z);
239*56952Sbostic 	else
240*56952Sbostic 		z = cos(M_PI*(.5-z));
241*56952Sbostic 	z = -log(z*x/M_PI);
242*56952Sbostic 
243*56952Sbostic 	if (x > 6.5)
244*56952Sbostic 		y -= large_lgam(x);
245*56952Sbostic 	else
246*56952Sbostic 		y = -small_lgam (x);
247*56952Sbostic 	return (y + z);
24824600Szliu }
249