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