1 /* e_lgammaf_r.c -- float version of e_lgamma_r.c.
2 */
3
4 /*
5 * ====================================================
6 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
7 *
8 * Developed at SunPro, a Sun Microsystems, Inc. business.
9 * Permission to use, copy, modify, and distribute this
10 * software is freely granted, provided that this notice
11 * is preserved.
12 * ====================================================
13 */
14
15 #include <math.h>
16 #include <math-narrow-eval.h>
17 #include <math_private.h>
18 #include <libc-diag.h>
19 #include <libm-alias-finite.h>
20
21 static const float
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 */
88
89 static const float zero= 0.0000000000e+00;
90
91 static float
92 sin_pif(float x)
93 {
94 float y,z;
95 int n,ix;
96
97 GET_FLOAT_WORD(ix,x);
98 ix &= 0x7fffffff;
99
100 if(ix<0x3e800000) return __sinf (pi*x);
101 y = -x; /* x is assume negative */
102
103 /*
104 * argument reduction, make sure inexact flag not raised if input
105 * is an integer
106 */
107 z = floorf(y);
108 if(z!=y) { /* inexact anyway */
109 y *= (float)0.5;
110 y = (float)2.0*(y - floorf(y)); /* y = |x| mod 2.0 */
111 n = (int) (y*(float)4.0);
112 } else {
113 if(ix>=0x4b800000) {
114 y = zero; n = 0; /* y must be even */
115 } else {
116 if(ix<0x4b000000) z = y+two23; /* exact */
117 GET_FLOAT_WORD(n,z);
118 n &= 1;
119 y = n;
120 n<<= 2;
121 }
122 }
123 switch (n) {
124 case 0: y = __sinf (pi*y); break;
125 case 1:
126 case 2: y = __cosf (pi*((float)0.5-y)); break;
127 case 3:
128 case 4: y = __sinf (pi*(one-y)); break;
129 case 5:
130 case 6: y = -__cosf (pi*(y-(float)1.5)); break;
131 default: y = __sinf (pi*(y-(float)2.0)); break;
132 }
133 return -y;
134 }
135
136
137 float
138 __ieee754_lgammaf_r(float x, int *signgamp)
139 {
140 float t,y,z,nadj,p,p1,p2,p3,q,r,w;
141 int i,hx,ix;
142
143 GET_FLOAT_WORD(hx,x);
144
145 /* purge off +-inf, NaN, +-0, and negative arguments */
146 *signgamp = 1;
147 ix = hx&0x7fffffff;
148 if(__builtin_expect(ix>=0x7f800000, 0)) return x*x;
149 if(__builtin_expect(ix==0, 0))
150 {
151 if (hx < 0)
152 *signgamp = -1;
153 return one/fabsf(x);
154 }
155 if(__builtin_expect(ix<0x30800000, 0)) {
156 /* |x|<2**-30, return -log(|x|) */
157 if(hx<0) {
158 *signgamp = -1;
159 return -__ieee754_logf(-x);
160 } else return -__ieee754_logf(x);
161 }
162 if(hx<0) {
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);
168 t = sin_pif(x);
169 if(t==zero) return one/fabsf(t); /* -integer */
170 nadj = __ieee754_logf(pi/fabsf(t*x));
171 if(t<zero) *signgamp = -1;
172 x = -x;
173 }
174
175 /* purge off 1 and 2 */
176 if (ix==0x3f800000||ix==0x40000000) r = 0;
177 /* for x < 2.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;}
183 else {y = x; i=2;}
184 } else {
185 r = zero;
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] */
188 else {y=x-one;i=2;}
189 }
190 switch(i) {
191 case 0:
192 z = y*y;
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)))));
195 p = y*p1+p2;
196 r += (p-(float)0.5*y); break;
197 case 1:
198 z = y*y;
199 w = z*y;
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;
205 case 2:
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);
209 }
210 }
211 else if(ix<0x41000000) { /* x < 8.0 */
212 i = (int)x;
213 t = zero;
214 y = x-(float)i;
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)))));
217 r = half*y+p/q;
218 z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
219 switch(i) {
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;
226 }
227 /* 8.0 <= x < 2**26 */
228 } else if (ix < 0x4c800000) {
229 t = __ieee754_logf(x);
230 z = one/x;
231 y = z*z;
232 w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
233 r = (x-half)*(t-one)+w;
234 } else
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
240 set. */
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;
245 return r;
246 }
247 libm_alias_finite (__ieee754_lgammaf_r, __lgammaf_r)