1 /* $NetBSD: n_lgamma.c,v 1.8 2024/06/09 14:09:27 riastradh Exp $ */
2 /*-
3 * Copyright (c) 1992, 1993
4 * The Regents of the University of California. All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 * 2. Redistributions in binary form must reproduce the above copyright
12 * notice, this list of conditions and the following disclaimer in the
13 * documentation and/or other materials provided with the distribution.
14 * 3. Neither the name of the University nor the names of its contributors
15 * may be used to endorse or promote products derived from this software
16 * without specific prior written permission.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
19 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
22 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
24 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
25 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
27 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
28 * SUCH DAMAGE.
29 */
30
31 #ifndef lint
32 #if 0
33 static char sccsid[] = "@(#)lgamma.c 8.2 (Berkeley) 11/30/93";
34 #endif
35 #endif /* not lint */
36
37 /*
38 * Coded by Peter McIlroy, Nov 1992;
39 *
40 * The financial support of UUNET Communications Services is gratefully
41 * acknowledged.
42 */
43
44 #include <math.h>
45 #include <errno.h>
46
47 #include "mathimpl.h"
48
49 /* Log gamma function.
50 * Error: x > 0 error < 1.3ulp.
51 * x > 4, error < 1ulp.
52 * x > 9, error < .6ulp.
53 * x < 0, all bets are off. (When G(x) ~ 1, log(G(x)) ~ 0)
54 * Method:
55 * x > 6:
56 * Use the asymptotic expansion (Stirling's Formula)
57 * 0 < x < 6:
58 * Use gamma(x+1) = x*gamma(x) for argument reduction.
59 * Use rational approximation in
60 * the range 1.2, 2.5
61 * Two approximations are used, one centered at the
62 * minimum to ensure monotonicity; one centered at 2
63 * to maintain small relative error.
64 * x < 0:
65 * Use the reflection formula,
66 * G(1-x)G(x) = PI/sin(PI*x)
67 * Special values:
68 * non-positive integer returns +Inf.
69 * NaN returns NaN
70 */
71 #if defined(__vax__) || defined(tahoe)
72 #define _IEEE 0
73 /* double and float have same size exponent field */
74 #define TRUNC(x) x = (double) (float) (x)
75 #else
76 static int endian;
77 #define _IEEE 1
78 #define TRUNC(x) *(((int *) &x) + endian) &= 0xf8000000
79 #define infnan(x) 0.0
80 #endif
81
82 static double small_lgam(double);
83 static double large_lgam(double);
84 static double neg_lgam(double, int *);
85 static const double one = 1.0;
86 int signgam;
87
88 #define UNDERFL (1e-1020 * 1e-1020)
89
90 #define LEFT (1.0 - (x0 + .25))
91 #define RIGHT (x0 - .218)
92 /*
93 * Constants for approximation in [1.244,1.712]
94 */
95 #define x0 0.461632144968362356785
96 #define x0_lo -.000000000000000015522348162858676890521
97 #define a0_hi -0.12148629128932952880859
98 #define a0_lo .0000000007534799204229502
99 #define r0 -2.771227512955130520e-002
100 #define r1 -2.980729795228150847e-001
101 #define r2 -3.257411333183093394e-001
102 #define r3 -1.126814387531706041e-001
103 #define r4 -1.129130057170225562e-002
104 #define r5 -2.259650588213369095e-005
105 #define s0 1.714457160001714442e+000
106 #define s1 2.786469504618194648e+000
107 #define s2 1.564546365519179805e+000
108 #define s3 3.485846389981109850e-001
109 #define s4 2.467759345363656348e-002
110 /*
111 * Constants for approximation in [1.71, 2.5]
112 */
113 #define a1_hi 4.227843350984671344505727574870e-01
114 #define a1_lo 4.670126436531227189e-18
115 #define p0 3.224670334241133695662995251041e-01
116 #define p1 3.569659696950364669021382724168e-01
117 #define p2 1.342918716072560025853732668111e-01
118 #define p3 1.950702176409779831089963408886e-02
119 #define p4 8.546740251667538090796227834289e-04
120 #define q0 1.000000000000000444089209850062e+00
121 #define q1 1.315850076960161985084596381057e+00
122 #define q2 6.274644311862156431658377186977e-01
123 #define q3 1.304706631926259297049597307705e-01
124 #define q4 1.102815279606722369265536798366e-02
125 #define q5 2.512690594856678929537585620579e-04
126 #define q6 -1.003597548112371003358107325598e-06
127 /*
128 * Stirling's Formula, adjusted for equal-ripple. x in [6,Inf].
129 */
130 #define lns2pi .418938533204672741780329736405
131 #define pb0 8.33333333333333148296162562474e-02
132 #define pb1 -2.77777777774548123579378966497e-03
133 #define pb2 7.93650778754435631476282786423e-04
134 #define pb3 -5.95235082566672847950717262222e-04
135 #define pb4 8.41428560346653702135821806252e-04
136 #define pb5 -1.89773526463879200348872089421e-03
137 #define pb6 5.69394463439411649408050664078e-03
138 #define pb7 -1.44705562421428915453880392761e-02
139
__weak_alias(lgammal,lgamma)140 __weak_alias(lgammal, lgamma)
141 __weak_alias(lgammal_r, lgamma_r)
142
143 double
144 lgamma(double x)
145 {
146
147 return lgamma_r(x, &signgam);
148 }
149
150 double
lgamma_r(double x,int * signgamp)151 lgamma_r(double x, int *signgamp)
152 {
153 double r;
154
155 *signgamp = 1;
156 #if _IEEE
157 endian = ((*(int *) &one)) ? 1 : 0;
158 #endif
159
160 if (!finite(x)) {
161 if (_IEEE)
162 return (x+x);
163 else return (infnan(EDOM));
164 }
165
166 if (x > 6 + RIGHT) {
167 r = large_lgam(x);
168 return (r);
169 } else if (x > 1e-16)
170 return (small_lgam(x));
171 else if (x > -1e-16) {
172 if (x < 0)
173 *signgamp = -1, x = -x;
174 return (-log(x));
175 } else
176 return (neg_lgam(x, signgamp));
177 }
178
179 static double
large_lgam(double x)180 large_lgam(double x)
181 {
182 double z, p, x1;
183 struct Double t, u, v;
184 u = __log__D(x);
185 u.a -= 1.0;
186 if (x > 1e15) {
187 v.a = x - 0.5;
188 TRUNC(v.a);
189 v.b = (x - v.a) - 0.5;
190 t.a = u.a*v.a;
191 t.b = x*u.b + v.b*u.a;
192 if (_IEEE == 0 && !finite(t.a))
193 return(infnan(ERANGE));
194 return(t.a + t.b);
195 }
196 x1 = 1./x;
197 z = x1*x1;
198 p = pb0+z*(pb1+z*(pb2+z*(pb3+z*(pb4+z*(pb5+z*(pb6+z*pb7))))));
199 /* error in approximation = 2.8e-19 */
200
201 p = p*x1; /* error < 2.3e-18 absolute */
202 /* 0 < p < 1/64 (at x = 5.5) */
203 v.a = x = x - 0.5;
204 TRUNC(v.a); /* truncate v.a to 26 bits. */
205 v.b = x - v.a;
206 t.a = v.a*u.a; /* t = (x-.5)*(log(x)-1) */
207 t.b = v.b*u.a + x*u.b;
208 t.b += p; t.b += lns2pi; /* return t + lns2pi + p */
209 return (t.a + t.b);
210 }
211
212 static double
small_lgam(double x)213 small_lgam(double x)
214 {
215 int x_int;
216 double y, z, t, r = 0, p, q, hi, lo;
217 struct Double rr;
218 x_int = (x + .5);
219 y = x - x_int;
220 if (x_int <= 2 && y > RIGHT) {
221 t = y - x0;
222 y--; x_int++;
223 goto CONTINUE;
224 } else if (y < -LEFT) {
225 t = y +(1.0-x0);
226 CONTINUE:
227 z = t - x0_lo;
228 p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5))));
229 q = s0+z*(s1+z*(s2+z*(s3+z*s4)));
230 r = t*(z*(p/q) - x0_lo);
231 t = .5*t*t;
232 z = 1.0;
233 switch (x_int) {
234 case 6: z = (y + 5); /* FALLTHROUGH */
235 case 5: z *= (y + 4); /* FALLTHROUGH */
236 case 4: z *= (y + 3); /* FALLTHROUGH */
237 case 3: z *= (y + 2);
238 rr = __log__D(z);
239 rr.b += a0_lo; rr.a += a0_hi;
240 return(((r+rr.b)+t+rr.a));
241 case 2: return(((r+a0_lo)+t)+a0_hi);
242 case 0: r -= log1p(x); /* FALLTHROUGH */
243 default: rr = __log__D(x);
244 rr.a -= a0_hi; rr.b -= a0_lo;
245 return(((r - rr.b) + t) - rr.a);
246 }
247 } else {
248 p = p0+y*(p1+y*(p2+y*(p3+y*p4)));
249 q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*q6)))));
250 p = p*(y/q);
251 t = (double)(float) y;
252 z = y-t;
253 hi = (double)(float) (p+a1_hi);
254 lo = a1_hi - hi; lo += p; lo += a1_lo;
255 r = lo*y + z*hi; /* q + r = y*(a0+p/q) */
256 q = hi*t;
257 z = 1.0;
258 switch (x_int) {
259 case 6: z = (y + 5); /* FALLTHROUGH */
260 case 5: z *= (y + 4); /* FALLTHROUGH */
261 case 4: z *= (y + 3); /* FALLTHROUGH */
262 case 3: z *= (y + 2);
263 rr = __log__D(z);
264 r += rr.b; r += q;
265 return(rr.a + r);
266 case 2: return (q+ r);
267 case 0: rr = __log__D(x);
268 r -= rr.b; r -= log1p(x);
269 r += q; r-= rr.a;
270 return(r);
271 default: rr = __log__D(x);
272 r -= rr.b;
273 q -= rr.a;
274 return (r+q);
275 }
276 }
277 }
278
279 static double
neg_lgam(double x,int * signgamp)280 neg_lgam(double x, int *signgamp)
281 {
282 int xi;
283 double y, z, zero = 0.0;
284
285 /* avoid destructive cancellation as much as possible */
286 if (x > -170) {
287 xi = x;
288 if (xi == x) {
289 if (_IEEE)
290 return(one/zero);
291 else
292 return(infnan(ERANGE));
293 }
294 y = gamma(x);
295 if (y < 0)
296 y = -y, *signgamp = -1;
297 return (log(y));
298 }
299 z = floor(x + .5);
300 if (z == x) { /* convention: G(-(integer)) -> +Inf */
301 if (_IEEE)
302 return (one/zero);
303 else
304 return (infnan(ERANGE));
305 }
306 y = .5*ceil(x);
307 if (y == ceil(y))
308 *signgamp = -1;
309 x = -x;
310 z = fabs(x + z); /* 0 < z <= .5 */
311 if (z < .25)
312 z = sin(M_PI*z);
313 else
314 z = cos(M_PI*(0.5-z));
315 z = log(M_PI/(z*x));
316 y = large_lgam(x);
317 return (z - y);
318 }
319