aboutsummaryrefslogtreecommitdiffstats
path: root/libm
diff options
context:
space:
mode:
authorElliott Hughes2014-10-09 18:20:37 -0500
committerElliott Hughes2014-10-09 18:20:37 -0500
commit488268b134723c7a6598338bb253be5f64d53be4 (patch)
treee2e3adbcb193d4f052907f7be991621469b5dd7c /libm
parente9c216fca56e84b0d0a96f7d5e3c99d3276ef071 (diff)
downloadplatform-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.c73
-rw-r--r--libm/upstream-freebsd/lib/msun/src/e_acoshl.c4
-rw-r--r--libm/upstream-freebsd/lib/msun/src/e_atanhl.c4
-rw-r--r--libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c77
-rw-r--r--libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c43
-rw-r--r--libm/upstream-freebsd/lib/msun/src/fenv-softfloat.h184
-rw-r--r--libm/upstream-freebsd/lib/msun/src/s_remquol.c2
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
225long double 224long double
226lgammal_r(long double x, int *signgamp) 225lgammal_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)
204double 201double
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)
122float 122float
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
55extern int __softfloat_float_exception_flags;
56extern int __softfloat_float_exception_mask;
57extern int __softfloat_float_rounding_mode;
58void __softfloat_float_raise(int);
59
60__fenv_static inline int
61feclearexcept(int __excepts)
62{
63
64 __softfloat_float_exception_flags &= ~__excepts;
65 return (0);
66}
67
68__fenv_static inline int
69fegetexceptflag(fexcept_t *__flagp, int __excepts)
70{
71
72 *__flagp = __softfloat_float_exception_flags & __excepts;
73 return (0);
74}
75
76__fenv_static inline int
77fesetexceptflag(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
86feraiseexcept(int __excepts)
87{
88
89 __softfloat_float_raise(__excepts);
90 return (0);
91}
92
93__fenv_static inline int
94fetestexcept(int __excepts)
95{
96
97 return (__softfloat_float_exception_flags & __excepts);
98}
99
100__fenv_static inline int
101fegetround(void)
102{
103
104 return (__softfloat_float_rounding_mode);
105}
106
107__fenv_static inline int
108fesetround(int __round)
109{
110
111 __softfloat_float_rounding_mode = __round;
112 return (0);
113}
114
115__fenv_static inline int
116fegetenv(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
125feholdexcept(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
136fesetenv(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
146feupdateenv(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
159static inline int
160feenableexcept(int __mask)
161{
162 int __omask = __softfloat_float_exception_mask;
163
164 __softfloat_float_exception_mask |= __mask;
165 return (__omask);
166}
167
168static inline int
169fedisableexcept(int __mask)
170{
171 int __omask = __softfloat_float_exception_mask;
172
173 __softfloat_float_exception_mask &= ~__mask;
174 return (__omask);
175}
176
177static inline int
178fegetexcept(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 */