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