1/* @(#)er_lgamma.c 5.1 93/09/24 */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/* __ieee754_lgamma_r(x, signgamp)
14 * Reentrant version of the logarithm of the Gamma function
15 * with user provide pointer for the sign of Gamma(x).
16 *
17 * Method:
18 * 1. Argument Reduction for 0 < x <= 8
19 * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
20 * reduce x to a number in [1.5,2.5] by
21 * lgamma(1+s) = log(s) + lgamma(s)
22 * for example,
23 * lgamma(7.3) = log(6.3) + lgamma(6.3)
24 * = log(6.3*5.3) + lgamma(5.3)
25 * = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
26 * 2. Polynomial approximation of lgamma around its
27 * minimun ymin=1.461632144968362245 to maintain monotonicity.
28 * On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
29 * Let z = x-ymin;
30 * lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
31 * where
32 * poly(z) is a 14 degree polynomial.
33 * 2. Rational approximation in the primary interval [2,3]
34 * We use the following approximation:
35 * s = x-2.0;
36 * lgamma(x) = 0.5*s + s*P(s)/Q(s)
37 * with accuracy
38 * |P/Q - (lgamma(x)-0.5s)| < 2**-61.71
39 * Our algorithms are based on the following observation
40 *
41 * zeta(2)-1 2 zeta(3)-1 3
42 * lgamma(2+s) = s*(1-Euler) + --------- * s - --------- * s + ...
43 * 2 3
44 *
45 * where Euler = 0.5771... is the Euler constant, which is very
46 * close to 0.5.
47 *
48 * 3. For x>=8, we have
49 * lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
50 * (better formula:
51 * lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
52 * Let z = 1/x, then we approximation
53 * f(z) = lgamma(x) - (x-0.5)(log(x)-1)
54 * by
55 * 3 5 11
56 * w = w0 + w1*z + w2*z + w3*z + ... + w6*z
57 * where
58 * |w - f(z)| < 2**-58.74
59 *
60 * 4. For negative x, since (G is gamma function)
61 * -x*G(-x)*G(x) = pi/sin(pi*x),
62 * we have
63 * G(x) = pi/(sin(pi*x)*(-x)*G(-x))
64 * since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
65 * Hence, for x<0, signgam = sign(sin(pi*x)) and
66 * lgamma(x) = log(|Gamma(x)|)
67 * = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
68 * Note: one should avoid compute pi*(-x) directly in the
69 * computation of sin(pi*(-x)).
70 *
71 * 5. Special Cases
72 * lgamma(2+s) ~ s*(1-Euler) for tiny s
73 * lgamma(1)=lgamma(2)=0
74 * lgamma(x) ~ -log(x) for tiny x
75 * lgamma(0) = lgamma(inf) = inf
76 * lgamma(-integer) = +-inf
77 *
78 */
79
80#include <math.h>
81#include <math_private.h>
82#include <libc-diag.h>
83
84static const double
85two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
86half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
87one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
88pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
89a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */
90a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */
91a2 = 6.73523010531292681824e-02, /* 0x3FB13E00, 0x1A5562A7 */
92a3 = 2.05808084325167332806e-02, /* 0x3F951322, 0xAC92547B */
93a4 = 7.38555086081402883957e-03, /* 0x3F7E404F, 0xB68FEFE8 */
94a5 = 2.89051383673415629091e-03, /* 0x3F67ADD8, 0xCCB7926B */
95a6 = 1.19270763183362067845e-03, /* 0x3F538A94, 0x116F3F5D */
96a7 = 5.10069792153511336608e-04, /* 0x3F40B6C6, 0x89B99C00 */
97a8 = 2.20862790713908385557e-04, /* 0x3F2CF2EC, 0xED10E54D */
98a9 = 1.08011567247583939954e-04, /* 0x3F1C5088, 0x987DFB07 */
99a10 = 2.52144565451257326939e-05, /* 0x3EFA7074, 0x428CFA52 */
100a11 = 4.48640949618915160150e-05, /* 0x3F07858E, 0x90A45837 */
101tc = 1.46163214496836224576e+00, /* 0x3FF762D8, 0x6356BE3F */
102tf = -1.21486290535849611461e-01, /* 0xBFBF19B9, 0xBCC38A42 */
103/* tt = -(tail of tf) */
104tt = -3.63867699703950536541e-18, /* 0xBC50C7CA, 0xA48A971F */
105t0 = 4.83836122723810047042e-01, /* 0x3FDEF72B, 0xC8EE38A2 */
106t1 = -1.47587722994593911752e-01, /* 0xBFC2E427, 0x8DC6C509 */
107t2 = 6.46249402391333854778e-02, /* 0x3FB08B42, 0x94D5419B */
108t3 = -3.27885410759859649565e-02, /* 0xBFA0C9A8, 0xDF35B713 */
109t4 = 1.79706750811820387126e-02, /* 0x3F9266E7, 0x970AF9EC */
110t5 = -1.03142241298341437450e-02, /* 0xBF851F9F, 0xBA91EC6A */
111t6 = 6.10053870246291332635e-03, /* 0x3F78FCE0, 0xE370E344 */
112t7 = -3.68452016781138256760e-03, /* 0xBF6E2EFF, 0xB3E914D7 */
113t8 = 2.25964780900612472250e-03, /* 0x3F6282D3, 0x2E15C915 */
114t9 = -1.40346469989232843813e-03, /* 0xBF56FE8E, 0xBF2D1AF1 */
115t10 = 8.81081882437654011382e-04, /* 0x3F4CDF0C, 0xEF61A8E9 */
116t11 = -5.38595305356740546715e-04, /* 0xBF41A610, 0x9C73E0EC */
117t12 = 3.15632070903625950361e-04, /* 0x3F34AF6D, 0x6C0EBBF7 */
118t13 = -3.12754168375120860518e-04, /* 0xBF347F24, 0xECC38C38 */
119t14 = 3.35529192635519073543e-04, /* 0x3F35FD3E, 0xE8C2D3F4 */
120u0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
121u1 = 6.32827064025093366517e-01, /* 0x3FE4401E, 0x8B005DFF */
122u2 = 1.45492250137234768737e+00, /* 0x3FF7475C, 0xD119BD6F */
123u3 = 9.77717527963372745603e-01, /* 0x3FEF4976, 0x44EA8450 */
124u4 = 2.28963728064692451092e-01, /* 0x3FCD4EAE, 0xF6010924 */
125u5 = 1.33810918536787660377e-02, /* 0x3F8B678B, 0xBF2BAB09 */
126v1 = 2.45597793713041134822e+00, /* 0x4003A5D7, 0xC2BD619C */
127v2 = 2.12848976379893395361e+00, /* 0x40010725, 0xA42B18F5 */
128v3 = 7.69285150456672783825e-01, /* 0x3FE89DFB, 0xE45050AF */
129v4 = 1.04222645593369134254e-01, /* 0x3FBAAE55, 0xD6537C88 */
130v5 = 3.21709242282423911810e-03, /* 0x3F6A5ABB, 0x57D0CF61 */
131s0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
132s1 = 2.14982415960608852501e-01, /* 0x3FCB848B, 0x36E20878 */
133s2 = 3.25778796408930981787e-01, /* 0x3FD4D98F, 0x4F139F59 */
134s3 = 1.46350472652464452805e-01, /* 0x3FC2BB9C, 0xBEE5F2F7 */
135s4 = 2.66422703033638609560e-02, /* 0x3F9B481C, 0x7E939961 */
136s5 = 1.84028451407337715652e-03, /* 0x3F5E26B6, 0x7368F239 */
137s6 = 3.19475326584100867617e-05, /* 0x3F00BFEC, 0xDD17E945 */
138r1 = 1.39200533467621045958e+00, /* 0x3FF645A7, 0x62C4AB74 */
139r2 = 7.21935547567138069525e-01, /* 0x3FE71A18, 0x93D3DCDC */
140r3 = 1.71933865632803078993e-01, /* 0x3FC601ED, 0xCCFBDF27 */
141r4 = 1.86459191715652901344e-02, /* 0x3F9317EA, 0x742ED475 */
142r5 = 7.77942496381893596434e-04, /* 0x3F497DDA, 0xCA41A95B */
143r6 = 7.32668430744625636189e-06, /* 0x3EDEBAF7, 0xA5B38140 */
144w0 = 4.18938533204672725052e-01, /* 0x3FDACFE3, 0x90C97D69 */
145w1 = 8.33333333333329678849e-02, /* 0x3FB55555, 0x5555553B */
146w2 = -2.77777777728775536470e-03, /* 0xBF66C16C, 0x16B02E5C */
147w3 = 7.93650558643019558500e-04, /* 0x3F4A019F, 0x98CF38B6 */
148w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
149w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
150w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
151
152static const double zero= 0.00000000000000000000e+00;
153
154static double
155sin_pi(double x)
156{
157 double y,z;
158 int n,ix;
159
160 GET_HIGH_WORD(ix,x);
161 ix &= 0x7fffffff;
162
163 if(ix<0x3fd00000) return __sin(pi*x);
164 y = -x; /* x is assume negative */
165
166 /*
167 * argument reduction, make sure inexact flag not raised if input
168 * is an integer
169 */
170 z = __floor(y);
171 if(z!=y) { /* inexact anyway */
172 y *= 0.5;
173 y = 2.0*(y - __floor(y)); /* y = |x| mod 2.0 */
174 n = (int) (y*4.0);
175 } else {
176 if(ix>=0x43400000) {
177 y = zero; n = 0; /* y must be even */
178 } else {
179 if(ix<0x43300000) z = y+two52; /* exact */
180 GET_LOW_WORD(n,z);
181 n &= 1;
182 y = n;
183 n<<= 2;
184 }
185 }
186 switch (n) {
187 case 0: y = __sin(pi*y); break;
188 case 1:
189 case 2: y = __cos(pi*(0.5-y)); break;
190 case 3:
191 case 4: y = __sin(pi*(one-y)); break;
192 case 5:
193 case 6: y = -__cos(pi*(y-1.5)); break;
194 default: y = __sin(pi*(y-2.0)); break;
195 }
196 return -y;
197}
198
199
200double
201__ieee754_lgamma_r(double x, int *signgamp)
202{
203 double t,y,z,nadj,p,p1,p2,p3,q,r,w;
204 int i,hx,lx,ix;
205
206 EXTRACT_WORDS(hx,lx,x);
207
208 /* purge off +-inf, NaN, +-0, and negative arguments */
209 *signgamp = 1;
210 ix = hx&0x7fffffff;
211 if(__builtin_expect(ix>=0x7ff00000, 0)) return x*x;
212 if(__builtin_expect((ix|lx)==0, 0))
213 {
214 if (hx < 0)
215 *signgamp = -1;
216 return one/fabs(x);
217 }
218 if(__builtin_expect(ix<0x3b900000, 0)) {
219 /* |x|<2**-70, return -log(|x|) */
220 if(hx<0) {
221 *signgamp = -1;
222 return -__ieee754_log(-x);
223 } else return -__ieee754_log(x);
224 }
225 if(hx<0) {
226 if(__builtin_expect(ix>=0x43300000, 0))
227 /* |x|>=2**52, must be -integer */
228 return fabs (x)/zero;
229 if (x < -2.0 && x > -28.0)
230 return __lgamma_neg (x, signgamp);
231 t = sin_pi(x);
232 if(t==zero) return one/fabsf(t); /* -integer */
233 nadj = __ieee754_log(pi/fabs(t*x));
234 if(t<zero) *signgamp = -1;
235 x = -x;
236 }
237
238 /* purge off 1 and 2 */
239 if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0;
240 /* for x < 2.0 */
241 else if(ix<0x40000000) {
242 if(ix<=0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */
243 r = -__ieee754_log(x);
244 if(ix>=0x3FE76944) {y = one-x; i= 0;}
245 else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;}
246 else {y = x; i=2;}
247 } else {
248 r = zero;
249 if(ix>=0x3FFBB4C3) {y=2.0-x;i=0;} /* [1.7316,2] */
250 else if(ix>=0x3FF3B4C4) {y=x-tc;i=1;} /* [1.23,1.73] */
251 else {y=x-one;i=2;}
252 }
253 switch(i) {
254 case 0:
255 z = y*y;
256 p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
257 p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
258 p = y*p1+p2;
259 r += (p-0.5*y); break;
260 case 1:
261 z = y*y;
262 w = z*y;
263 p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */
264 p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
265 p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
266 p = z*p1-(tt-w*(p2+y*p3));
267 r += (tf + p); break;
268 case 2:
269 p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
270 p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
271 r += (-0.5*y + p1/p2);
272 }
273 }
274 else if(ix<0x40200000) { /* x < 8.0 */
275 i = (int)x;
276 t = zero;
277 y = x-(double)i;
278 p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
279 q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
280 r = half*y+p/q;
281 z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
282 switch(i) {
283 case 7: z *= (y+6.0); /* FALLTHRU */
284 case 6: z *= (y+5.0); /* FALLTHRU */
285 case 5: z *= (y+4.0); /* FALLTHRU */
286 case 4: z *= (y+3.0); /* FALLTHRU */
287 case 3: z *= (y+2.0); /* FALLTHRU */
288 r += __ieee754_log(z); break;
289 }
290 /* 8.0 <= x < 2**58 */
291 } else if (ix < 0x43900000) {
292 t = __ieee754_log(x);
293 z = one/x;
294 y = z*z;
295 w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
296 r = (x-half)*(t-one)+w;
297 } else
298 /* 2**58 <= x <= inf */
299 r = math_narrow_eval (x*(__ieee754_log(x)-one));
300 /* NADJ is set for negative arguments but not otherwise,
301 resulting in warnings that it may be used uninitialized
302 although in the cases where it is used it has always been
303 set. */
304 DIAG_PUSH_NEEDS_COMMENT;
305 DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wmaybe-uninitialized");
306 if(hx<0) r = nadj - r;
307 DIAG_POP_NEEDS_COMMENT;
308 return r;
309}
310strong_alias (__ieee754_lgamma_r, __lgamma_r_finite)
311