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 * ====================================================
17 #include <math-narrow-eval.h>
18 #include <math_private.h>
19 #include <libc-diag.h>
22 two23
= 8.3886080000e+06, /* 0x4b000000 */
23 half
= 5.0000000000e-01, /* 0x3f000000 */
24 one
= 1.0000000000e+00, /* 0x3f800000 */
25 pi
= 3.1415927410e+00, /* 0x40490fdb */
26 a0
= 7.7215664089e-02, /* 0x3d9e233f */
27 a1
= 3.2246702909e-01, /* 0x3ea51a66 */
28 a2
= 6.7352302372e-02, /* 0x3d89f001 */
29 a3
= 2.0580807701e-02, /* 0x3ca89915 */
30 a4
= 7.3855509982e-03, /* 0x3bf2027e */
31 a5
= 2.8905137442e-03, /* 0x3b3d6ec6 */
32 a6
= 1.1927076848e-03, /* 0x3a9c54a1 */
33 a7
= 5.1006977446e-04, /* 0x3a05b634 */
34 a8
= 2.2086278477e-04, /* 0x39679767 */
35 a9
= 1.0801156895e-04, /* 0x38e28445 */
36 a10
= 2.5214456400e-05, /* 0x37d383a2 */
37 a11
= 4.4864096708e-05, /* 0x383c2c75 */
38 tc
= 1.4616321325e+00, /* 0x3fbb16c3 */
39 tf
= -1.2148628384e-01, /* 0xbdf8cdcd */
40 /* tt = -(tail of tf) */
41 tt
= 6.6971006518e-09, /* 0x31e61c52 */
42 t0
= 4.8383611441e-01, /* 0x3ef7b95e */
43 t1
= -1.4758771658e-01, /* 0xbe17213c */
44 t2
= 6.4624942839e-02, /* 0x3d845a15 */
45 t3
= -3.2788541168e-02, /* 0xbd064d47 */
46 t4
= 1.7970675603e-02, /* 0x3c93373d */
47 t5
= -1.0314224288e-02, /* 0xbc28fcfe */
48 t6
= 6.1005386524e-03, /* 0x3bc7e707 */
49 t7
= -3.6845202558e-03, /* 0xbb7177fe */
50 t8
= 2.2596477065e-03, /* 0x3b141699 */
51 t9
= -1.4034647029e-03, /* 0xbab7f476 */
52 t10
= 8.8108185446e-04, /* 0x3a66f867 */
53 t11
= -5.3859531181e-04, /* 0xba0d3085 */
54 t12
= 3.1563205994e-04, /* 0x39a57b6b */
55 t13
= -3.1275415677e-04, /* 0xb9a3f927 */
56 t14
= 3.3552918467e-04, /* 0x39afe9f7 */
57 u0
= -7.7215664089e-02, /* 0xbd9e233f */
58 u1
= 6.3282704353e-01, /* 0x3f2200f4 */
59 u2
= 1.4549225569e+00, /* 0x3fba3ae7 */
60 u3
= 9.7771751881e-01, /* 0x3f7a4bb2 */
61 u4
= 2.2896373272e-01, /* 0x3e6a7578 */
62 u5
= 1.3381091878e-02, /* 0x3c5b3c5e */
63 v1
= 2.4559779167e+00, /* 0x401d2ebe */
64 v2
= 2.1284897327e+00, /* 0x4008392d */
65 v3
= 7.6928514242e-01, /* 0x3f44efdf */
66 v4
= 1.0422264785e-01, /* 0x3dd572af */
67 v5
= 3.2170924824e-03, /* 0x3b52d5db */
68 s0
= -7.7215664089e-02, /* 0xbd9e233f */
69 s1
= 2.1498242021e-01, /* 0x3e5c245a */
70 s2
= 3.2577878237e-01, /* 0x3ea6cc7a */
71 s3
= 1.4635047317e-01, /* 0x3e15dce6 */
72 s4
= 2.6642270386e-02, /* 0x3cda40e4 */
73 s5
= 1.8402845599e-03, /* 0x3af135b4 */
74 s6
= 3.1947532989e-05, /* 0x3805ff67 */
75 r1
= 1.3920053244e+00, /* 0x3fb22d3b */
76 r2
= 7.2193557024e-01, /* 0x3f38d0c5 */
77 r3
= 1.7193385959e-01, /* 0x3e300f6e */
78 r4
= 1.8645919859e-02, /* 0x3c98bf54 */
79 r5
= 7.7794247773e-04, /* 0x3a4beed6 */
80 r6
= 7.3266842264e-06, /* 0x36f5d7bd */
81 w0
= 4.1893854737e-01, /* 0x3ed67f1d */
82 w1
= 8.3333335817e-02, /* 0x3daaaaab */
83 w2
= -2.7777778450e-03, /* 0xbb360b61 */
84 w3
= 7.9365057172e-04, /* 0x3a500cfd */
85 w4
= -5.9518753551e-04, /* 0xba1c065c */
86 w5
= 8.3633989561e-04, /* 0x3a5b3dd2 */
87 w6
= -1.6309292987e-03; /* 0xbad5c4e8 */
89 static const float zero
= 0.0000000000e+00;
100 if(ix
<0x3e800000) return __sinf (pi
*x
);
101 y
= -x
; /* x is assume negative */
104 * argument reduction, make sure inexact flag not raised if input
108 if(z
!=y
) { /* inexact anyway */
110 y
= (float)2.0*(y
- __floorf(y
)); /* y = |x| mod 2.0 */
111 n
= (int) (y
*(float)4.0);
114 y
= zero
; n
= 0; /* y must be even */
116 if(ix
<0x4b000000) z
= y
+two23
; /* exact */
124 case 0: y
= __sinf (pi
*y
); break;
126 case 2: y
= __cosf (pi
*((float)0.5-y
)); break;
128 case 4: y
= __sinf (pi
*(one
-y
)); break;
130 case 6: y
= -__cosf (pi
*(y
-(float)1.5)); break;
131 default: y
= __sinf (pi
*(y
-(float)2.0)); break;
138 __ieee754_lgammaf_r(float x
, int *signgamp
)
140 float t
,y
,z
,nadj
,p
,p1
,p2
,p3
,q
,r
,w
;
143 GET_FLOAT_WORD(hx
,x
);
145 /* purge off +-inf, NaN, +-0, and negative arguments */
148 if(__builtin_expect(ix
>=0x7f800000, 0)) return x
*x
;
149 if(__builtin_expect(ix
==0, 0))
155 if(__builtin_expect(ix
<0x30800000, 0)) {
156 /* |x|<2**-30, return -log(|x|) */
159 return -__ieee754_logf(-x
);
160 } else return -__ieee754_logf(x
);
163 if(ix
>=0x4b000000) /* |x|>=2**23, must be -integer */
164 return fabsf (x
)/zero
;
165 if (ix
> 0x40000000 /* X < 2.0f. */
166 && ix
< 0x41700000 /* X > -15.0f. */)
167 return __lgamma_negf (x
, signgamp
);
169 if(t
==zero
) return one
/fabsf(t
); /* -integer */
170 nadj
= __ieee754_logf(pi
/fabsf(t
*x
));
171 if(t
<zero
) *signgamp
= -1;
175 /* purge off 1 and 2 */
176 if (ix
==0x3f800000||ix
==0x40000000) r
= 0;
178 else if(ix
<0x40000000) {
179 if(ix
<=0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */
180 r
= -__ieee754_logf(x
);
181 if(ix
>=0x3f3b4a20) {y
= one
-x
; i
= 0;}
182 else if(ix
>=0x3e6d3308) {y
= x
-(tc
-one
); i
=1;}
186 if(ix
>=0x3fdda618) {y
=(float)2.0-x
;i
=0;} /* [1.7316,2] */
187 else if(ix
>=0x3F9da620) {y
=x
-tc
;i
=1;} /* [1.23,1.73] */
193 p1
= a0
+z
*(a2
+z
*(a4
+z
*(a6
+z
*(a8
+z
*a10
))));
194 p2
= z
*(a1
+z
*(a3
+z
*(a5
+z
*(a7
+z
*(a9
+z
*a11
)))));
196 r
+= (p
-(float)0.5*y
); break;
200 p1
= t0
+w
*(t3
+w
*(t6
+w
*(t9
+w
*t12
))); /* parallel comp */
201 p2
= t1
+w
*(t4
+w
*(t7
+w
*(t10
+w
*t13
)));
202 p3
= t2
+w
*(t5
+w
*(t8
+w
*(t11
+w
*t14
)));
203 p
= z
*p1
-(tt
-w
*(p2
+y
*p3
));
204 r
+= (tf
+ p
); break;
206 p1
= y
*(u0
+y
*(u1
+y
*(u2
+y
*(u3
+y
*(u4
+y
*u5
)))));
207 p2
= one
+y
*(v1
+y
*(v2
+y
*(v3
+y
*(v4
+y
*v5
))));
208 r
+= (-(float)0.5*y
+ p1
/p2
);
211 else if(ix
<0x41000000) { /* x < 8.0 */
215 p
= y
*(s0
+y
*(s1
+y
*(s2
+y
*(s3
+y
*(s4
+y
*(s5
+y
*s6
))))));
216 q
= one
+y
*(r1
+y
*(r2
+y
*(r3
+y
*(r4
+y
*(r5
+y
*r6
)))));
218 z
= one
; /* lgamma(1+s) = log(s) + lgamma(s) */
220 case 7: z
*= (y
+(float)6.0); /* FALLTHRU */
221 case 6: z
*= (y
+(float)5.0); /* FALLTHRU */
222 case 5: z
*= (y
+(float)4.0); /* FALLTHRU */
223 case 4: z
*= (y
+(float)3.0); /* FALLTHRU */
224 case 3: z
*= (y
+(float)2.0); /* FALLTHRU */
225 r
+= __ieee754_logf(z
); break;
227 /* 8.0 <= x < 2**26 */
228 } else if (ix
< 0x4c800000) {
229 t
= __ieee754_logf(x
);
232 w
= w0
+z
*(w1
+y
*(w2
+y
*(w3
+y
*(w4
+y
*(w5
+y
*w6
)))));
233 r
= (x
-half
)*(t
-one
)+w
;
235 /* 2**26 <= x <= inf */
236 r
= math_narrow_eval (x
*(__ieee754_logf(x
)-one
));
237 /* NADJ is set for negative arguments but not otherwise,
238 resulting in warnings that it may be used uninitialized
239 although in the cases where it is used it has always been
241 DIAG_PUSH_NEEDS_COMMENT
;
242 DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wmaybe-uninitialized");
243 if(hx
<0) r
= nadj
- r
;
244 DIAG_POP_NEEDS_COMMENT
;
247 strong_alias (__ieee754_lgammaf_r
, __lgammaf_r_finite
)