diff options
author | Elliott Hughes | 2014-10-09 18:20:37 -0500 |
---|---|---|
committer | Elliott Hughes | 2014-10-09 18:20:37 -0500 |
commit | 488268b134723c7a6598338bb253be5f64d53be4 (patch) | |
tree | e2e3adbcb193d4f052907f7be991621469b5dd7c /libm | |
parent | e9c216fca56e84b0d0a96f7d5e3c99d3276ef071 (diff) | |
download | platform-bionic-488268b134723c7a6598338bb253be5f64d53be4.tar.gz platform-bionic-488268b134723c7a6598338bb253be5f64d53be4.tar.xz platform-bionic-488268b134723c7a6598338bb253be5f64d53be4.zip |
Sync libm with upstream.
Change-Id: I3b4e2c9c6ce6c5934f270a51ce5eb9154c5805d5
Diffstat (limited to 'libm')
-rw-r--r-- | libm/upstream-freebsd/lib/msun/ld128/e_lgammal_r.c | 73 | ||||
-rw-r--r-- | libm/upstream-freebsd/lib/msun/src/e_acoshl.c | 4 | ||||
-rw-r--r-- | libm/upstream-freebsd/lib/msun/src/e_atanhl.c | 4 | ||||
-rw-r--r-- | libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c | 77 | ||||
-rw-r--r-- | libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c | 43 | ||||
-rw-r--r-- | libm/upstream-freebsd/lib/msun/src/fenv-softfloat.h | 184 | ||||
-rw-r--r-- | libm/upstream-freebsd/lib/msun/src/s_remquol.c | 2 |
7 files changed, 99 insertions, 288 deletions
diff --git a/libm/upstream-freebsd/lib/msun/ld128/e_lgammal_r.c b/libm/upstream-freebsd/lib/msun/ld128/e_lgammal_r.c index 13f8f75a..53d3af17 100644 --- a/libm/upstream-freebsd/lib/msun/ld128/e_lgammal_r.c +++ b/libm/upstream-freebsd/lib/msun/ld128/e_lgammal_r.c | |||
@@ -206,13 +206,13 @@ sin_pil(long double x) | |||
206 | n--; | 206 | n--; |
207 | } | 207 | } |
208 | n &= 7; | 208 | n &= 7; |
209 | y = y - z + n * 0.25L; | 209 | y = y - z + n * 0.25; |
210 | 210 | ||
211 | switch (n) { | 211 | switch (n) { |
212 | case 0: y = __kernel_sinl(pi*y,zero,0); break; | 212 | case 0: y = __kernel_sinl(pi*y,zero,0); break; |
213 | case 1: | 213 | case 1: |
214 | case 2: y = __kernel_cosl(pi*(0.5-y),zero); break; | 214 | case 2: y = __kernel_cosl(pi*(0.5-y),zero); break; |
215 | case 3: | 215 | case 3: |
216 | case 4: y = __kernel_sinl(pi*(one-y),zero,0); break; | 216 | case 4: y = __kernel_sinl(pi*(one-y),zero,0); break; |
217 | case 5: | 217 | case 5: |
218 | case 6: y = -__kernel_cosl(pi*(y-1.5),zero); break; | 218 | case 6: y = -__kernel_cosl(pi*(y-1.5),zero); break; |
@@ -221,41 +221,33 @@ sin_pil(long double x) | |||
221 | return -y; | 221 | return -y; |
222 | } | 222 | } |
223 | 223 | ||
224 | |||
225 | long double | 224 | long double |
226 | lgammal_r(long double x, int *signgamp) | 225 | lgammal_r(long double x, int *signgamp) |
227 | { | 226 | { |
228 | long double nadj,p,p1,p2,p3,q,r,t,w,y,z; | 227 | long double nadj,p,p1,p2,p3,q,r,t,w,y,z; |
229 | uint64_t llx,lx; | 228 | uint64_t llx,lx; |
230 | int i; | 229 | int i; |
231 | uint16_t hx; | 230 | uint16_t hx,ix; |
232 | 231 | ||
233 | EXTRACT_LDBL128_WORDS(hx, lx, llx, x); | 232 | EXTRACT_LDBL128_WORDS(hx,lx,llx,x); |
234 | |||
235 | if((hx & 0x7fff) == 0x7fff) { /* erfl(nan)=nan */ | ||
236 | i = (hx>>15)<<1; | ||
237 | return (1-i)+one/x; /* erfl(+-inf)=+-1 */ | ||
238 | } | ||
239 | 233 | ||
240 | /* purge off +-inf, NaN, +-0, tiny and negative arguments */ | 234 | /* purge +-Inf and NaNs */ |
241 | *signgamp = 1; | 235 | *signgamp = 1; |
242 | if((hx & 0x7fff) == 0x7fff) /* x is +-Inf or NaN */ | 236 | ix = hx&0x7fff; |
243 | return x*x; | 237 | if(ix==0x7fff) return x*x; |
244 | if((hx==0||hx==0x8000)&&lx==0) { | ||
245 | if (hx&0x8000) | ||
246 | *signgamp = -1; | ||
247 | return one/vzero; | ||
248 | } | ||
249 | 238 | ||
250 | /* purge off tiny and negative arguments */ | 239 | /* purge +-0 and tiny arguments */ |
251 | if(fabsl(x)<0x1p-119L) { | 240 | *signgamp = 1-2*(hx>>15); |
252 | if(hx&0x8000) { | 241 | if(ix<0x3fff-116) { /* |x|<2**-(p+3), return -log(|x|) */ |
253 | *signgamp = -1; | 242 | if((ix|lx|llx)==0) |
254 | return -logl(-x); | 243 | return one/vzero; |
255 | } else return -logl(x); | 244 | return -logl(fabsl(x)); |
256 | } | 245 | } |
246 | |||
247 | /* purge negative integers and start evaluation for other x < 0 */ | ||
257 | if(hx&0x8000) { | 248 | if(hx&0x8000) { |
258 | if(fabsl(x)>=0x1p112) | 249 | *signgamp = 1; |
250 | if(ix>=0x3fff+112) /* |x|>=2**(p-1), must be -integer */ | ||
259 | return one/vzero; | 251 | return one/vzero; |
260 | t = sin_pil(x); | 252 | t = sin_pil(x); |
261 | if(t==zero) return one/vzero; | 253 | if(t==zero) return one/vzero; |
@@ -264,17 +256,19 @@ lgammal_r(long double x, int *signgamp) | |||
264 | x = -x; | 256 | x = -x; |
265 | } | 257 | } |
266 | 258 | ||
267 | if(x == 1 || x ==2) r = 0; | 259 | /* purge 1 and 2 */ |
268 | else if(x<2) { | 260 | if((ix==0x3fff || ix==0x4000) && (lx|llx)==0) r = 0; |
269 | if(x<=0.8999996185302734) { | 261 | /* for x < 2.0 */ |
262 | else if(ix<0x4000) { | ||
263 | if(x<=8.9999961853027344e-01) { | ||
270 | r = -logl(x); | 264 | r = -logl(x); |
271 | if(x>=0.7315998077392578) {y = 1-x; i= 0;} | 265 | if(x>=7.3159980773925781e-01) {y = 1-x; i= 0;} |
272 | else if(x>=0.2316399812698364) {y= x-(tc-1); i=1;} | 266 | else if(x>=2.3163998126983643e-01) {y= x-(tc-1); i=1;} |
273 | else {y = x; i=2;} | 267 | else {y = x; i=2;} |
274 | } else { | 268 | } else { |
275 | r = 0; | 269 | r = 0; |
276 | if(x>=1.7316312789916992) {y=2-x;i=0;} | 270 | if(x>=1.7316312789916992e+00) {y=2-x;i=0;} |
277 | else if(x>=1.2316322326660156) {y=x-tc;i=1;} | 271 | else if(x>=1.2316322326660156e+00) {y=x-tc;i=1;} |
278 | else {y=x-1;i=2;} | 272 | else {y=x-1;i=2;} |
279 | } | 273 | } |
280 | switch(i) { | 274 | switch(i) { |
@@ -285,23 +279,24 @@ lgammal_r(long double x, int *signgamp) | |||
285 | p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*(a13+z*(a15+ | 279 | p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*(a13+z*(a15+ |
286 | z*(a17+z*(a19+z*(a21+z*a23))))))))))); | 280 | z*(a17+z*(a19+z*(a21+z*a23))))))))))); |
287 | p = y*p1+p2; | 281 | p = y*p1+p2; |
288 | r += (p-y/2); break; | 282 | r += p-y/2; break; |
289 | case 1: | 283 | case 1: |
290 | p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+ | 284 | p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+ |
291 | y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+ | 285 | y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+ |
292 | y*(t17+y*(t18+y*(t19+y*(t20+y*(t21+y*(t22+y*(t23+ | 286 | y*(t17+y*(t18+y*(t19+y*(t20+y*(t21+y*(t22+y*(t23+ |
293 | y*(t24+y*(t25+y*(t26+y*(t27+y*(t28+y*(t29+y*(t30+ | 287 | y*(t24+y*(t25+y*(t26+y*(t27+y*(t28+y*(t29+y*(t30+ |
294 | y*(t31+y*t32)))))))))))))))))))))))))))))); | 288 | y*(t31+y*t32)))))))))))))))))))))))))))))); |
295 | r += (tf + p); break; | 289 | r += tf + p; break; |
296 | case 2: | 290 | case 2: |
297 | p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*(u6+y*(u7+ | 291 | p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*(u6+y*(u7+ |
298 | y*(u8+y*(u9+y*u10)))))))))); | 292 | y*(u8+y*(u9+y*u10)))))))))); |
299 | p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*(v6+y*(v7+ | 293 | p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*(v6+y*(v7+ |
300 | y*(v8+y*(v9+y*(v10+y*v11)))))))))); | 294 | y*(v8+y*(v9+y*(v10+y*v11)))))))))); |
301 | r += (-y/2 + p1/p2); | 295 | r += p1/p2-y/2; |
302 | } | 296 | } |
303 | } | 297 | } |
304 | else if(x<8) { | 298 | /* x < 8.0 */ |
299 | else if(ix<0x4002) { | ||
305 | i = x; | 300 | i = x; |
306 | y = x-i; | 301 | y = x-i; |
307 | p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*(s6+y*(s7+y*(s8+ | 302 | p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*(s6+y*(s7+y*(s8+ |
@@ -318,7 +313,8 @@ lgammal_r(long double x, int *signgamp) | |||
318 | case 3: z *= (y+2); /* FALLTHRU */ | 313 | case 3: z *= (y+2); /* FALLTHRU */ |
319 | r += logl(z); break; | 314 | r += logl(z); break; |
320 | } | 315 | } |
321 | } else if (x < 0x1p119L) { | 316 | /* 8.0 <= x < 2**(p+3) */ |
317 | } else if (ix<0x3fff+116) { | ||
322 | t = logl(x); | 318 | t = logl(x); |
323 | z = one/x; | 319 | z = one/x; |
324 | y = z*z; | 320 | y = z*z; |
@@ -326,6 +322,7 @@ lgammal_r(long double x, int *signgamp) | |||
326 | y*(w9+y*(w10+y*(w11+y*(w12+y*(w13+y*(w14+y*(w15+y*(w16+ | 322 | y*(w9+y*(w10+y*(w11+y*(w12+y*(w13+y*(w14+y*(w15+y*(w16+ |
327 | y*(w17+y*w18))))))))))))))))); | 323 | y*(w17+y*w18))))))))))))))))); |
328 | r = (x-half)*(t-one)+w; | 324 | r = (x-half)*(t-one)+w; |
325 | /* 2**(p+3) <= x <= inf */ | ||
329 | } else | 326 | } else |
330 | r = x*(logl(x)-1); | 327 | r = x*(logl(x)-1); |
331 | if(hx&0x8000) r = nadj - r; | 328 | if(hx&0x8000) r = nadj - r; |
diff --git a/libm/upstream-freebsd/lib/msun/src/e_acoshl.c b/libm/upstream-freebsd/lib/msun/src/e_acoshl.c index 59faeb0f..b9f3aed6 100644 --- a/libm/upstream-freebsd/lib/msun/src/e_acoshl.c +++ b/libm/upstream-freebsd/lib/msun/src/e_acoshl.c | |||
@@ -7,7 +7,7 @@ | |||
7 | * | 7 | * |
8 | * Developed at SunSoft, a Sun Microsystems, Inc. business. | 8 | * Developed at SunSoft, a Sun Microsystems, Inc. business. |
9 | * Permission to use, copy, modify, and distribute this | 9 | * Permission to use, copy, modify, and distribute this |
10 | * software is freely granted, provided that this notice | 10 | * software is freely granted, provided that this notice |
11 | * is preserved. | 11 | * is preserved. |
12 | * ==================================================== | 12 | * ==================================================== |
13 | * | 13 | * |
@@ -75,7 +75,7 @@ acoshl(long double x) | |||
75 | } else if (hx >= BIAS + EXP_LARGE) { /* x >= LARGE */ | 75 | } else if (hx >= BIAS + EXP_LARGE) { /* x >= LARGE */ |
76 | if (hx >= 0x7fff) { /* x is inf, NaN or misnormal */ | 76 | if (hx >= 0x7fff) { /* x is inf, NaN or misnormal */ |
77 | RETURNI(x+x); | 77 | RETURNI(x+x); |
78 | } else | 78 | } else |
79 | RETURNI(logl(x)+ln2); /* acosh(huge)=log(2x), or misnormal */ | 79 | RETURNI(logl(x)+ln2); /* acosh(huge)=log(2x), or misnormal */ |
80 | } else if (hx == 0x3fff && x == 1) { | 80 | } else if (hx == 0x3fff && x == 1) { |
81 | RETURNI(0.0); /* acosh(1) = 0 */ | 81 | RETURNI(0.0); /* acosh(1) = 0 */ |
diff --git a/libm/upstream-freebsd/lib/msun/src/e_atanhl.c b/libm/upstream-freebsd/lib/msun/src/e_atanhl.c index a888426d..11d56ea5 100644 --- a/libm/upstream-freebsd/lib/msun/src/e_atanhl.c +++ b/libm/upstream-freebsd/lib/msun/src/e_atanhl.c | |||
@@ -7,7 +7,7 @@ | |||
7 | * | 7 | * |
8 | * Developed at SunSoft, a Sun Microsystems, Inc. business. | 8 | * Developed at SunSoft, a Sun Microsystems, Inc. business. |
9 | * Permission to use, copy, modify, and distribute this | 9 | * Permission to use, copy, modify, and distribute this |
10 | * software is freely granted, provided that this notice | 10 | * software is freely granted, provided that this notice |
11 | * is preserved. | 11 | * is preserved. |
12 | * ==================================================== | 12 | * ==================================================== |
13 | * | 13 | * |
@@ -68,7 +68,7 @@ atanhl(long double x) | |||
68 | if (ix < 0x3ffe) { /* |x| < 0.5, or misnormal */ | 68 | if (ix < 0x3ffe) { /* |x| < 0.5, or misnormal */ |
69 | t = x+x; | 69 | t = x+x; |
70 | t = 0.5*log1pl(t+t*x/(one-x)); | 70 | t = 0.5*log1pl(t+t*x/(one-x)); |
71 | } else | 71 | } else |
72 | t = 0.5*log1pl((x+x)/(one-x)); | 72 | t = 0.5*log1pl((x+x)/(one-x)); |
73 | RETURNI((hx & 0x8000) == 0 ? t : -t); | 73 | RETURNI((hx & 0x8000) == 0 ? t : -t); |
74 | } | 74 | } |
diff --git a/libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c b/libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c index a6fd691f..be70767e 100644 --- a/libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c +++ b/libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c | |||
@@ -1,4 +1,3 @@ | |||
1 | |||
2 | /* @(#)e_lgamma_r.c 1.3 95/01/18 */ | 1 | /* @(#)e_lgamma_r.c 1.3 95/01/18 */ |
3 | /* | 2 | /* |
4 | * ==================================================== | 3 | * ==================================================== |
@@ -6,22 +5,21 @@ | |||
6 | * | 5 | * |
7 | * Developed at SunSoft, a Sun Microsystems, Inc. business. | 6 | * Developed at SunSoft, a Sun Microsystems, Inc. business. |
8 | * Permission to use, copy, modify, and distribute this | 7 | * Permission to use, copy, modify, and distribute this |
9 | * software is freely granted, provided that this notice | 8 | * software is freely granted, provided that this notice |
10 | * is preserved. | 9 | * is preserved. |
11 | * ==================================================== | 10 | * ==================================================== |
12 | * | ||
13 | */ | 11 | */ |
14 | 12 | ||
15 | #include <sys/cdefs.h> | 13 | #include <sys/cdefs.h> |
16 | __FBSDID("$FreeBSD$"); | 14 | __FBSDID("$FreeBSD$"); |
17 | 15 | ||
18 | /* __ieee754_lgamma_r(x, signgamp) | 16 | /* __ieee754_lgamma_r(x, signgamp) |
19 | * Reentrant version of the logarithm of the Gamma function | 17 | * Reentrant version of the logarithm of the Gamma function |
20 | * with user provide pointer for the sign of Gamma(x). | 18 | * with user provide pointer for the sign of Gamma(x). |
21 | * | 19 | * |
22 | * Method: | 20 | * Method: |
23 | * 1. Argument Reduction for 0 < x <= 8 | 21 | * 1. Argument Reduction for 0 < x <= 8 |
24 | * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may | 22 | * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may |
25 | * reduce x to a number in [1.5,2.5] by | 23 | * reduce x to a number in [1.5,2.5] by |
26 | * lgamma(1+s) = log(s) + lgamma(s) | 24 | * lgamma(1+s) = log(s) + lgamma(s) |
27 | * for example, | 25 | * for example, |
@@ -59,20 +57,20 @@ __FBSDID("$FreeBSD$"); | |||
59 | * by | 57 | * by |
60 | * 3 5 11 | 58 | * 3 5 11 |
61 | * w = w0 + w1*z + w2*z + w3*z + ... + w6*z | 59 | * w = w0 + w1*z + w2*z + w3*z + ... + w6*z |
62 | * where | 60 | * where |
63 | * |w - f(z)| < 2**-58.74 | 61 | * |w - f(z)| < 2**-58.74 |
64 | * | 62 | * |
65 | * 4. For negative x, since (G is gamma function) | 63 | * 4. For negative x, since (G is gamma function) |
66 | * -x*G(-x)*G(x) = pi/sin(pi*x), | 64 | * -x*G(-x)*G(x) = pi/sin(pi*x), |
67 | * we have | 65 | * we have |
68 | * G(x) = pi/(sin(pi*x)*(-x)*G(-x)) | 66 | * G(x) = pi/(sin(pi*x)*(-x)*G(-x)) |
69 | * since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0 | 67 | * since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0 |
70 | * Hence, for x<0, signgam = sign(sin(pi*x)) and | 68 | * Hence, for x<0, signgam = sign(sin(pi*x)) and |
71 | * lgamma(x) = log(|Gamma(x)|) | 69 | * lgamma(x) = log(|Gamma(x)|) |
72 | * = log(pi/(|x*sin(pi*x)|)) - lgamma(-x); | 70 | * = log(pi/(|x*sin(pi*x)|)) - lgamma(-x); |
73 | * Note: one should avoid compute pi*(-x) directly in the | 71 | * Note: one should avoid compute pi*(-x) directly in the |
74 | * computation of sin(pi*(-x)). | 72 | * computation of sin(pi*(-x)). |
75 | * | 73 | * |
76 | * 5. Special Cases | 74 | * 5. Special Cases |
77 | * lgamma(2+s) ~ s*(1-Euler) for tiny s | 75 | * lgamma(2+s) ~ s*(1-Euler) for tiny s |
78 | * lgamma(1) = lgamma(2) = 0 | 76 | * lgamma(1) = lgamma(2) = 0 |
@@ -80,7 +78,6 @@ __FBSDID("$FreeBSD$"); | |||
80 | * lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero | 78 | * lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero |
81 | * lgamma(inf) = inf | 79 | * lgamma(inf) = inf |
82 | * lgamma(-inf) = inf (bug for bug compatible with C99!?) | 80 | * lgamma(-inf) = inf (bug for bug compatible with C99!?) |
83 | * | ||
84 | */ | 81 | */ |
85 | 82 | ||
86 | #include <float.h> | 83 | #include <float.h> |
@@ -189,9 +186,9 @@ sin_pi(double x) | |||
189 | 186 | ||
190 | switch (n) { | 187 | switch (n) { |
191 | case 0: y = __kernel_sin(pi*y,zero,0); break; | 188 | case 0: y = __kernel_sin(pi*y,zero,0); break; |
192 | case 1: | 189 | case 1: |
193 | case 2: y = __kernel_cos(pi*(0.5-y),zero); break; | 190 | case 2: y = __kernel_cos(pi*(0.5-y),zero); break; |
194 | case 3: | 191 | case 3: |
195 | case 4: y = __kernel_sin(pi*(one-y),zero,0); break; | 192 | case 4: y = __kernel_sin(pi*(one-y),zero,0); break; |
196 | case 5: | 193 | case 5: |
197 | case 6: y = -__kernel_cos(pi*(y-1.5),zero); break; | 194 | case 6: y = -__kernel_cos(pi*(y-1.5),zero); break; |
@@ -204,28 +201,28 @@ sin_pi(double x) | |||
204 | double | 201 | double |
205 | __ieee754_lgamma_r(double x, int *signgamp) | 202 | __ieee754_lgamma_r(double x, int *signgamp) |
206 | { | 203 | { |
207 | double t,y,z,nadj,p,p1,p2,p3,q,r,w; | 204 | double nadj,p,p1,p2,p3,q,r,t,w,y,z; |
208 | int32_t hx; | 205 | int32_t hx; |
209 | int i,ix,lx; | 206 | int i,ix,lx; |
210 | 207 | ||
211 | EXTRACT_WORDS(hx,lx,x); | 208 | EXTRACT_WORDS(hx,lx,x); |
212 | 209 | ||
213 | /* purge off +-inf, NaN, +-0, tiny and negative arguments */ | 210 | /* purge +-Inf and NaNs */ |
214 | *signgamp = 1; | 211 | *signgamp = 1; |
215 | ix = hx&0x7fffffff; | 212 | ix = hx&0x7fffffff; |
216 | if(ix>=0x7ff00000) return x*x; | 213 | if(ix>=0x7ff00000) return x*x; |
217 | if((ix|lx)==0) { | 214 | |
218 | if(hx<0) | 215 | /* purge +-0 and tiny arguments */ |
219 | *signgamp = -1; | 216 | *signgamp = 1-2*((uint32_t)hx>>31); |
220 | return one/vzero; | 217 | if(ix<0x3c700000) { /* |x|<2**-56, return -log(|x|) */ |
221 | } | 218 | if((ix|lx)==0) |
222 | if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */ | 219 | return one/vzero; |
223 | if(hx<0) { | 220 | return -__ieee754_log(fabs(x)); |
224 | *signgamp = -1; | ||
225 | return -__ieee754_log(-x); | ||
226 | } else return -__ieee754_log(x); | ||
227 | } | 221 | } |
222 | |||
223 | /* purge negative integers and start evaluation for other x < 0 */ | ||
228 | if(hx<0) { | 224 | if(hx<0) { |
225 | *signgamp = 1; | ||
229 | if(ix>=0x43300000) /* |x|>=2**52, must be -integer */ | 226 | if(ix>=0x43300000) /* |x|>=2**52, must be -integer */ |
230 | return one/vzero; | 227 | return one/vzero; |
231 | t = sin_pi(x); | 228 | t = sin_pi(x); |
@@ -235,7 +232,7 @@ __ieee754_lgamma_r(double x, int *signgamp) | |||
235 | x = -x; | 232 | x = -x; |
236 | } | 233 | } |
237 | 234 | ||
238 | /* purge off 1 and 2 */ | 235 | /* purge 1 and 2 */ |
239 | if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0; | 236 | if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0; |
240 | /* for x < 2.0 */ | 237 | /* for x < 2.0 */ |
241 | else if(ix<0x40000000) { | 238 | else if(ix<0x40000000) { |
@@ -256,7 +253,7 @@ __ieee754_lgamma_r(double x, int *signgamp) | |||
256 | p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); | 253 | p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); |
257 | p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); | 254 | p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); |
258 | p = y*p1+p2; | 255 | p = y*p1+p2; |
259 | r += (p-y/2); break; | 256 | r += p-y/2; break; |
260 | case 1: | 257 | case 1: |
261 | z = y*y; | 258 | z = y*y; |
262 | w = z*y; | 259 | w = z*y; |
@@ -264,19 +261,20 @@ __ieee754_lgamma_r(double x, int *signgamp) | |||
264 | p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); | 261 | p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); |
265 | p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); | 262 | p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); |
266 | p = z*p1-(tt-w*(p2+y*p3)); | 263 | p = z*p1-(tt-w*(p2+y*p3)); |
267 | r += (tf + p); break; | 264 | r += tf + p; break; |
268 | case 2: | 265 | case 2: |
269 | p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); | 266 | p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); |
270 | p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); | 267 | p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); |
271 | r += (-0.5*y + p1/p2); | 268 | r += p1/p2-y/2; |
272 | } | 269 | } |
273 | } | 270 | } |
274 | else if(ix<0x40200000) { /* x < 8.0 */ | 271 | /* x < 8.0 */ |
275 | i = (int)x; | 272 | else if(ix<0x40200000) { |
276 | y = x-(double)i; | 273 | i = x; |
274 | y = x-i; | ||
277 | p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); | 275 | p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); |
278 | q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); | 276 | q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); |
279 | r = half*y+p/q; | 277 | r = y/2+p/q; |
280 | z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ | 278 | z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ |
281 | switch(i) { | 279 | switch(i) { |
282 | case 7: z *= (y+6); /* FALLTHRU */ | 280 | case 7: z *= (y+6); /* FALLTHRU */ |
@@ -286,15 +284,15 @@ __ieee754_lgamma_r(double x, int *signgamp) | |||
286 | case 3: z *= (y+2); /* FALLTHRU */ | 284 | case 3: z *= (y+2); /* FALLTHRU */ |
287 | r += __ieee754_log(z); break; | 285 | r += __ieee754_log(z); break; |
288 | } | 286 | } |
289 | /* 8.0 <= x < 2**58 */ | 287 | /* 8.0 <= x < 2**56 */ |
290 | } else if (ix < 0x43900000) { | 288 | } else if (ix < 0x43700000) { |
291 | t = __ieee754_log(x); | 289 | t = __ieee754_log(x); |
292 | z = one/x; | 290 | z = one/x; |
293 | y = z*z; | 291 | y = z*z; |
294 | w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); | 292 | w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); |
295 | r = (x-half)*(t-one)+w; | 293 | r = (x-half)*(t-one)+w; |
296 | } else | 294 | } else |
297 | /* 2**58 <= x <= inf */ | 295 | /* 2**56 <= x <= inf */ |
298 | r = x*(__ieee754_log(x)-one); | 296 | r = x*(__ieee754_log(x)-one); |
299 | if(hx<0) r = nadj - r; | 297 | if(hx<0) r = nadj - r; |
300 | return r; | 298 | return r; |
@@ -303,4 +301,3 @@ __ieee754_lgamma_r(double x, int *signgamp) | |||
303 | #if (LDBL_MANT_DIG == 53) | 301 | #if (LDBL_MANT_DIG == 53) |
304 | __weak_reference(lgamma_r, lgammal_r); | 302 | __weak_reference(lgamma_r, lgammal_r); |
305 | #endif | 303 | #endif |
306 | |||
diff --git a/libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c b/libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c index 9d23053b..9084e18d 100644 --- a/libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c +++ b/libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c | |||
@@ -122,29 +122,29 @@ sin_pif(float x) | |||
122 | float | 122 | float |
123 | __ieee754_lgammaf_r(float x, int *signgamp) | 123 | __ieee754_lgammaf_r(float x, int *signgamp) |
124 | { | 124 | { |
125 | float t,y,z,nadj,p,p1,p2,p3,q,r,w; | 125 | float nadj,p,p1,p2,p3,q,r,t,w,y,z; |
126 | int32_t hx; | 126 | int32_t hx; |
127 | int i,ix; | 127 | int i,ix; |
128 | 128 | ||
129 | GET_FLOAT_WORD(hx,x); | 129 | GET_FLOAT_WORD(hx,x); |
130 | 130 | ||
131 | /* purge off +-inf, NaN, +-0, tiny and negative arguments */ | 131 | /* purge +-Inf and NaNs */ |
132 | *signgamp = 1; | 132 | *signgamp = 1; |
133 | ix = hx&0x7fffffff; | 133 | ix = hx&0x7fffffff; |
134 | if(ix>=0x7f800000) return x*x; | 134 | if(ix>=0x7f800000) return x*x; |
135 | if(ix==0) { | 135 | |
136 | if(hx<0) | 136 | /* purge +-0 and tiny arguments */ |
137 | *signgamp = -1; | 137 | *signgamp = 1-2*((uint32_t)hx>>31); |
138 | return one/vzero; | 138 | if(ix<0x32000000) { /* |x|<2**-27, return -log(|x|) */ |
139 | } | 139 | if(ix==0) |
140 | if(ix<0x35000000) { /* |x|<2**-21, return -log(|x|) */ | 140 | return one/vzero; |
141 | if(hx<0) { | 141 | return -__ieee754_logf(fabsf(x)); |
142 | *signgamp = -1; | ||
143 | return -__ieee754_logf(-x); | ||
144 | } else return -__ieee754_logf(x); | ||
145 | } | 142 | } |
143 | |||
144 | /* purge negative integers and start evaluation for other x < 0 */ | ||
146 | if(hx<0) { | 145 | if(hx<0) { |
147 | if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */ | 146 | *signgamp = 1; |
147 | if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */ | ||
148 | return one/vzero; | 148 | return one/vzero; |
149 | t = sin_pif(x); | 149 | t = sin_pif(x); |
150 | if(t==zero) return one/vzero; /* -integer */ | 150 | if(t==zero) return one/vzero; /* -integer */ |
@@ -153,7 +153,7 @@ __ieee754_lgammaf_r(float x, int *signgamp) | |||
153 | x = -x; | 153 | x = -x; |
154 | } | 154 | } |
155 | 155 | ||
156 | /* purge off 1 and 2 */ | 156 | /* purge 1 and 2 */ |
157 | if (ix==0x3f800000||ix==0x40000000) r = 0; | 157 | if (ix==0x3f800000||ix==0x40000000) r = 0; |
158 | /* for x < 2.0 */ | 158 | /* for x < 2.0 */ |
159 | else if(ix<0x40000000) { | 159 | else if(ix<0x40000000) { |
@@ -174,17 +174,18 @@ __ieee754_lgammaf_r(float x, int *signgamp) | |||
174 | p1 = a0+z*(a2+z*a4); | 174 | p1 = a0+z*(a2+z*a4); |
175 | p2 = z*(a1+z*(a3+z*a5)); | 175 | p2 = z*(a1+z*(a3+z*a5)); |
176 | p = y*p1+p2; | 176 | p = y*p1+p2; |
177 | r += (p-y/2); break; | 177 | r += p-y/2; break; |
178 | case 1: | 178 | case 1: |
179 | p = t0+y*t1+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*t7))))); | 179 | p = t0+y*t1+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*t7))))); |
180 | r += (tf + p); break; | 180 | r += tf + p; break; |
181 | case 2: | 181 | case 2: |
182 | p1 = y*(u0+y*(u1+y*u2)); | 182 | p1 = y*(u0+y*(u1+y*u2)); |
183 | p2 = one+y*(v1+y*(v2+y*v3)); | 183 | p2 = one+y*(v1+y*(v2+y*v3)); |
184 | r += (p1/p2-y/2); | 184 | r += p1/p2-y/2; |
185 | } | 185 | } |
186 | } | 186 | } |
187 | else if(ix<0x41000000) { /* x < 8.0 */ | 187 | /* x < 8.0 */ |
188 | else if(ix<0x41000000) { | ||
188 | i = x; | 189 | i = x; |
189 | y = x-i; | 190 | y = x-i; |
190 | p = y*(s0+y*(s1+y*(s2+y*s3))); | 191 | p = y*(s0+y*(s1+y*(s2+y*s3))); |
@@ -199,15 +200,15 @@ __ieee754_lgammaf_r(float x, int *signgamp) | |||
199 | case 3: z *= (y+2); /* FALLTHRU */ | 200 | case 3: z *= (y+2); /* FALLTHRU */ |
200 | r += __ieee754_logf(z); break; | 201 | r += __ieee754_logf(z); break; |
201 | } | 202 | } |
202 | /* 8.0 <= x < 2**24 */ | 203 | /* 8.0 <= x < 2**27 */ |
203 | } else if (ix < 0x4b800000) { | 204 | } else if (ix < 0x4d000000) { |
204 | t = __ieee754_logf(x); | 205 | t = __ieee754_logf(x); |
205 | z = one/x; | 206 | z = one/x; |
206 | y = z*z; | 207 | y = z*z; |
207 | w = w0+z*(w1+y*w2); | 208 | w = w0+z*(w1+y*w2); |
208 | r = (x-half)*(t-one)+w; | 209 | r = (x-half)*(t-one)+w; |
209 | } else | 210 | } else |
210 | /* 2**24 <= x <= inf */ | 211 | /* 2**27 <= x <= inf */ |
211 | r = x*(__ieee754_logf(x)-one); | 212 | r = x*(__ieee754_logf(x)-one); |
212 | if(hx<0) r = nadj - r; | 213 | if(hx<0) r = nadj - r; |
213 | return r; | 214 | return r; |
diff --git a/libm/upstream-freebsd/lib/msun/src/fenv-softfloat.h b/libm/upstream-freebsd/lib/msun/src/fenv-softfloat.h deleted file mode 100644 index 02d2a2c6..00000000 --- a/libm/upstream-freebsd/lib/msun/src/fenv-softfloat.h +++ /dev/null | |||
@@ -1,184 +0,0 @@ | |||
1 | /*- | ||
2 | * Copyright (c) 2004-2011 David Schultz <das@FreeBSD.ORG> | ||
3 | * All rights reserved. | ||
4 | * | ||
5 | * Redistribution and use in source and binary forms, with or without | ||
6 | * modification, are permitted provided that the following conditions | ||
7 | * are met: | ||
8 | * 1. Redistributions of source code must retain the above copyright | ||
9 | * notice, this list of conditions and the following disclaimer. | ||
10 | * 2. Redistributions in binary form must reproduce the above copyright | ||
11 | * notice, this list of conditions and the following disclaimer in the | ||
12 | * documentation and/or other materials provided with the distribution. | ||
13 | * | ||
14 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND | ||
15 | * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | ||
16 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | ||
17 | * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE | ||
18 | * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL | ||
19 | * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS | ||
20 | * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) | ||
21 | * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT | ||
22 | * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY | ||
23 | * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF | ||
24 | * SUCH DAMAGE. | ||
25 | * | ||
26 | * $FreeBSD$ | ||
27 | */ | ||
28 | |||
29 | #ifndef _FENV_H_ | ||
30 | #error "This file is meant to be included only by <fenv.h>." | ||
31 | #endif | ||
32 | |||
33 | /* | ||
34 | * This file implements the functionality of <fenv.h> on platforms that | ||
35 | * lack an FPU and use softfloat in libc for floating point. To use it, | ||
36 | * you must write an <fenv.h> that provides the following: | ||
37 | * | ||
38 | * - a typedef for fenv_t, which may be an integer or struct type | ||
39 | * - a typedef for fexcept_t (XXX This file assumes fexcept_t is a | ||
40 | * simple integer type containing the exception mask.) | ||
41 | * - definitions of FE_* constants for the five exceptions and four | ||
42 | * rounding modes in IEEE 754, as described in fenv(3) | ||
43 | * - a definition, and the corresponding external symbol, for FE_DFL_ENV | ||
44 | * - a macro __set_env(env, flags, mask, rnd), which sets the given fenv_t | ||
45 | * from the exception flags, mask, and rounding mode | ||
46 | * - macros __env_flags(env), __env_mask(env), and __env_round(env), which | ||
47 | * extract fields from an fenv_t | ||
48 | * - a definition of __fenv_static | ||
49 | * | ||
50 | * If the architecture supports an optional FPU, it's recommended that you | ||
51 | * define fenv_t and fexcept_t to match the hardware ABI. Otherwise, it | ||
52 | * doesn't matter how you define them. | ||
53 | */ | ||
54 | |||
55 | extern int __softfloat_float_exception_flags; | ||
56 | extern int __softfloat_float_exception_mask; | ||
57 | extern int __softfloat_float_rounding_mode; | ||
58 | void __softfloat_float_raise(int); | ||
59 | |||
60 | __fenv_static inline int | ||
61 | feclearexcept(int __excepts) | ||
62 | { | ||
63 | |||
64 | __softfloat_float_exception_flags &= ~__excepts; | ||
65 | return (0); | ||
66 | } | ||
67 | |||
68 | __fenv_static inline int | ||
69 | fegetexceptflag(fexcept_t *__flagp, int __excepts) | ||
70 | { | ||
71 | |||
72 | *__flagp = __softfloat_float_exception_flags & __excepts; | ||
73 | return (0); | ||
74 | } | ||
75 | |||
76 | __fenv_static inline int | ||
77 | fesetexceptflag(const fexcept_t *__flagp, int __excepts) | ||
78 | { | ||
79 | |||
80 | __softfloat_float_exception_flags &= ~__excepts; | ||
81 | __softfloat_float_exception_flags |= *__flagp & __excepts; | ||
82 | return (0); | ||
83 | } | ||
84 | |||
85 | __fenv_static inline int | ||
86 | feraiseexcept(int __excepts) | ||
87 | { | ||
88 | |||
89 | __softfloat_float_raise(__excepts); | ||
90 | return (0); | ||
91 | } | ||
92 | |||
93 | __fenv_static inline int | ||
94 | fetestexcept(int __excepts) | ||
95 | { | ||
96 | |||
97 | return (__softfloat_float_exception_flags & __excepts); | ||
98 | } | ||
99 | |||
100 | __fenv_static inline int | ||
101 | fegetround(void) | ||
102 | { | ||
103 | |||
104 | return (__softfloat_float_rounding_mode); | ||
105 | } | ||
106 | |||
107 | __fenv_static inline int | ||
108 | fesetround(int __round) | ||
109 | { | ||
110 | |||
111 | __softfloat_float_rounding_mode = __round; | ||
112 | return (0); | ||
113 | } | ||
114 | |||
115 | __fenv_static inline int | ||
116 | fegetenv(fenv_t *__envp) | ||
117 | { | ||
118 | |||
119 | __set_env(*__envp, __softfloat_float_exception_flags, | ||
120 | __softfloat_float_exception_mask, __softfloat_float_rounding_mode); | ||
121 | return (0); | ||
122 | } | ||
123 | |||
124 | __fenv_static inline int | ||
125 | feholdexcept(fenv_t *__envp) | ||
126 | { | ||
127 | fenv_t __env; | ||
128 | |||
129 | fegetenv(__envp); | ||
130 | __softfloat_float_exception_flags = 0; | ||
131 | __softfloat_float_exception_mask = 0; | ||
132 | return (0); | ||
133 | } | ||
134 | |||
135 | __fenv_static inline int | ||
136 | fesetenv(const fenv_t *__envp) | ||
137 | { | ||
138 | |||
139 | __softfloat_float_exception_flags = __env_flags(*__envp); | ||
140 | __softfloat_float_exception_mask = __env_mask(*__envp); | ||
141 | __softfloat_float_rounding_mode = __env_round(*__envp); | ||
142 | return (0); | ||
143 | } | ||
144 | |||
145 | __fenv_static inline int | ||
146 | feupdateenv(const fenv_t *__envp) | ||
147 | { | ||
148 | int __oflags = __softfloat_float_exception_flags; | ||
149 | |||
150 | fesetenv(__envp); | ||
151 | feraiseexcept(__oflags); | ||
152 | return (0); | ||
153 | } | ||
154 | |||
155 | #if __BSD_VISIBLE | ||
156 | |||
157 | /* We currently provide no external definitions of the functions below. */ | ||
158 | |||
159 | static inline int | ||
160 | feenableexcept(int __mask) | ||
161 | { | ||
162 | int __omask = __softfloat_float_exception_mask; | ||
163 | |||
164 | __softfloat_float_exception_mask |= __mask; | ||
165 | return (__omask); | ||
166 | } | ||
167 | |||
168 | static inline int | ||
169 | fedisableexcept(int __mask) | ||
170 | { | ||
171 | int __omask = __softfloat_float_exception_mask; | ||
172 | |||
173 | __softfloat_float_exception_mask &= ~__mask; | ||
174 | return (__omask); | ||
175 | } | ||
176 | |||
177 | static inline int | ||
178 | fegetexcept(void) | ||
179 | { | ||
180 | |||
181 | return (__softfloat_float_exception_mask); | ||
182 | } | ||
183 | |||
184 | #endif /* __BSD_VISIBLE */ | ||
diff --git a/libm/upstream-freebsd/lib/msun/src/s_remquol.c b/libm/upstream-freebsd/lib/msun/src/s_remquol.c index 8899293b..712651c9 100644 --- a/libm/upstream-freebsd/lib/msun/src/s_remquol.c +++ b/libm/upstream-freebsd/lib/msun/src/s_remquol.c | |||
@@ -5,7 +5,7 @@ | |||
5 | * | 5 | * |
6 | * Developed at SunSoft, a Sun Microsystems, Inc. business. | 6 | * Developed at SunSoft, a Sun Microsystems, Inc. business. |
7 | * Permission to use, copy, modify, and distribute this | 7 | * Permission to use, copy, modify, and distribute this |
8 | * software is freely granted, provided that this notice | 8 | * software is freely granted, provided that this notice |
9 | * is preserved. | 9 | * is preserved. |
10 | * ==================================================== | 10 | * ==================================================== |
11 | */ | 11 | */ |