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-narrow-eval.h>
82 #include <math_private.h>
83 #include <libc-diag.h>
84 #include <libm-alias-finite.h>
85
86 static const double
87 two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
88 half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
89 one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
90 pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
91 a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */
92 a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */
93 a2 = 6.73523010531292681824e-02, /* 0x3FB13E00, 0x1A5562A7 */
94 a3 = 2.05808084325167332806e-02, /* 0x3F951322, 0xAC92547B */
95 a4 = 7.38555086081402883957e-03, /* 0x3F7E404F, 0xB68FEFE8 */
96 a5 = 2.89051383673415629091e-03, /* 0x3F67ADD8, 0xCCB7926B */
97 a6 = 1.19270763183362067845e-03, /* 0x3F538A94, 0x116F3F5D */
98 a7 = 5.10069792153511336608e-04, /* 0x3F40B6C6, 0x89B99C00 */
99 a8 = 2.20862790713908385557e-04, /* 0x3F2CF2EC, 0xED10E54D */
100 a9 = 1.08011567247583939954e-04, /* 0x3F1C5088, 0x987DFB07 */
101 a10 = 2.52144565451257326939e-05, /* 0x3EFA7074, 0x428CFA52 */
102 a11 = 4.48640949618915160150e-05, /* 0x3F07858E, 0x90A45837 */
103 tc = 1.46163214496836224576e+00, /* 0x3FF762D8, 0x6356BE3F */
104 tf = -1.21486290535849611461e-01, /* 0xBFBF19B9, 0xBCC38A42 */
105 /* tt = -(tail of tf) */
106 tt = -3.63867699703950536541e-18, /* 0xBC50C7CA, 0xA48A971F */
107 t0 = 4.83836122723810047042e-01, /* 0x3FDEF72B, 0xC8EE38A2 */
108 t1 = -1.47587722994593911752e-01, /* 0xBFC2E427, 0x8DC6C509 */
109 t2 = 6.46249402391333854778e-02, /* 0x3FB08B42, 0x94D5419B */
110 t3 = -3.27885410759859649565e-02, /* 0xBFA0C9A8, 0xDF35B713 */
111 t4 = 1.79706750811820387126e-02, /* 0x3F9266E7, 0x970AF9EC */
112 t5 = -1.03142241298341437450e-02, /* 0xBF851F9F, 0xBA91EC6A */
113 t6 = 6.10053870246291332635e-03, /* 0x3F78FCE0, 0xE370E344 */
114 t7 = -3.68452016781138256760e-03, /* 0xBF6E2EFF, 0xB3E914D7 */
115 t8 = 2.25964780900612472250e-03, /* 0x3F6282D3, 0x2E15C915 */
116 t9 = -1.40346469989232843813e-03, /* 0xBF56FE8E, 0xBF2D1AF1 */
117 t10 = 8.81081882437654011382e-04, /* 0x3F4CDF0C, 0xEF61A8E9 */
118 t11 = -5.38595305356740546715e-04, /* 0xBF41A610, 0x9C73E0EC */
119 t12 = 3.15632070903625950361e-04, /* 0x3F34AF6D, 0x6C0EBBF7 */
120 t13 = -3.12754168375120860518e-04, /* 0xBF347F24, 0xECC38C38 */
121 t14 = 3.35529192635519073543e-04, /* 0x3F35FD3E, 0xE8C2D3F4 */
122 u0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
123 u1 = 6.32827064025093366517e-01, /* 0x3FE4401E, 0x8B005DFF */
124 u2 = 1.45492250137234768737e+00, /* 0x3FF7475C, 0xD119BD6F */
125 u3 = 9.77717527963372745603e-01, /* 0x3FEF4976, 0x44EA8450 */
126 u4 = 2.28963728064692451092e-01, /* 0x3FCD4EAE, 0xF6010924 */
127 u5 = 1.33810918536787660377e-02, /* 0x3F8B678B, 0xBF2BAB09 */
128 v1 = 2.45597793713041134822e+00, /* 0x4003A5D7, 0xC2BD619C */
129 v2 = 2.12848976379893395361e+00, /* 0x40010725, 0xA42B18F5 */
130 v3 = 7.69285150456672783825e-01, /* 0x3FE89DFB, 0xE45050AF */
131 v4 = 1.04222645593369134254e-01, /* 0x3FBAAE55, 0xD6537C88 */
132 v5 = 3.21709242282423911810e-03, /* 0x3F6A5ABB, 0x57D0CF61 */
133 s0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */
134 s1 = 2.14982415960608852501e-01, /* 0x3FCB848B, 0x36E20878 */
135 s2 = 3.25778796408930981787e-01, /* 0x3FD4D98F, 0x4F139F59 */
136 s3 = 1.46350472652464452805e-01, /* 0x3FC2BB9C, 0xBEE5F2F7 */
137 s4 = 2.66422703033638609560e-02, /* 0x3F9B481C, 0x7E939961 */
138 s5 = 1.84028451407337715652e-03, /* 0x3F5E26B6, 0x7368F239 */
139 s6 = 3.19475326584100867617e-05, /* 0x3F00BFEC, 0xDD17E945 */
140 r1 = 1.39200533467621045958e+00, /* 0x3FF645A7, 0x62C4AB74 */
141 r2 = 7.21935547567138069525e-01, /* 0x3FE71A18, 0x93D3DCDC */
142 r3 = 1.71933865632803078993e-01, /* 0x3FC601ED, 0xCCFBDF27 */
143 r4 = 1.86459191715652901344e-02, /* 0x3F9317EA, 0x742ED475 */
144 r5 = 7.77942496381893596434e-04, /* 0x3F497DDA, 0xCA41A95B */
145 r6 = 7.32668430744625636189e-06, /* 0x3EDEBAF7, 0xA5B38140 */
146 w0 = 4.18938533204672725052e-01, /* 0x3FDACFE3, 0x90C97D69 */
147 w1 = 8.33333333333329678849e-02, /* 0x3FB55555, 0x5555553B */
148 w2 = -2.77777777728775536470e-03, /* 0xBF66C16C, 0x16B02E5C */
149 w3 = 7.93650558643019558500e-04, /* 0x3F4A019F, 0x98CF38B6 */
150 w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
151 w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
152 w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
153
154 static const double zero= 0.00000000000000000000e+00;
155
156 static double
sin_pi(double x)157 sin_pi(double x)
158 {
159 double y,z;
160 int n,ix;
161
162 GET_HIGH_WORD(ix,x);
163 ix &= 0x7fffffff;
164
165 if(ix<0x3fd00000) return __sin(pi*x);
166 y = -x; /* x is assume negative */
167
168 /*
169 * argument reduction, make sure inexact flag not raised if input
170 * is an integer
171 */
172 z = floor(y);
173 if(z!=y) { /* inexact anyway */
174 y *= 0.5;
175 y = 2.0*(y - floor(y)); /* y = |x| mod 2.0 */
176 n = (int) (y*4.0);
177 } else {
178 if(ix>=0x43400000) {
179 y = zero; n = 0; /* y must be even */
180 } else {
181 if(ix<0x43300000) z = y+two52; /* exact */
182 GET_LOW_WORD(n,z);
183 n &= 1;
184 y = n;
185 n<<= 2;
186 }
187 }
188 switch (n) {
189 case 0: y = __sin(pi*y); break;
190 case 1:
191 case 2: y = __cos(pi*(0.5-y)); break;
192 case 3:
193 case 4: y = __sin(pi*(one-y)); break;
194 case 5:
195 case 6: y = -__cos(pi*(y-1.5)); break;
196 default: y = __sin(pi*(y-2.0)); break;
197 }
198 return -y;
199 }
200
201
202 double
__ieee754_lgamma_r(double x,int * signgamp)203 __ieee754_lgamma_r(double x, int *signgamp)
204 {
205 double t,y,z,nadj,p,p1,p2,p3,q,r,w;
206 int i,hx,lx,ix;
207
208 EXTRACT_WORDS(hx,lx,x);
209
210 /* purge off +-inf, NaN, +-0, and negative arguments */
211 *signgamp = 1;
212 ix = hx&0x7fffffff;
213 if(__builtin_expect(ix>=0x7ff00000, 0)) return x*x;
214 if(__builtin_expect((ix|lx)==0, 0))
215 {
216 if (hx < 0)
217 *signgamp = -1;
218 return one/fabs(x);
219 }
220 if(__builtin_expect(ix<0x3b900000, 0)) {
221 /* |x|<2**-70, return -log(|x|) */
222 if(hx<0) {
223 *signgamp = -1;
224 return -__ieee754_log(-x);
225 } else return -__ieee754_log(x);
226 }
227 if(hx<0) {
228 if(__builtin_expect(ix>=0x43300000, 0))
229 /* |x|>=2**52, must be -integer */
230 return fabs (x)/zero;
231 if (x < -2.0 && x > -28.0)
232 return __lgamma_neg (x, signgamp);
233 t = sin_pi(x);
234 if(t==zero) return one/fabsf(t); /* -integer */
235 nadj = __ieee754_log(pi/fabs(t*x));
236 if(t<zero) *signgamp = -1;
237 x = -x;
238 }
239
240 /* purge off 1 and 2 */
241 if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0;
242 /* for x < 2.0 */
243 else if(ix<0x40000000) {
244 if(ix<=0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */
245 r = -__ieee754_log(x);
246 if(ix>=0x3FE76944) {y = one-x; i= 0;}
247 else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;}
248 else {y = x; i=2;}
249 } else {
250 r = zero;
251 if(ix>=0x3FFBB4C3) {y=2.0-x;i=0;} /* [1.7316,2] */
252 else if(ix>=0x3FF3B4C4) {y=x-tc;i=1;} /* [1.23,1.73] */
253 else {y=x-one;i=2;}
254 }
255 switch(i) {
256 case 0:
257 z = y*y;
258 p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
259 p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
260 p = y*p1+p2;
261 r += (p-0.5*y); break;
262 case 1:
263 z = y*y;
264 w = z*y;
265 p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */
266 p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
267 p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
268 p = z*p1-(tt-w*(p2+y*p3));
269 r += (tf + p); break;
270 case 2:
271 p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
272 p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
273 r += (-0.5*y + p1/p2);
274 }
275 }
276 else if(ix<0x40200000) { /* x < 8.0 */
277 i = (int)x;
278 t = zero;
279 y = x-(double)i;
280 p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
281 q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
282 r = half*y+p/q;
283 z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
284 switch(i) {
285 case 7: z *= (y+6.0); /* FALLTHRU */
286 case 6: z *= (y+5.0); /* FALLTHRU */
287 case 5: z *= (y+4.0); /* FALLTHRU */
288 case 4: z *= (y+3.0); /* FALLTHRU */
289 case 3: z *= (y+2.0); /* FALLTHRU */
290 r += __ieee754_log(z); break;
291 }
292 /* 8.0 <= x < 2**58 */
293 } else if (ix < 0x43900000) {
294 t = __ieee754_log(x);
295 z = one/x;
296 y = z*z;
297 w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
298 r = (x-half)*(t-one)+w;
299 } else
300 /* 2**58 <= x <= inf */
301 r = math_narrow_eval (x*(__ieee754_log(x)-one));
302 /* NADJ is set for negative arguments but not otherwise,
303 resulting in warnings that it may be used uninitialized
304 although in the cases where it is used it has always been
305 set. */
306 DIAG_PUSH_NEEDS_COMMENT;
307 DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wmaybe-uninitialized");
308 if(hx<0) r = nadj - r;
309 DIAG_POP_NEEDS_COMMENT;
310 return r;
311 }
312 libm_alias_finite (__ieee754_lgamma_r, __lgamma_r)
313