1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_lgammal.c */
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
10 * ====================================================
13 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
15 * Permission to use, copy, modify, and distribute this software for any
16 * purpose with or without fee is hereby granted, provided that the above
17 * copyright notice and this permission notice appear in all copies.
19 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
20 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
21 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
22 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
23 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
24 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
25 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
28 * Reentrant version of the logarithm of the Gamma function
29 * with user provide pointer for the sign of Gamma(x).
32 * 1. Argument Reduction for 0 < x <= 8
33 * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
34 * reduce x to a number in [1.5,2.5] by
35 * lgamma(1+s) = log(s) + lgamma(s)
37 * lgamma(7.3) = log(6.3) + lgamma(6.3)
38 * = log(6.3*5.3) + lgamma(5.3)
39 * = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
40 * 2. Polynomial approximation of lgamma around its
41 * minimun ymin=1.461632144968362245 to maintain monotonicity.
42 * On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
44 * lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
45 * 2. Rational approximation in the primary interval [2,3]
46 * We use the following approximation:
48 * lgamma(x) = 0.5*s + s*P(s)/Q(s)
49 * Our algorithms are based on the following observation
51 * zeta(2)-1 2 zeta(3)-1 3
52 * lgamma(2+s) = s*(1-Euler) + --------- * s - --------- * s + ...
55 * where Euler = 0.5771... is the Euler constant, which is very
58 * 3. For x>=8, we have
59 * lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
61 * lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
62 * Let z = 1/x, then we approximation
63 * f(z) = lgamma(x) - (x-0.5)(log(x)-1)
66 * w = w0 + w1*z + w2*z + w3*z + ... + w6*z
68 * 4. For negative x, since (G is gamma function)
69 * -x*G(-x)*G(x) = pi/sin(pi*x),
71 * G(x) = pi/(sin(pi*x)*(-x)*G(-x))
72 * since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
73 * Hence, for x<0, signgam = sign(sin(pi*x)) and
74 * lgamma(x) = log(|Gamma(x)|)
75 * = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
76 * Note: one should avoid compute pi*(-x) directly in the
77 * computation of sin(pi*(-x)).
80 * lgamma(2+s) ~ s*(1-Euler) for tiny s
81 * lgamma(1)=lgamma(2)=0
82 * lgamma(x) ~ -log(x) for tiny x
83 * lgamma(0) = lgamma(inf) = inf
84 * lgamma(-integer) = +-inf
92 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
93 double __lgamma_r(double x
, int *sg
);
95 long double __lgammal_r(long double x
, int *sg
)
97 return __lgamma_r(x
, sg
);
99 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
100 static const long double
101 pi
= 3.14159265358979323846264L,
103 /* lgam(1+x) = 0.5 x + x a(x)/b(x)
104 -0.268402099609375 <= x <= 0
105 peak relative error 6.6e-22 */
106 a0
= -6.343246574721079391729402781192128239938E2L
,
107 a1
= 1.856560238672465796768677717168371401378E3L
,
108 a2
= 2.404733102163746263689288466865843408429E3L
,
109 a3
= 8.804188795790383497379532868917517596322E2L
,
110 a4
= 1.135361354097447729740103745999661157426E2L
,
111 a5
= 3.766956539107615557608581581190400021285E0L
,
113 b0
= 8.214973713960928795704317259806842490498E3L
,
114 b1
= 1.026343508841367384879065363925870888012E4L
,
115 b2
= 4.553337477045763320522762343132210919277E3L
,
116 b3
= 8.506975785032585797446253359230031874803E2L
,
117 b4
= 6.042447899703295436820744186992189445813E1L
,
118 /* b5 = 1.000000000000000000000000000000000000000E0 */
121 tc
= 1.4616321449683623412626595423257213284682E0L
,
122 tf
= -1.2148629053584961146050602565082954242826E-1, /* double precision */
123 /* tt = (tail of tf), i.e. tf + tt has extended precision. */
124 tt
= 3.3649914684731379602768989080467587736363E-18L,
125 /* lgam ( 1.4616321449683623412626595423257213284682E0 ) =
126 -1.2148629053584960809551455717769158215135617312999903886372437313313530E-1 */
128 /* lgam (x + tc) = tf + tt + x g(x)/h(x)
129 -0.230003726999612341262659542325721328468 <= x
130 <= 0.2699962730003876587373404576742786715318
131 peak relative error 2.1e-21 */
132 g0
= 3.645529916721223331888305293534095553827E-18L,
133 g1
= 5.126654642791082497002594216163574795690E3L
,
134 g2
= 8.828603575854624811911631336122070070327E3L
,
135 g3
= 5.464186426932117031234820886525701595203E3L
,
136 g4
= 1.455427403530884193180776558102868592293E3L
,
137 g5
= 1.541735456969245924860307497029155838446E2L
,
138 g6
= 4.335498275274822298341872707453445815118E0L
,
140 h0
= 1.059584930106085509696730443974495979641E4L
,
141 h1
= 2.147921653490043010629481226937850618860E4L
,
142 h2
= 1.643014770044524804175197151958100656728E4L
,
143 h3
= 5.869021995186925517228323497501767586078E3L
,
144 h4
= 9.764244777714344488787381271643502742293E2L
,
145 h5
= 6.442485441570592541741092969581997002349E1L
,
146 /* h6 = 1.000000000000000000000000000000000000000E0 */
149 /* lgam (x+1) = -0.5 x + x u(x)/v(x)
150 -0.100006103515625 <= x <= 0.231639862060546875
151 peak relative error 1.3e-21 */
152 u0
= -8.886217500092090678492242071879342025627E1L
,
153 u1
= 6.840109978129177639438792958320783599310E2L
,
154 u2
= 2.042626104514127267855588786511809932433E3L
,
155 u3
= 1.911723903442667422201651063009856064275E3L
,
156 u4
= 7.447065275665887457628865263491667767695E2L
,
157 u5
= 1.132256494121790736268471016493103952637E2L
,
158 u6
= 4.484398885516614191003094714505960972894E0L
,
160 v0
= 1.150830924194461522996462401210374632929E3L
,
161 v1
= 3.399692260848747447377972081399737098610E3L
,
162 v2
= 3.786631705644460255229513563657226008015E3L
,
163 v3
= 1.966450123004478374557778781564114347876E3L
,
164 v4
= 4.741359068914069299837355438370682773122E2L
,
165 v5
= 4.508989649747184050907206782117647852364E1L
,
166 /* v6 = 1.000000000000000000000000000000000000000E0 */
169 /* lgam (x+2) = .5 x + x s(x)/r(x)
171 peak relative error 7.2e-22 */
172 s0
= 1.454726263410661942989109455292824853344E6L
,
173 s1
= -3.901428390086348447890408306153378922752E6L
,
174 s2
= -6.573568698209374121847873064292963089438E6L
,
175 s3
= -3.319055881485044417245964508099095984643E6L
,
176 s4
= -7.094891568758439227560184618114707107977E5L
,
177 s5
= -6.263426646464505837422314539808112478303E4L
,
178 s6
= -1.684926520999477529949915657519454051529E3L
,
180 r0
= -1.883978160734303518163008696712983134698E7L
,
181 r1
= -2.815206082812062064902202753264922306830E7L
,
182 r2
= -1.600245495251915899081846093343626358398E7L
,
183 r3
= -4.310526301881305003489257052083370058799E6L
,
184 r4
= -5.563807682263923279438235987186184968542E5L
,
185 r5
= -3.027734654434169996032905158145259713083E4L
,
186 r6
= -4.501995652861105629217250715790764371267E2L
,
187 /* r6 = 1.000000000000000000000000000000000000000E0 */
190 /* lgam(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x w(1/x^2)
192 Peak relative error 1.51e-21
194 w0
= 4.189385332046727417803e-1L,
195 w1
= 8.333333333333331447505E-2L,
196 w2
= -2.777777777750349603440E-3L,
197 w3
= 7.936507795855070755671E-4L,
198 w4
= -5.952345851765688514613E-4L,
199 w5
= 8.412723297322498080632E-4L,
200 w6
= -1.880801938119376907179E-3L,
201 w7
= 4.885026142432270781165E-3L;
203 /* sin(pi*x) assuming x > 2^-1000, if sin(pi*x)==0 the sign is arbitrary */
204 static long double sin_pi(long double x
)
208 /* spurious inexact if odd int */
210 x
= 2.0*(x
- floorl(x
)); /* x mod 2.0 */
218 default: /* case 4: */
219 case 0: return __sinl(x
, 0.0, 0);
220 case 1: return __cosl(x
, 0.0);
221 case 2: return __sinl(-x
, 0.0, 0);
222 case 3: return -__cosl(x
, 0.0);
226 long double __lgammal_r(long double x
, int *sg
) {
227 long double t
, y
, z
, nadj
, p
, p1
, p2
, q
, r
, w
;
228 union ldshape u
= {x
};
229 uint32_t ix
= (u
.i
.se
& 0x7fffU
)<<16 | u
.i
.m
>>48;
230 int sign
= u
.i
.se
>> 15;
235 /* purge off +-inf, NaN, +-0, tiny and negative arguments */
236 if (ix
>= 0x7fff0000)
238 if (ix
< 0x3fc08000) { /* |x|<2**-63, return -log(|x|) */
249 return 1.0 / (x
-x
); /* -integer */
254 nadj
= logl(pi
/ (t
* x
));
257 /* purge off 1 and 2 (so the sign is ok with downward rounding) */
258 if ((ix
== 0x3fff8000 || ix
== 0x40008000) && u
.i
.m
== 0) {
260 } else if (ix
< 0x40008000) { /* x < 2.0 */
261 if (ix
<= 0x3ffee666) { /* 8.99993896484375e-1 */
262 /* lgamma(x) = lgamma(x+1) - log(x) */
264 if (ix
>= 0x3ffebb4a) { /* 7.31597900390625e-1 */
267 } else if (ix
>= 0x3ffced33) { /* 2.31639862060546875e-1 */
270 } else { /* x < 0.23 */
276 if (ix
>= 0x3fffdda6) { /* 1.73162841796875 */
280 } else if (ix
>= 0x3fff9da6) { /* 1.23162841796875 */
292 p1
= a0
+ y
* (a1
+ y
* (a2
+ y
* (a3
+ y
* (a4
+ y
* a5
))));
293 p2
= b0
+ y
* (b1
+ y
* (b2
+ y
* (b3
+ y
* (b4
+ y
))));
294 r
+= 0.5 * y
+ y
* p1
/p2
;
297 p1
= g0
+ y
* (g1
+ y
* (g2
+ y
* (g3
+ y
* (g4
+ y
* (g5
+ y
* g6
)))));
298 p2
= h0
+ y
* (h1
+ y
* (h2
+ y
* (h3
+ y
* (h4
+ y
* (h5
+ y
)))));
303 p1
= y
* (u0
+ y
* (u1
+ y
* (u2
+ y
* (u3
+ y
* (u4
+ y
* (u5
+ y
* u6
))))));
304 p2
= v0
+ y
* (v1
+ y
* (v2
+ y
* (v3
+ y
* (v4
+ y
* (v5
+ y
)))));
305 r
+= (-0.5 * y
+ p1
/ p2
);
307 } else if (ix
< 0x40028000) { /* 8.0 */
311 p
= y
* (s0
+ y
* (s1
+ y
* (s2
+ y
* (s3
+ y
* (s4
+ y
* (s5
+ y
* s6
))))));
312 q
= r0
+ y
* (r1
+ y
* (r2
+ y
* (r3
+ y
* (r4
+ y
* (r5
+ y
* (r6
+ y
))))));
315 /* lgamma(1+s) = log(s) + lgamma(s) */
318 z
*= (y
+ 6.0); /* FALLTHRU */
320 z
*= (y
+ 5.0); /* FALLTHRU */
322 z
*= (y
+ 4.0); /* FALLTHRU */
324 z
*= (y
+ 3.0); /* FALLTHRU */
326 z
*= (y
+ 2.0); /* FALLTHRU */
330 } else if (ix
< 0x40418000) { /* 2^66 */
331 /* 8.0 <= x < 2**66 */
335 w
= w0
+ z
* (w1
+ y
* (w2
+ y
* (w3
+ y
* (w4
+ y
* (w5
+ y
* (w6
+ y
* w7
))))));
336 r
= (x
- 0.5) * (t
- 1.0) + w
;
337 } else /* 2**66 <= x <= inf */
338 r
= x
* (logl(x
) - 1.0);
343 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
344 // TODO: broken implementation to make things compile
345 double __lgamma_r(double x
, int *sg
);
347 long double __lgammal_r(long double x
, int *sg
)
349 return __lgamma_r(x
, sg
);
353 extern int __signgam
;
355 long double lgammal(long double x
)
357 return __lgammal_r(x
, &__signgam
);
360 weak_alias(__lgammal_r
, lgammal_r
);