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 static char rcsid
[] = "$FreeBSD: src/lib/msun/src/e_lgammaf_r.c,v 1.10 2007/05/02 16:54:22 bde Exp $";
21 #include "math_private.h"
24 two23
= 8.3886080000e+06, /* 0x4b000000 */
25 half
= 5.0000000000e-01, /* 0x3f000000 */
26 one
= 1.0000000000e+00, /* 0x3f800000 */
27 pi
= 3.1415927410e+00, /* 0x40490fdb */
28 a0
= 7.7215664089e-02, /* 0x3d9e233f */
29 a1
= 3.2246702909e-01, /* 0x3ea51a66 */
30 a2
= 6.7352302372e-02, /* 0x3d89f001 */
31 a3
= 2.0580807701e-02, /* 0x3ca89915 */
32 a4
= 7.3855509982e-03, /* 0x3bf2027e */
33 a5
= 2.8905137442e-03, /* 0x3b3d6ec6 */
34 a6
= 1.1927076848e-03, /* 0x3a9c54a1 */
35 a7
= 5.1006977446e-04, /* 0x3a05b634 */
36 a8
= 2.2086278477e-04, /* 0x39679767 */
37 a9
= 1.0801156895e-04, /* 0x38e28445 */
38 a10
= 2.5214456400e-05, /* 0x37d383a2 */
39 a11
= 4.4864096708e-05, /* 0x383c2c75 */
40 tc
= 1.4616321325e+00, /* 0x3fbb16c3 */
41 tf
= -1.2148628384e-01, /* 0xbdf8cdcd */
42 /* tt = -(tail of tf) */
43 tt
= 6.6971006518e-09, /* 0x31e61c52 */
44 t0
= 4.8383611441e-01, /* 0x3ef7b95e */
45 t1
= -1.4758771658e-01, /* 0xbe17213c */
46 t2
= 6.4624942839e-02, /* 0x3d845a15 */
47 t3
= -3.2788541168e-02, /* 0xbd064d47 */
48 t4
= 1.7970675603e-02, /* 0x3c93373d */
49 t5
= -1.0314224288e-02, /* 0xbc28fcfe */
50 t6
= 6.1005386524e-03, /* 0x3bc7e707 */
51 t7
= -3.6845202558e-03, /* 0xbb7177fe */
52 t8
= 2.2596477065e-03, /* 0x3b141699 */
53 t9
= -1.4034647029e-03, /* 0xbab7f476 */
54 t10
= 8.8108185446e-04, /* 0x3a66f867 */
55 t11
= -5.3859531181e-04, /* 0xba0d3085 */
56 t12
= 3.1563205994e-04, /* 0x39a57b6b */
57 t13
= -3.1275415677e-04, /* 0xb9a3f927 */
58 t14
= 3.3552918467e-04, /* 0x39afe9f7 */
59 u0
= -7.7215664089e-02, /* 0xbd9e233f */
60 u1
= 6.3282704353e-01, /* 0x3f2200f4 */
61 u2
= 1.4549225569e+00, /* 0x3fba3ae7 */
62 u3
= 9.7771751881e-01, /* 0x3f7a4bb2 */
63 u4
= 2.2896373272e-01, /* 0x3e6a7578 */
64 u5
= 1.3381091878e-02, /* 0x3c5b3c5e */
65 v1
= 2.4559779167e+00, /* 0x401d2ebe */
66 v2
= 2.1284897327e+00, /* 0x4008392d */
67 v3
= 7.6928514242e-01, /* 0x3f44efdf */
68 v4
= 1.0422264785e-01, /* 0x3dd572af */
69 v5
= 3.2170924824e-03, /* 0x3b52d5db */
70 s0
= -7.7215664089e-02, /* 0xbd9e233f */
71 s1
= 2.1498242021e-01, /* 0x3e5c245a */
72 s2
= 3.2577878237e-01, /* 0x3ea6cc7a */
73 s3
= 1.4635047317e-01, /* 0x3e15dce6 */
74 s4
= 2.6642270386e-02, /* 0x3cda40e4 */
75 s5
= 1.8402845599e-03, /* 0x3af135b4 */
76 s6
= 3.1947532989e-05, /* 0x3805ff67 */
77 r1
= 1.3920053244e+00, /* 0x3fb22d3b */
78 r2
= 7.2193557024e-01, /* 0x3f38d0c5 */
79 r3
= 1.7193385959e-01, /* 0x3e300f6e */
80 r4
= 1.8645919859e-02, /* 0x3c98bf54 */
81 r5
= 7.7794247773e-04, /* 0x3a4beed6 */
82 r6
= 7.3266842264e-06, /* 0x36f5d7bd */
83 w0
= 4.1893854737e-01, /* 0x3ed67f1d */
84 w1
= 8.3333335817e-02, /* 0x3daaaaab */
85 w2
= -2.7777778450e-03, /* 0xbb360b61 */
86 w3
= 7.9365057172e-04, /* 0x3a500cfd */
87 w4
= -5.9518753551e-04, /* 0xba1c065c */
88 w5
= 8.3633989561e-04, /* 0x3a5b3dd2 */
89 w6
= -1.6309292987e-03; /* 0xbad5c4e8 */
91 static const float zero
= 0.0000000000e+00;
93 static float sin_pif(float x
)
101 if(ix
<0x3e800000) return __kernel_sindf(pi
*x
);
102 y
= -x
; /* x is assume negative */
105 * argument reduction, make sure inexact flag not raised if input
109 if(z
!=y
) { /* inexact anyway */
111 y
= (float)2.0*(y
- floorf(y
)); /* y = |x| mod 2.0 */
112 n
= (int) (y
*(float)4.0);
115 y
= zero
; n
= 0; /* y must be even */
117 if(ix
<0x4b000000) z
= y
+two23
; /* exact */
125 case 0: y
= __kernel_sindf(pi
*y
); break;
127 case 2: y
= __kernel_cosdf(pi
*((float)0.5-y
)); break;
129 case 4: y
= __kernel_sindf(pi
*(one
-y
)); break;
131 case 6: y
= -__kernel_cosdf(pi
*(y
-(float)1.5)); break;
132 default: y
= __kernel_sindf(pi
*(y
-(float)2.0)); break;
139 __ieee754_lgammaf_r(float x
, int *signgamp
)
141 float t
,y
,z
,nadj
,p
,p1
,p2
,p3
,q
,r
,w
;
145 GET_FLOAT_WORD(hx
,x
);
147 /* purge off +-inf, NaN, +-0, tiny and negative arguments */
150 if(ix
>=0x7f800000) return x
*x
;
151 if(ix
==0) return one
/zero
;
152 if(ix
<0x35000000) { /* |x|<2**-21, return -log(|x|) */
155 return -__ieee754_logf(-x
);
156 } else return -__ieee754_logf(x
);
159 if(ix
>=0x4b000000) /* |x|>=2**23, must be -integer */
162 if(t
==zero
) return one
/zero
; /* -integer */
163 nadj
= __ieee754_logf(pi
/fabsf(t
*x
));
164 if(t
<zero
) *signgamp
= -1;
168 /* purge off 1 and 2 */
169 if (ix
==0x3f800000||ix
==0x40000000) r
= 0;
171 else if(ix
<0x40000000) {
172 if(ix
<=0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */
173 r
= -__ieee754_logf(x
);
174 if(ix
>=0x3f3b4a20) {y
= one
-x
; i
= 0;}
175 else if(ix
>=0x3e6d3308) {y
= x
-(tc
-one
); i
=1;}
179 if(ix
>=0x3fdda618) {y
=(float)2.0-x
;i
=0;} /* [1.7316,2] */
180 else if(ix
>=0x3F9da620) {y
=x
-tc
;i
=1;} /* [1.23,1.73] */
186 p1
= a0
+z
*(a2
+z
*(a4
+z
*(a6
+z
*(a8
+z
*a10
))));
187 p2
= z
*(a1
+z
*(a3
+z
*(a5
+z
*(a7
+z
*(a9
+z
*a11
)))));
189 r
+= (p
-(float)0.5*y
); break;
193 p1
= t0
+w
*(t3
+w
*(t6
+w
*(t9
+w
*t12
))); /* parallel comp */
194 p2
= t1
+w
*(t4
+w
*(t7
+w
*(t10
+w
*t13
)));
195 p3
= t2
+w
*(t5
+w
*(t8
+w
*(t11
+w
*t14
)));
196 p
= z
*p1
-(tt
-w
*(p2
+y
*p3
));
197 r
+= (tf
+ p
); break;
199 p1
= y
*(u0
+y
*(u1
+y
*(u2
+y
*(u3
+y
*(u4
+y
*u5
)))));
200 p2
= one
+y
*(v1
+y
*(v2
+y
*(v3
+y
*(v4
+y
*v5
))));
201 r
+= (-(float)0.5*y
+ p1
/p2
);
204 else if(ix
<0x41000000) { /* x < 8.0 */
208 p
= y
*(s0
+y
*(s1
+y
*(s2
+y
*(s3
+y
*(s4
+y
*(s5
+y
*s6
))))));
209 q
= one
+y
*(r1
+y
*(r2
+y
*(r3
+y
*(r4
+y
*(r5
+y
*r6
)))));
211 z
= one
; /* lgamma(1+s) = log(s) + lgamma(s) */
213 case 7: z
*= (y
+(float)6.0); /* FALLTHRU */
214 case 6: z
*= (y
+(float)5.0); /* FALLTHRU */
215 case 5: z
*= (y
+(float)4.0); /* FALLTHRU */
216 case 4: z
*= (y
+(float)3.0); /* FALLTHRU */
217 case 3: z
*= (y
+(float)2.0); /* FALLTHRU */
218 r
+= __ieee754_logf(z
); break;
220 /* 8.0 <= x < 2**58 */
221 } else if (ix
< 0x5c800000) {
222 t
= __ieee754_logf(x
);
225 w
= w0
+z
*(w1
+y
*(w2
+y
*(w3
+y
*(w4
+y
*(w5
+y
*w6
)))));
226 r
= (x
-half
)*(t
-one
)+w
;
228 /* 2**58 <= x <= inf */
229 r
= x
*(__ieee754_logf(x
)-one
);
230 if(hx
<0) r
= nadj
- r
;