1 /* @(#)e_lgamma_r.c 1.3 95/01/18 */
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 * Developed at SunSoft, 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 //__FBSDID("$FreeBSD$");
16 * See e_lgamma_r.c for complete comments.
18 * Converted to long double by Steven G. Kargl.
23 #include "math_private.h"
25 #if defined(__x86_64__)
26 #include "x86_64/ieeefp.h"
27 #elif defined(__i386__)
28 #include "i386/ieeefp.h"
31 static const volatile double vzero
__attribute__ ((__section__(".rodata"))) = 0;
38 static const union IEEEl2bits
39 piu
= LD80C(0xc90fdaa22168c235, 1, 3.14159265358979323851e+00L);
42 * Domain y in [0x1p-70, 0.27], range ~[-4.5264e-22, 4.5264e-22]:
43 * |(lgamma(2 - y) + y / 2) / y - a(y)| < 2**-70.9
45 static const union IEEEl2bits
46 a0u
= LD80C(0x9e233f1bed863d26, -4, 7.72156649015328606028e-02L),
47 a1u
= LD80C(0xa51a6625307d3249, -2, 3.22467033424113218889e-01L),
48 a2u
= LD80C(0x89f000d2abafda8c, -4, 6.73523010531979398946e-02L),
49 a3u
= LD80C(0xa8991563eca75f26, -6, 2.05808084277991211934e-02L),
50 a4u
= LD80C(0xf2027e10634ce6b6, -8, 7.38555102796070454026e-03L),
51 a5u
= LD80C(0xbd6eb76dd22187f4, -9, 2.89051035162703932972e-03L),
52 a6u
= LD80C(0x9c562ab05e0458ed, -10, 1.19275351624639999297e-03L),
53 a7u
= LD80C(0x859baed93ee48e46, -11, 5.09674593842117925320e-04L),
54 a8u
= LD80C(0xe9f28a4432949af2, -13, 2.23109648015769155122e-04L),
55 a9u
= LD80C(0xd12ad0d9b93c6bb0, -14, 9.97387167479808509830e-05L),
56 a10u
= LD80C(0xb7522643c78a219b, -15, 4.37071076331030136818e-05L),
57 a11u
= LD80C(0xca024dcdece2cb79, -16, 2.40813493372040143061e-05L),
58 a12u
= LD80C(0xbb90fb6968ebdbf9, -19, 2.79495621083634031729e-06L),
59 a13u
= LD80C(0xba1c9ffeeae07b37, -17, 1.10931287015513924136e-05L);
75 * Domain x in [tc-0.24, tc+0.28], range ~[-6.1205e-22, 6.1205e-22]:
76 * |(lgamma(x) - tf) - t(x - tc)| < 2**-70.5
78 static const union IEEEl2bits
79 tcu
= LD80C(0xbb16c31ab5f1fb71, 0, 1.46163214496836234128e+00L),
80 tfu
= LD80C(0xf8cdcde61c520e0f, -4, -1.21486290535849608093e-01L),
81 ttu
= LD80C(0xd46ee54b27d4de99, -69, -2.81152980996018785880e-21L),
82 t0u
= LD80C(0x80b9406556a62a6b, -68, 3.40728634996055147231e-21L),
83 t1u
= LD80C(0xc7e9c6f6df3f8c39, -67, -1.05833162742737073665e-20L),
84 t2u
= LD80C(0xf7b95e4771c55d51, -2, 4.83836122723810583532e-01L),
85 t3u
= LD80C(0x97213c6e35e119ff, -3, -1.47587722994530691476e-01L),
86 t4u
= LD80C(0x845a14a6a81dc94b, -4, 6.46249402389135358063e-02L),
87 t5u
= LD80C(0x864d46fa89997796, -5, -3.27885410884846056084e-02L),
88 t6u
= LD80C(0x93373cbd00297438, -6, 1.79706751150707171293e-02L),
89 t7u
= LD80C(0xa8fcfca7eddc8d1d, -7, -1.03142230361450732547e-02L),
90 t8u
= LD80C(0xc7e7015ff4bc45af, -8, 6.10053603296546099193e-03L),
91 t9u
= LD80C(0xf178d2247adc5093, -9, -3.68456964904901200152e-03L),
92 t10u
= LD80C(0x94188d58f12e5e9f, -9, 2.25976420273774583089e-03L),
93 t11u
= LD80C(0xb7cbaef14e1406f1, -10, -1.40224943666225639823e-03L),
94 t12u
= LD80C(0xe63a671e6704ea4d, -11, 8.78250640744776944887e-04L),
95 t13u
= LD80C(0x914b6c9cae61783e, -11, -5.54255012657716808811e-04L),
96 t14u
= LD80C(0xb858f5bdb79276fe, -12, 3.51614951536825927370e-04L),
97 t15u
= LD80C(0xea73e744c34b9591, -13, -2.23591563824520112236e-04L),
98 t16u
= LD80C(0x99aeabb0d67ba835, -13, 1.46562869351659194136e-04L),
99 t17u
= LD80C(0xd7c6938325db2024, -14, -1.02889866046435680588e-04L),
100 t18u
= LD80C(0xe24cb1e3b0474775, -15, 5.39540265505221957652e-05L);
124 * Domain y in [-0.1, 0.232], range ~[-8.1938e-22, 8.3815e-22]:
125 * |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-71.2
127 static const union IEEEl2bits
128 u0u
= LD80C(0x9e233f1bed863d27, -4, -7.72156649015328606095e-02L),
129 u1u
= LD80C(0x98280ee45e4ddd3d, -1, 5.94361239198682739769e-01L),
130 u2u
= LD80C(0xe330c8ead4130733, 0, 1.77492629495841234275e+00L),
131 u3u
= LD80C(0xd4a213f1a002ec52, 0, 1.66119622514818078064e+00L),
132 u4u
= LD80C(0xa5a9ca6f5bc62163, -1, 6.47122051417476492989e-01L),
133 u5u
= LD80C(0xc980e49cd5b019e6, -4, 9.83903751718671509455e-02L),
134 u6u
= LD80C(0xff636a8bdce7025b, -9, 3.89691687802305743450e-03L),
135 v1u
= LD80C(0xbd109c533a19fbf5, 1, 2.95413883330948556544e+00L),
136 v2u
= LD80C(0xd295cbf96f31f099, 1, 3.29039286955665403176e+00L),
137 v3u
= LD80C(0xdab8bcfee40496cb, 0, 1.70876276441416471410e+00L),
138 v4u
= LD80C(0xd2f2dc3638567e9f, -2, 4.12009126299534668571e-01L),
139 v5u
= LD80C(0xa07d9b0851070f41, -5, 3.91822868305682491442e-02L),
140 v6u
= LD80C(0xe3cd8318f7adb2c4, -11, 8.68998648222144351114e-04L);
155 * Domain x in (2, 3], range ~[-3.3648e-22, 3.4416e-22]:
156 * |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-72.3
159 static const union IEEEl2bits
160 s0u
= LD80C(0x9e233f1bed863d27, -4, -7.72156649015328606095e-02L),
161 s1u
= LD80C(0xd3ff0dcc7fa91f94, -3, 2.07027640921219389860e-01L),
162 s2u
= LD80C(0xb2bb62782478ef31, -2, 3.49085881391362090549e-01L),
163 s3u
= LD80C(0xb49f7438c4611a74, -3, 1.76389518704213357954e-01L),
164 s4u
= LD80C(0x9a957008fa27ecf9, -5, 3.77401710862930008071e-02L),
165 s5u
= LD80C(0xda9b389a6ca7a7ac, -9, 3.33566791452943399399e-03L),
166 s6u
= LD80C(0xbc7a2263faf59c14, -14, 8.98728786745638844395e-05L),
167 r1u
= LD80C(0xbf5cff5b11477d4d, 0, 1.49502555796294337722e+00L),
168 r2u
= LD80C(0xd9aec89de08e3da6, -1, 8.50323236984473285866e-01L),
169 r3u
= LD80C(0xeab7ae5057c443f9, -3, 2.29216312078225806131e-01L),
170 r4u
= LD80C(0xf29707d9bd2b1e37, -6, 2.96130326586640089145e-02L),
171 r5u
= LD80C(0xd376c2f09736c5a3, -10, 1.61334161411590662495e-03L),
172 r6u
= LD80C(0xc985983d0cd34e3d, -16, 2.40232770710953450636e-05L),
173 r7u
= LD80C(0xe5c7a4f7fc2ef13d, -25, -5.34997929289167573510e-08L);
189 * Domain z in [8, 0x1p70], range ~[-3.0235e-22, 3.0563e-22]:
190 * |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-71.7
192 static const union IEEEl2bits
193 w0u
= LD80C(0xd67f1c864beb4a69, -2, 4.18938533204672741776e-01L),
194 w1u
= LD80C(0xaaaaaaaaaaaaaaa1, -4, 8.33333333333333332678e-02L),
195 w2u
= LD80C(0xb60b60b60b5491c9, -9, -2.77777777777760927870e-03L),
196 w3u
= LD80C(0xd00d00cf58aede4c, -11, 7.93650793490637233668e-04L),
197 w4u
= LD80C(0x9c09bf626783d4a5, -11, -5.95238023926039051268e-04L),
198 w5u
= LD80C(0xdca7cadc5baa517b, -11, 8.41733700408000822962e-04L),
199 w6u
= LD80C(0xfb060e361e1ffd07, -10, -1.91515849570245136604e-03L),
200 w7u
= LD80C(0xcbd5101bb58d1f2b, -8, 6.22046743903262649294e-03L),
201 w8u
= LD80C(0xad27a668d32c821b, -6, -2.11370706734662081843e-02L);
213 sin_pil(long double x
)
215 volatile long double vz
;
218 __unused
uint16_t hx
;
228 EXTRACT_LDBL80_WORDS(hx
,n
,vz
);
231 z
-= 0.25; /* adjust to round down */
234 n
&= 7; /* octant of y mod 2 */
235 y
= y
- z
+ n
* 0.25; /* y mod 2 */
238 case 0: y
= __kernel_sinl(pi
*y
,zero
,0); break;
240 case 2: y
= __kernel_cosl(pi
*(0.5-y
),zero
); break;
242 case 4: y
= __kernel_sinl(pi
*(one
-y
),zero
,0); break;
244 case 6: y
= -__kernel_cosl(pi
*(y
-1.5),zero
); break;
245 default: y
= __kernel_sinl(pi
*(y
-2.0),zero
,0); break;
251 lgammal_r(long double x
, int *signgamp
)
253 long double nadj
,p
,p1
,p2
,q
,r
,t
,w
,y
,z
;
258 EXTRACT_LDBL80_WORDS(hx
,lx
,x
);
260 /* purge +-Inf and NaNs */
263 if(ix
==0x7fff) return x
*x
;
267 /* purge +-0 and tiny arguments */
268 *signgamp
= 1-2*(hx
>>15);
269 if(ix
<0x3fff-67) { /* |x|<2**-(p+3), return -log(|x|) */
272 RETURNI(-logl(fabsl(x
)));
275 /* purge negative integers and start evaluation for other x < 0 */
278 if(ix
>=0x3fff+63) /* |x|>=2**(p-1), must be -integer */
281 if(t
==zero
) RETURNI(one
/vzero
); /* -integer */
282 nadj
= logl(pi
/fabsl(t
*x
));
283 if(t
<zero
) *signgamp
= -1;
288 if((ix
==0x3fff || ix
==0x4000) && lx
==0x8000000000000000ULL
) r
= 0;
292 * XXX Supposedly, one can use the following information to replace the
293 * XXX FP rational expressions. A similar approach is appropriate
294 * XXX for ld128, but one (may need?) needs to consider llx, too.
296 * 8.9999961853027344e-01 3ffe e666600000000000
297 * 7.3159980773925781e-01 3ffe bb4a200000000000
298 * 2.3163998126983643e-01 3ffc ed33080000000000
299 * 1.7316312789916992e+00 3fff dda6180000000000
300 * 1.2316322326660156e+00 3fff 9da6200000000000
302 if(x
<8.9999961853027344e-01) {
304 if(x
>=7.3159980773925781e-01) {y
= 1-x
; i
= 0;}
305 else if(x
>=2.3163998126983643e-01) {y
= x
-(tc
-1); i
=1;}
309 if(x
>=1.7316312789916992e+00) {y
=2-x
;i
=0;}
310 else if(x
>=1.2316322326660156e+00) {y
=x
-tc
;i
=1;}
316 p1
= a0
+z
*(a2
+z
*(a4
+z
*(a6
+z
*(a8
+z
*(a10
+z
*a12
)))));
317 p2
= z
*(a1
+z
*(a3
+z
*(a5
+z
*(a7
+z
*(a9
+z
*(a11
+z
*a13
))))));
321 p
= t0
+y
*t1
+tt
+y
*y
*(t2
+y
*(t3
+y
*(t4
+y
*(t5
+y
*(t6
+y
*(t7
+y
*(t8
+
322 y
*(t9
+y
*(t10
+y
*(t11
+y
*(t12
+y
*(t13
+y
*(t14
+y
*(t15
+y
*(t16
+
323 y
*(t17
+y
*t18
))))))))))))))));
326 p1
= y
*(u0
+y
*(u1
+y
*(u2
+y
*(u3
+y
*(u4
+y
*(u5
+y
*u6
))))));
327 p2
= 1+y
*(v1
+y
*(v2
+y
*(v3
+y
*(v4
+y
*(v5
+y
*v6
)))));
335 p
= y
*(s0
+y
*(s1
+y
*(s2
+y
*(s3
+y
*(s4
+y
*(s5
+y
*s6
))))));
336 q
= 1+y
*(r1
+y
*(r2
+y
*(r3
+y
*(r4
+y
*(r5
+y
*(r6
+y
*r7
))))));
338 z
= 1; /* lgamma(1+s) = log(s) + lgamma(s) */
340 case 7: z
*= (y
+6); /* FALLTHRU */
341 case 6: z
*= (y
+5); /* FALLTHRU */
342 case 5: z
*= (y
+4); /* FALLTHRU */
343 case 4: z
*= (y
+3); /* FALLTHRU */
344 case 3: z
*= (y
+2); /* FALLTHRU */
347 /* 8.0 <= x < 2**(p+3) */
348 } else if (ix
<0x3fff+67) {
352 w
= w0
+z
*(w1
+y
*(w2
+y
*(w3
+y
*(w4
+y
*(w5
+y
*(w6
+y
*(w7
+y
*w8
)))))));
353 r
= (x
-half
)*(t
-one
)+w
;
354 /* 2**(p+3) <= x <= inf */
357 if(hx
&0x8000) r
= nadj
- r
;