1 /*-
2 * Copyright (c) 1992, 1993
3 * The Regents of the University of California. All rights reserved.
4 *
5 * %sccs.include.redist.c%
6 */
7
8 #ifndef lint
9 static char sccsid[] = "@(#)lgamma.c 8.2 (Berkeley) 11/30/93";
10 #endif /* not lint */
11
12 /*
13 * Coded by Peter McIlroy, Nov 1992;
14 *
15 * The financial support of UUNET Communications Services is greatfully
16 * acknowledged.
17 */
18
19 #include <math.h>
20 #include <errno.h>
21
22 #include "mathimpl.h"
23
24 /* Log gamma function.
25 * Error: x > 0 error < 1.3ulp.
26 * x > 4, error < 1ulp.
27 * x > 9, error < .6ulp.
28 * x < 0, all bets are off. (When G(x) ~ 1, log(G(x)) ~ 0)
29 * Method:
30 * x > 6:
31 * Use the asymptotic expansion (Stirling's Formula)
32 * 0 < x < 6:
33 * Use gamma(x+1) = x*gamma(x) for argument reduction.
34 * Use rational approximation in
35 * the range 1.2, 2.5
36 * Two approximations are used, one centered at the
37 * minimum to ensure monotonicity; one centered at 2
38 * to maintain small relative error.
39 * x < 0:
40 * Use the reflection formula,
41 * G(1-x)G(x) = PI/sin(PI*x)
42 * Special values:
43 * non-positive integer returns +Inf.
44 * NaN returns NaN
45 */
46 static int endian;
47 #if defined(vax) || defined(tahoe)
48 #define _IEEE 0
49 /* double and float have same size exponent field */
50 #define TRUNC(x) x = (double) (float) (x)
51 #else
52 #define _IEEE 1
53 #define TRUNC(x) *(((int *) &x) + endian) &= 0xf8000000
54 #define infnan(x) 0.0
55 #endif
56
57 static double small_lgam(double);
58 static double large_lgam(double);
59 static double neg_lgam(double);
60 static double zero = 0.0, one = 1.0;
61 int signgam;
62
63 #define UNDERFL (1e-1020 * 1e-1020)
64
65 #define LEFT (1.0 - (x0 + .25))
66 #define RIGHT (x0 - .218)
67 /*
68 /* Constants for approximation in [1.244,1.712]
69 */
70 #define x0 0.461632144968362356785
71 #define x0_lo -.000000000000000015522348162858676890521
72 #define a0_hi -0.12148629128932952880859
73 #define a0_lo .0000000007534799204229502
74 #define r0 -2.771227512955130520e-002
75 #define r1 -2.980729795228150847e-001
76 #define r2 -3.257411333183093394e-001
77 #define r3 -1.126814387531706041e-001
78 #define r4 -1.129130057170225562e-002
79 #define r5 -2.259650588213369095e-005
80 #define s0 1.714457160001714442e+000
81 #define s1 2.786469504618194648e+000
82 #define s2 1.564546365519179805e+000
83 #define s3 3.485846389981109850e-001
84 #define s4 2.467759345363656348e-002
85 /*
86 * Constants for approximation in [1.71, 2.5]
87 */
88 #define a1_hi 4.227843350984671344505727574870e-01
89 #define a1_lo 4.670126436531227189e-18
90 #define p0 3.224670334241133695662995251041e-01
91 #define p1 3.569659696950364669021382724168e-01
92 #define p2 1.342918716072560025853732668111e-01
93 #define p3 1.950702176409779831089963408886e-02
94 #define p4 8.546740251667538090796227834289e-04
95 #define q0 1.000000000000000444089209850062e+00
96 #define q1 1.315850076960161985084596381057e+00
97 #define q2 6.274644311862156431658377186977e-01
98 #define q3 1.304706631926259297049597307705e-01
99 #define q4 1.102815279606722369265536798366e-02
100 #define q5 2.512690594856678929537585620579e-04
101 #define q6 -1.003597548112371003358107325598e-06
102 /*
103 * Stirling's Formula, adjusted for equal-ripple. x in [6,Inf].
104 */
105 #define lns2pi .418938533204672741780329736405
106 #define pb0 8.33333333333333148296162562474e-02
107 #define pb1 -2.77777777774548123579378966497e-03
108 #define pb2 7.93650778754435631476282786423e-04
109 #define pb3 -5.95235082566672847950717262222e-04
110 #define pb4 8.41428560346653702135821806252e-04
111 #define pb5 -1.89773526463879200348872089421e-03
112 #define pb6 5.69394463439411649408050664078e-03
113 #define pb7 -1.44705562421428915453880392761e-02
114
115 __pure double
lgamma(double x)116 lgamma(double x)
117 {
118 double r;
119
120 signgam = 1;
121 endian = ((*(int *) &one)) ? 1 : 0;
122
123 if (!finite(x))
124 if (_IEEE)
125 return (x+x);
126 else return (infnan(EDOM));
127
128 if (x > 6 + RIGHT) {
129 r = large_lgam(x);
130 return (r);
131 } else if (x > 1e-16)
132 return (small_lgam(x));
133 else if (x > -1e-16) {
134 if (x < 0)
135 signgam = -1, x = -x;
136 return (-log(x));
137 } else
138 return (neg_lgam(x));
139 }
140
141 static double
large_lgam(double x)142 large_lgam(double x)
143 {
144 double z, p, x1;
145 int i;
146 struct Double t, u, v;
147 u = __log__D(x);
148 u.a -= 1.0;
149 if (x > 1e15) {
150 v.a = x - 0.5;
151 TRUNC(v.a);
152 v.b = (x - v.a) - 0.5;
153 t.a = u.a*v.a;
154 t.b = x*u.b + v.b*u.a;
155 if (_IEEE == 0 && !finite(t.a))
156 return(infnan(ERANGE));
157 return(t.a + t.b);
158 }
159 x1 = 1./x;
160 z = x1*x1;
161 p = pb0+z*(pb1+z*(pb2+z*(pb3+z*(pb4+z*(pb5+z*(pb6+z*pb7))))));
162 /* error in approximation = 2.8e-19 */
163
164 p = p*x1; /* error < 2.3e-18 absolute */
165 /* 0 < p < 1/64 (at x = 5.5) */
166 v.a = x = x - 0.5;
167 TRUNC(v.a); /* truncate v.a to 26 bits. */
168 v.b = x - v.a;
169 t.a = v.a*u.a; /* t = (x-.5)*(log(x)-1) */
170 t.b = v.b*u.a + x*u.b;
171 t.b += p; t.b += lns2pi; /* return t + lns2pi + p */
172 return (t.a + t.b);
173 }
174
175 static double
small_lgam(double x)176 small_lgam(double x)
177 {
178 int x_int;
179 double y, z, t, r = 0, p, q, hi, lo;
180 struct Double rr;
181 x_int = (x + .5);
182 y = x - x_int;
183 if (x_int <= 2 && y > RIGHT) {
184 t = y - x0;
185 y--; x_int++;
186 goto CONTINUE;
187 } else if (y < -LEFT) {
188 t = y +(1.0-x0);
189 CONTINUE:
190 z = t - x0_lo;
191 p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5))));
192 q = s0+z*(s1+z*(s2+z*(s3+z*s4)));
193 r = t*(z*(p/q) - x0_lo);
194 t = .5*t*t;
195 z = 1.0;
196 switch (x_int) {
197 case 6: z = (y + 5);
198 case 5: z *= (y + 4);
199 case 4: z *= (y + 3);
200 case 3: z *= (y + 2);
201 rr = __log__D(z);
202 rr.b += a0_lo; rr.a += a0_hi;
203 return(((r+rr.b)+t+rr.a));
204 case 2: return(((r+a0_lo)+t)+a0_hi);
205 case 0: r -= log1p(x);
206 default: rr = __log__D(x);
207 rr.a -= a0_hi; rr.b -= a0_lo;
208 return(((r - rr.b) + t) - rr.a);
209 }
210 } else {
211 p = p0+y*(p1+y*(p2+y*(p3+y*p4)));
212 q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*q6)))));
213 p = p*(y/q);
214 t = (double)(float) y;
215 z = y-t;
216 hi = (double)(float) (p+a1_hi);
217 lo = a1_hi - hi; lo += p; lo += a1_lo;
218 r = lo*y + z*hi; /* q + r = y*(a0+p/q) */
219 q = hi*t;
220 z = 1.0;
221 switch (x_int) {
222 case 6: z = (y + 5);
223 case 5: z *= (y + 4);
224 case 4: z *= (y + 3);
225 case 3: z *= (y + 2);
226 rr = __log__D(z);
227 r += rr.b; r += q;
228 return(rr.a + r);
229 case 2: return (q+ r);
230 case 0: rr = __log__D(x);
231 r -= rr.b; r -= log1p(x);
232 r += q; r-= rr.a;
233 return(r);
234 default: rr = __log__D(x);
235 r -= rr.b;
236 q -= rr.a;
237 return (r+q);
238 }
239 }
240 }
241
242 static double
neg_lgam(double x)243 neg_lgam(double x)
244 {
245 int xi;
246 double y, z, one = 1.0, zero = 0.0;
247 extern double gamma();
248
249 /* avoid destructive cancellation as much as possible */
250 if (x > -170) {
251 xi = x;
252 if (xi == x)
253 if (_IEEE)
254 return(one/zero);
255 else
256 return(infnan(ERANGE));
257 y = gamma(x);
258 if (y < 0)
259 y = -y, signgam = -1;
260 return (log(y));
261 }
262 z = floor(x + .5);
263 if (z == x) { /* convention: G(-(integer)) -> +Inf */
264 if (_IEEE)
265 return (one/zero);
266 else
267 return (infnan(ERANGE));
268 }
269 y = .5*ceil(x);
270 if (y == ceil(y))
271 signgam = -1;
272 x = -x;
273 z = fabs(x + z); /* 0 < z <= .5 */
274 if (z < .25)
275 z = sin(M_PI*z);
276 else
277 z = cos(M_PI*(0.5-z));
278 z = log(M_PI/(z*x));
279 y = large_lgam(x);
280 return (z - y);
281 }
282