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