1 /* e_lgammaf_r.c -- float version of e_lgamma_r.c.
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
9 * Developed at SunPro, a Sun Microsystems, Inc. business.
10 * Permission to use, copy, modify, and distribute this
11 * software is freely granted, provided that this notice
13 * ====================================================
16 #include <libc-internal.h>
18 #include <math_private.h>
21 two23
= 8.3886080000e+06, /* 0x4b000000 */
22 half
= 5.0000000000e-01, /* 0x3f000000 */
23 one
= 1.0000000000e+00, /* 0x3f800000 */
24 pi
= 3.1415927410e+00, /* 0x40490fdb */
25 a0
= 7.7215664089e-02, /* 0x3d9e233f */
26 a1
= 3.2246702909e-01, /* 0x3ea51a66 */
27 a2
= 6.7352302372e-02, /* 0x3d89f001 */
28 a3
= 2.0580807701e-02, /* 0x3ca89915 */
29 a4
= 7.3855509982e-03, /* 0x3bf2027e */
30 a5
= 2.8905137442e-03, /* 0x3b3d6ec6 */
31 a6
= 1.1927076848e-03, /* 0x3a9c54a1 */
32 a7
= 5.1006977446e-04, /* 0x3a05b634 */
33 a8
= 2.2086278477e-04, /* 0x39679767 */
34 a9
= 1.0801156895e-04, /* 0x38e28445 */
35 a10
= 2.5214456400e-05, /* 0x37d383a2 */
36 a11
= 4.4864096708e-05, /* 0x383c2c75 */
37 tc
= 1.4616321325e+00, /* 0x3fbb16c3 */
38 tf
= -1.2148628384e-01, /* 0xbdf8cdcd */
39 /* tt = -(tail of tf) */
40 tt
= 6.6971006518e-09, /* 0x31e61c52 */
41 t0
= 4.8383611441e-01, /* 0x3ef7b95e */
42 t1
= -1.4758771658e-01, /* 0xbe17213c */
43 t2
= 6.4624942839e-02, /* 0x3d845a15 */
44 t3
= -3.2788541168e-02, /* 0xbd064d47 */
45 t4
= 1.7970675603e-02, /* 0x3c93373d */
46 t5
= -1.0314224288e-02, /* 0xbc28fcfe */
47 t6
= 6.1005386524e-03, /* 0x3bc7e707 */
48 t7
= -3.6845202558e-03, /* 0xbb7177fe */
49 t8
= 2.2596477065e-03, /* 0x3b141699 */
50 t9
= -1.4034647029e-03, /* 0xbab7f476 */
51 t10
= 8.8108185446e-04, /* 0x3a66f867 */
52 t11
= -5.3859531181e-04, /* 0xba0d3085 */
53 t12
= 3.1563205994e-04, /* 0x39a57b6b */
54 t13
= -3.1275415677e-04, /* 0xb9a3f927 */
55 t14
= 3.3552918467e-04, /* 0x39afe9f7 */
56 u0
= -7.7215664089e-02, /* 0xbd9e233f */
57 u1
= 6.3282704353e-01, /* 0x3f2200f4 */
58 u2
= 1.4549225569e+00, /* 0x3fba3ae7 */
59 u3
= 9.7771751881e-01, /* 0x3f7a4bb2 */
60 u4
= 2.2896373272e-01, /* 0x3e6a7578 */
61 u5
= 1.3381091878e-02, /* 0x3c5b3c5e */
62 v1
= 2.4559779167e+00, /* 0x401d2ebe */
63 v2
= 2.1284897327e+00, /* 0x4008392d */
64 v3
= 7.6928514242e-01, /* 0x3f44efdf */
65 v4
= 1.0422264785e-01, /* 0x3dd572af */
66 v5
= 3.2170924824e-03, /* 0x3b52d5db */
67 s0
= -7.7215664089e-02, /* 0xbd9e233f */
68 s1
= 2.1498242021e-01, /* 0x3e5c245a */
69 s2
= 3.2577878237e-01, /* 0x3ea6cc7a */
70 s3
= 1.4635047317e-01, /* 0x3e15dce6 */
71 s4
= 2.6642270386e-02, /* 0x3cda40e4 */
72 s5
= 1.8402845599e-03, /* 0x3af135b4 */
73 s6
= 3.1947532989e-05, /* 0x3805ff67 */
74 r1
= 1.3920053244e+00, /* 0x3fb22d3b */
75 r2
= 7.2193557024e-01, /* 0x3f38d0c5 */
76 r3
= 1.7193385959e-01, /* 0x3e300f6e */
77 r4
= 1.8645919859e-02, /* 0x3c98bf54 */
78 r5
= 7.7794247773e-04, /* 0x3a4beed6 */
79 r6
= 7.3266842264e-06, /* 0x36f5d7bd */
80 w0
= 4.1893854737e-01, /* 0x3ed67f1d */
81 w1
= 8.3333335817e-02, /* 0x3daaaaab */
82 w2
= -2.7777778450e-03, /* 0xbb360b61 */
83 w3
= 7.9365057172e-04, /* 0x3a500cfd */
84 w4
= -5.9518753551e-04, /* 0xba1c065c */
85 w5
= 8.3633989561e-04, /* 0x3a5b3dd2 */
86 w6
= -1.6309292987e-03; /* 0xbad5c4e8 */
88 static const float zero
= 0.0000000000e+00;
99 if(ix
<0x3e800000) return __kernel_sinf(pi
*x
,zero
,0);
100 y
= -x
; /* x is assume negative */
103 * argument reduction, make sure inexact flag not raised if input
107 if(z
!=y
) { /* inexact anyway */
109 y
= (float)2.0*(y
- __floorf(y
)); /* y = |x| mod 2.0 */
110 n
= (int) (y
*(float)4.0);
113 y
= zero
; n
= 0; /* y must be even */
115 if(ix
<0x4b000000) z
= y
+two23
; /* exact */
123 case 0: y
= __kernel_sinf(pi
*y
,zero
,0); break;
125 case 2: y
= __kernel_cosf(pi
*((float)0.5-y
),zero
); break;
127 case 4: y
= __kernel_sinf(pi
*(one
-y
),zero
,0); break;
129 case 6: y
= -__kernel_cosf(pi
*(y
-(float)1.5),zero
); break;
130 default: y
= __kernel_sinf(pi
*(y
-(float)2.0),zero
,0); break;
137 __ieee754_lgammaf_r(float x
, int *signgamp
)
139 float t
,y
,z
,nadj
,p
,p1
,p2
,p3
,q
,r
,w
;
142 GET_FLOAT_WORD(hx
,x
);
144 /* purge off +-inf, NaN, +-0, and negative arguments */
147 if(__builtin_expect(ix
>=0x7f800000, 0)) return x
*x
;
148 if(__builtin_expect(ix
==0, 0))
154 if(__builtin_expect(ix
<0x30800000, 0)) {
155 /* |x|<2**-30, return -log(|x|) */
158 return -__ieee754_logf(-x
);
159 } else return -__ieee754_logf(x
);
162 if(ix
>=0x4b000000) /* |x|>=2**23, must be -integer */
165 if(t
==zero
) return one
/fabsf(t
); /* -integer */
166 nadj
= __ieee754_logf(pi
/fabsf(t
*x
));
167 if(t
<zero
) *signgamp
= -1;
171 /* purge off 1 and 2 */
172 if (ix
==0x3f800000||ix
==0x40000000) r
= 0;
174 else if(ix
<0x40000000) {
175 if(ix
<=0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */
176 r
= -__ieee754_logf(x
);
177 if(ix
>=0x3f3b4a20) {y
= one
-x
; i
= 0;}
178 else if(ix
>=0x3e6d3308) {y
= x
-(tc
-one
); i
=1;}
182 if(ix
>=0x3fdda618) {y
=(float)2.0-x
;i
=0;} /* [1.7316,2] */
183 else if(ix
>=0x3F9da620) {y
=x
-tc
;i
=1;} /* [1.23,1.73] */
189 p1
= a0
+z
*(a2
+z
*(a4
+z
*(a6
+z
*(a8
+z
*a10
))));
190 p2
= z
*(a1
+z
*(a3
+z
*(a5
+z
*(a7
+z
*(a9
+z
*a11
)))));
192 r
+= (p
-(float)0.5*y
); break;
196 p1
= t0
+w
*(t3
+w
*(t6
+w
*(t9
+w
*t12
))); /* parallel comp */
197 p2
= t1
+w
*(t4
+w
*(t7
+w
*(t10
+w
*t13
)));
198 p3
= t2
+w
*(t5
+w
*(t8
+w
*(t11
+w
*t14
)));
199 p
= z
*p1
-(tt
-w
*(p2
+y
*p3
));
200 r
+= (tf
+ p
); break;
202 p1
= y
*(u0
+y
*(u1
+y
*(u2
+y
*(u3
+y
*(u4
+y
*u5
)))));
203 p2
= one
+y
*(v1
+y
*(v2
+y
*(v3
+y
*(v4
+y
*v5
))));
204 r
+= (-(float)0.5*y
+ p1
/p2
);
207 else if(ix
<0x41000000) { /* x < 8.0 */
211 p
= y
*(s0
+y
*(s1
+y
*(s2
+y
*(s3
+y
*(s4
+y
*(s5
+y
*s6
))))));
212 q
= one
+y
*(r1
+y
*(r2
+y
*(r3
+y
*(r4
+y
*(r5
+y
*r6
)))));
214 z
= one
; /* lgamma(1+s) = log(s) + lgamma(s) */
216 case 7: z
*= (y
+(float)6.0); /* FALLTHRU */
217 case 6: z
*= (y
+(float)5.0); /* FALLTHRU */
218 case 5: z
*= (y
+(float)4.0); /* FALLTHRU */
219 case 4: z
*= (y
+(float)3.0); /* FALLTHRU */
220 case 3: z
*= (y
+(float)2.0); /* FALLTHRU */
221 r
+= __ieee754_logf(z
); break;
223 /* 8.0 <= x < 2**26 */
224 } else if (ix
< 0x4c800000) {
225 t
= __ieee754_logf(x
);
228 w
= w0
+z
*(w1
+y
*(w2
+y
*(w3
+y
*(w4
+y
*(w5
+y
*w6
)))));
229 r
= (x
-half
)*(t
-one
)+w
;
231 /* 2**26 <= x <= inf */
232 r
= x
*(__ieee754_logf(x
)-one
);
233 /* NADJ is set for negative arguments but not otherwise,
234 resulting in warnings that it may be used uninitialized
235 although in the cases where it is used it has always been
237 DIAG_PUSH_NEEDS_COMMENT
;
238 #if __GNUC_PREREQ (4, 7)
239 DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wmaybe-uninitialized");
241 DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wuninitialized");
243 if(hx
<0) r
= nadj
- r
;
244 DIAG_POP_NEEDS_COMMENT
;
247 strong_alias (__ieee754_lgammaf_r
, __lgammaf_r_finite
)