2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
9 * ====================================================
12 /* Long double expansions are
13 Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
14 and are incorporated herein by permission of the author. The author
15 reserves the right to distribute this material elsewhere under different
16 copying permissions. These modifications are distributed here under
19 This library is free software; you can redistribute it and/or
20 modify it under the terms of the GNU Lesser General Public
21 License as published by the Free Software Foundation; either
22 version 2.1 of the License, or (at your option) any later version.
24 This library is distributed in the hope that it will be useful,
25 but WITHOUT ANY WARRANTY; without even the implied warranty of
26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27 Lesser General Public License for more details.
29 You should have received a copy of the GNU Lesser General Public
30 License along with this library; if not, write to the Free Software
31 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
33 /* __ieee754_lgammal_r(x, signgamp)
34 * Reentrant version of the logarithm of the Gamma function
35 * with user provide pointer for the sign of Gamma(x).
38 * 1. Argument Reduction for 0 < x <= 8
39 * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
40 * reduce x to a number in [1.5,2.5] by
41 * lgamma(1+s) = log(s) + lgamma(s)
43 * lgamma(7.3) = log(6.3) + lgamma(6.3)
44 * = log(6.3*5.3) + lgamma(5.3)
45 * = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
46 * 2. Polynomial approximation of lgamma around its
47 * minimun ymin=1.461632144968362245 to maintain monotonicity.
48 * On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
50 * lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
51 * 2. Rational approximation in the primary interval [2,3]
52 * We use the following approximation:
54 * lgamma(x) = 0.5*s + s*P(s)/Q(s)
55 * Our algorithms are based on the following observation
57 * zeta(2)-1 2 zeta(3)-1 3
58 * lgamma(2+s) = s*(1-Euler) + --------- * s - --------- * s + ...
61 * where Euler = 0.5771... is the Euler constant, which is very
64 * 3. For x>=8, we have
65 * lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
67 * lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
68 * Let z = 1/x, then we approximation
69 * f(z) = lgamma(x) - (x-0.5)(log(x)-1)
72 * w = w0 + w1*z + w2*z + w3*z + ... + w6*z
74 * 4. For negative x, since (G is gamma function)
75 * -x*G(-x)*G(x) = pi/sin(pi*x),
77 * G(x) = pi/(sin(pi*x)*(-x)*G(-x))
78 * since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
79 * Hence, for x<0, signgam = sign(sin(pi*x)) and
80 * lgamma(x) = log(|Gamma(x)|)
81 * = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
82 * Note: one should avoid compute pi*(-x) directly in the
83 * computation of sin(pi*(-x)).
86 * lgamma(2+s) ~ s*(1-Euler) for tiny s
87 * lgamma(1)=lgamma(2)=0
88 * lgamma(x) ~ -log(x) for tiny x
89 * lgamma(0) = lgamma(inf) = inf
90 * lgamma(-integer) = +-inf
95 #include "math_private.h"
98 static const long double
104 pi
= 3.14159265358979323846264L,
105 two63
= 9.223372036854775808e18L
,
107 /* lgam(1+x) = 0.5 x + x a(x)/b(x)
108 -0.268402099609375 <= x <= 0
109 peak relative error 6.6e-22 */
110 a0
= -6.343246574721079391729402781192128239938E2L
,
111 a1
= 1.856560238672465796768677717168371401378E3L
,
112 a2
= 2.404733102163746263689288466865843408429E3L
,
113 a3
= 8.804188795790383497379532868917517596322E2L
,
114 a4
= 1.135361354097447729740103745999661157426E2L
,
115 a5
= 3.766956539107615557608581581190400021285E0L
,
117 b0
= 8.214973713960928795704317259806842490498E3L
,
118 b1
= 1.026343508841367384879065363925870888012E4L
,
119 b2
= 4.553337477045763320522762343132210919277E3L
,
120 b3
= 8.506975785032585797446253359230031874803E2L
,
121 b4
= 6.042447899703295436820744186992189445813E1L
,
122 /* b5 = 1.000000000000000000000000000000000000000E0 */
125 tc
= 1.4616321449683623412626595423257213284682E0L
,
126 tf
= -1.2148629053584961146050602565082954242826E-1,/* double precision */
127 /* tt = (tail of tf), i.e. tf + tt has extended precision. */
128 tt
= 3.3649914684731379602768989080467587736363E-18L,
129 /* lgam ( 1.4616321449683623412626595423257213284682E0 ) =
130 -1.2148629053584960809551455717769158215135617312999903886372437313313530E-1 */
132 /* lgam (x + tc) = tf + tt + x g(x)/h(x)
133 - 0.230003726999612341262659542325721328468 <= x
134 <= 0.2699962730003876587373404576742786715318
135 peak relative error 2.1e-21 */
136 g0
= 3.645529916721223331888305293534095553827E-18L,
137 g1
= 5.126654642791082497002594216163574795690E3L
,
138 g2
= 8.828603575854624811911631336122070070327E3L
,
139 g3
= 5.464186426932117031234820886525701595203E3L
,
140 g4
= 1.455427403530884193180776558102868592293E3L
,
141 g5
= 1.541735456969245924860307497029155838446E2L
,
142 g6
= 4.335498275274822298341872707453445815118E0L
,
144 h0
= 1.059584930106085509696730443974495979641E4L
,
145 h1
= 2.147921653490043010629481226937850618860E4L
,
146 h2
= 1.643014770044524804175197151958100656728E4L
,
147 h3
= 5.869021995186925517228323497501767586078E3L
,
148 h4
= 9.764244777714344488787381271643502742293E2L
,
149 h5
= 6.442485441570592541741092969581997002349E1L
,
150 /* h6 = 1.000000000000000000000000000000000000000E0 */
153 /* lgam (x+1) = -0.5 x + x u(x)/v(x)
154 -0.100006103515625 <= x <= 0.231639862060546875
155 peak relative error 1.3e-21 */
156 u0
= -8.886217500092090678492242071879342025627E1L
,
157 u1
= 6.840109978129177639438792958320783599310E2L
,
158 u2
= 2.042626104514127267855588786511809932433E3L
,
159 u3
= 1.911723903442667422201651063009856064275E3L
,
160 u4
= 7.447065275665887457628865263491667767695E2L
,
161 u5
= 1.132256494121790736268471016493103952637E2L
,
162 u6
= 4.484398885516614191003094714505960972894E0L
,
164 v0
= 1.150830924194461522996462401210374632929E3L
,
165 v1
= 3.399692260848747447377972081399737098610E3L
,
166 v2
= 3.786631705644460255229513563657226008015E3L
,
167 v3
= 1.966450123004478374557778781564114347876E3L
,
168 v4
= 4.741359068914069299837355438370682773122E2L
,
169 v5
= 4.508989649747184050907206782117647852364E1L
,
170 /* v6 = 1.000000000000000000000000000000000000000E0 */
173 /* lgam (x+2) = .5 x + x s(x)/r(x)
175 peak relative error 7.2e-22 */
176 s0
= 1.454726263410661942989109455292824853344E6L
,
177 s1
= -3.901428390086348447890408306153378922752E6L
,
178 s2
= -6.573568698209374121847873064292963089438E6L
,
179 s3
= -3.319055881485044417245964508099095984643E6L
,
180 s4
= -7.094891568758439227560184618114707107977E5L
,
181 s5
= -6.263426646464505837422314539808112478303E4L
,
182 s6
= -1.684926520999477529949915657519454051529E3L
,
184 r0
= -1.883978160734303518163008696712983134698E7L
,
185 r1
= -2.815206082812062064902202753264922306830E7L
,
186 r2
= -1.600245495251915899081846093343626358398E7L
,
187 r3
= -4.310526301881305003489257052083370058799E6L
,
188 r4
= -5.563807682263923279438235987186184968542E5L
,
189 r5
= -3.027734654434169996032905158145259713083E4L
,
190 r6
= -4.501995652861105629217250715790764371267E2L
,
191 /* r6 = 1.000000000000000000000000000000000000000E0 */
194 /* lgam(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x w(1/x^2)
196 Peak relative error 1.51e-21
198 w0
= 4.189385332046727417803e-1L,
199 w1
= 8.333333333333331447505E-2L,
200 w2
= -2.777777777750349603440E-3L,
201 w3
= 7.936507795855070755671E-4L,
202 w4
= -5.952345851765688514613E-4L,
203 w5
= 8.412723297322498080632E-4L,
204 w6
= -1.880801938119376907179E-3L,
205 w7
= 4.885026142432270781165E-3L;
208 static const long double zero
= 0.0L;
210 static long double zero
= 0.0L;
215 sin_pi (long double x
)
224 u_int32_t se
, i0
, i1
;
226 GET_LDOUBLE_WORDS (se
, i0
, i1
, x
);
228 ix
= (ix
<< 16) | (i0
>> 16);
229 if (ix
< 0x3ffd8000) /* 0.25 */
230 return __sinl (pi
* x
);
231 y
= -x
; /* x is assume negative */
234 * argument reduction, make sure inexact flag not raised if input
239 { /* inexact anyway */
241 y
= 2.0*(y
- __floorl(y
)); /* y = |x| mod 2.0 */
246 if (ix
>= 0x403f8000) /* 2^64 */
248 y
= zero
; n
= 0; /* y must be even */
252 if (ix
< 0x403e8000) /* 2^63 */
253 z
= y
+ two63
; /* exact */
254 GET_LDOUBLE_WORDS (se
, i0
, i1
, z
);
268 y
= __cosl (pi
* (half
- y
));
272 y
= __sinl (pi
* (one
- y
));
276 y
= -__cosl (pi
* (y
- 1.5));
279 y
= __sinl (pi
* (y
- 2.0));
288 __ieee754_lgammal_r (long double x
, int *signgamp
)
291 __ieee754_lgammal_r (x
, signgamp
)
296 long double t
, y
, z
, nadj
, p
, p1
, p2
, q
, r
, w
;
298 u_int32_t se
, i0
, i1
;
301 GET_LDOUBLE_WORDS (se
, i0
, i1
, x
);
304 if ((ix
| i0
| i1
) == 0)
308 return one
/ fabsl (x
);
311 ix
= (ix
<< 16) | (i0
>> 16);
313 /* purge off +-inf, NaN, +-0, and negative arguments */
314 if (ix
>= 0x7fff0000)
317 if (ix
< 0x3fc08000) /* 2^-63 */
318 { /* |x|<2**-63, return -log(|x|) */
322 return -__ieee754_logl (-x
);
325 return -__ieee754_logl (x
);
331 return one
/ fabsl (t
); /* -integer */
332 nadj
= __ieee754_logl (pi
/ fabsl (t
* x
));
338 /* purge off 1 and 2 */
339 if ((((ix
- 0x3fff8000) | i0
| i1
) == 0)
340 || (((ix
- 0x40008000) | i0
| i1
) == 0))
342 else if (ix
< 0x40008000) /* 2.0 */
345 if (ix
<= 0x3ffee666) /* 8.99993896484375e-1 */
347 /* lgamma(x) = lgamma(x+1) - log(x) */
348 r
= -__ieee754_logl (x
);
349 if (ix
>= 0x3ffebb4a) /* 7.31597900390625e-1 */
354 else if (ix
>= 0x3ffced33)/* 2.31639862060546875e-1 */
369 if (ix
>= 0x3fffdda6) /* 1.73162841796875 */
375 else if (ix
>= 0x3fff9da6)/* 1.23162841796875 */
391 p1
= a0
+ y
* (a1
+ y
* (a2
+ y
* (a3
+ y
* (a4
+ y
* a5
))));
392 p2
= b0
+ y
* (b1
+ y
* (b2
+ y
* (b3
+ y
* (b4
+ y
))));
393 r
+= half
* y
+ y
* p1
/p2
;
396 p1
= g0
+ y
* (g1
+ y
* (g2
+ y
* (g3
+ y
* (g4
+ y
* (g5
+ y
* g6
)))));
397 p2
= h0
+ y
* (h1
+ y
* (h2
+ y
* (h3
+ y
* (h4
+ y
* (h5
+ y
)))));
402 p1
= y
* (u0
+ y
* (u1
+ y
* (u2
+ y
* (u3
+ y
* (u4
+ y
* (u5
+ y
* u6
))))));
403 p2
= v0
+ y
* (v1
+ y
* (v2
+ y
* (v3
+ y
* (v4
+ y
* (v5
+ y
)))));
404 r
+= (-half
* y
+ p1
/ p2
);
407 else if (ix
< 0x40028000) /* 8.0 */
414 (s0
+ y
* (s1
+ y
* (s2
+ y
* (s3
+ y
* (s4
+ y
* (s5
+ y
* s6
))))));
415 q
= r0
+ y
* (r1
+ y
* (r2
+ y
* (r3
+ y
* (r4
+ y
* (r5
+ y
* (r6
+ y
))))));
416 r
= half
* y
+ p
/ q
;
417 z
= one
; /* lgamma(1+s) = log(s) + lgamma(s) */
421 z
*= (y
+ 6.0); /* FALLTHRU */
423 z
*= (y
+ 5.0); /* FALLTHRU */
425 z
*= (y
+ 4.0); /* FALLTHRU */
427 z
*= (y
+ 3.0); /* FALLTHRU */
429 z
*= (y
+ 2.0); /* FALLTHRU */
430 r
+= __ieee754_logl (z
);
434 else if (ix
< 0x40418000) /* 2^66 */
436 /* 8.0 <= x < 2**66 */
437 t
= __ieee754_logl (x
);
441 + y
* (w2
+ y
* (w3
+ y
* (w4
+ y
* (w5
+ y
* (w6
+ y
* w7
))))));
442 r
= (x
- half
) * (t
- one
) + w
;
445 /* 2**66 <= x <= inf */
446 r
= x
* (__ieee754_logl (x
) - one
);