am 8dd3b65f: Merge "Sync libm with upstream."

* commit '8dd3b65f2dd63df1baa83285f607cf110f21c27c':
  Sync libm with upstream.

libm/upstream-freebsd/lib/msun/ld128/e_lgammal_r.c

@@ -206,7 +206,7 @@ sin_pil(long double x)
 	    n--;
 	}
 	n &= 7;
-	y = y - z + n * 0.25L;
+	y = y - z + n * 0.25;
 
 	switch (n) {
 	    case 0:   y =  __kernel_sinl(pi*y,zero,0); break;
@@ -221,41 +221,33 @@ sin_pil(long double x)
 	return -y;
 }
 
-
 long double
 lgammal_r(long double x, int *signgamp)
 {
 	long double nadj,p,p1,p2,p3,q,r,t,w,y,z;
 	uint64_t llx,lx;
 	int i;
-	uint16_t hx;
+	uint16_t hx,ix;
 
-	EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
+	EXTRACT_LDBL128_WORDS(hx,lx,llx,x);
 
-	if((hx & 0x7fff) == 0x7fff) {	/* erfl(nan)=nan */
+    /* purge +-Inf and NaNs */
-		i = (hx>>15)<<1;
-		return (1-i)+one/x;	/* erfl(+-inf)=+-1 */
-	}
-
-    /* purge off +-inf, NaN, +-0, tiny and negative arguments */
 	*signgamp = 1;
-	if((hx & 0x7fff) == 0x7fff)	/* x is +-Inf or NaN */
+	ix = hx&0x7fff;
-		return x*x;
+	if(ix==0x7fff) return x*x;
-	if((hx==0||hx==0x8000)&&lx==0) {
+
-	    if (hx&0x8000)
+   /* purge +-0 and tiny arguments */
-		*signgamp = -1;
+	*signgamp = 1-2*(hx>>15);
-	    return one/vzero;
+	if(ix<0x3fff-116) {		/* |x|<2**-(p+3), return -log(|x|) */
+	    if((ix|lx|llx)==0)
+		return one/vzero;
+	    return -logl(fabsl(x));
 	}
 
-   /* purge off tiny and negative arguments */
+    /* purge negative integers and start evaluation for other x < 0 */
-	if(fabsl(x)<0x1p-119L) {
-	    if(hx&0x8000) {
-	        *signgamp = -1;
-	        return -logl(-x);
-	    } else return -logl(x);
-	}
 	if(hx&0x8000) {
-	    if(fabsl(x)>=0x1p112)
+	    *signgamp = 1;
+	    if(ix>=0x3fff+112) 		/* |x|>=2**(p-1), must be -integer */
 		return one/vzero;
 	    t = sin_pil(x);
 	    if(t==zero) return one/vzero;
@@ -264,17 +256,19 @@ lgammal_r(long double x, int *signgamp)
 	    x = -x;
 	}
 
-	if(x == 1 || x ==2) r = 0;
+    /* purge 1 and 2 */
-	else if(x<2) {
+	if((ix==0x3fff || ix==0x4000) && (lx|llx)==0) r = 0;
-	    if(x<=0.8999996185302734) {
+    /* for x < 2.0 */
+	else if(ix<0x4000) {
+	    if(x<=8.9999961853027344e-01) {
 		r = -logl(x);
-		if(x>=0.7315998077392578) {y = 1-x; i= 0;}
+		if(x>=7.3159980773925781e-01) {y = 1-x; i= 0;}
-		else if(x>=0.2316399812698364) {y= x-(tc-1); i=1;}
+		else if(x>=2.3163998126983643e-01) {y= x-(tc-1); i=1;}
 	  	else {y = x; i=2;}
 	    } else {
-	  	r = 0;
+		r = 0;
-	        if(x>=1.7316312789916992) {y=2-x;i=0;}
+	        if(x>=1.7316312789916992e+00) {y=2-x;i=0;}
-	        else if(x>=1.2316322326660156) {y=x-tc;i=1;}
+	        else if(x>=1.2316322326660156e+00) {y=x-tc;i=1;}
 		else {y=x-1;i=2;}
 	    }
 	    switch(i) {
@@ -285,23 +279,24 @@ lgammal_r(long double x, int *signgamp)
 		p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*(a13+z*(a15+
 		    z*(a17+z*(a19+z*(a21+z*a23)))))))))));
 		p  = y*p1+p2;
-		r  += (p-y/2); break;
+		r  += p-y/2; break;
 	      case 1:
 		p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+
 		    y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+
 		    y*(t17+y*(t18+y*(t19+y*(t20+y*(t21+y*(t22+y*(t23+
 		    y*(t24+y*(t25+y*(t26+y*(t27+y*(t28+y*(t29+y*(t30+
 		    y*(t31+y*t32))))))))))))))))))))))))))))));
-		r += (tf + p); break;
+		r += tf + p; break;
 	      case 2:
 		p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*(u6+y*(u7+
 		    y*(u8+y*(u9+y*u10))))))))));
 		p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*(v6+y*(v7+
 		    y*(v8+y*(v9+y*(v10+y*v11))))))))));
-		r += (-y/2 + p1/p2);
+		r += p1/p2-y/2;
 	    }
 	}
-	else if(x<8) {
+    /* x < 8.0 */
+	else if(ix<0x4002) {
 	    i = x;
 	    y = x-i;
 	    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)
 	    case 3: z *= (y+2);		/* FALLTHRU */
 		    r += logl(z); break;
 	    }
-	} else if (x < 0x1p119L) {
+    /* 8.0 <= x < 2**(p+3) */
+	} else if (ix<0x3fff+116) {
 	    t = logl(x);
 	    z = one/x;
 	    y = z*z;
@@ -326,6 +322,7 @@ lgammal_r(long double x, int *signgamp)
 		y*(w9+y*(w10+y*(w11+y*(w12+y*(w13+y*(w14+y*(w15+y*(w16+
 		y*(w17+y*w18)))))))))))))))));
 	    r = (x-half)*(t-one)+w;
+    /* 2**(p+3) <= x <= inf */
 	} else 
 	    r =  x*(logl(x)-1);
 	if(hx&0x8000) r = nadj - r;

libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c

@@ -1,4 +1,3 @@
-
 /* @(#)e_lgamma_r.c 1.3 95/01/18 */
 /*
  * ====================================================
@@ -9,7 +8,6 @@
  * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
- *
  */
 
 #include <sys/cdefs.h>
@@ -80,7 +78,6 @@ __FBSDID("$FreeBSD$");
  *		lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero
  *		lgamma(inf) = inf
  *		lgamma(-inf) = inf (bug for bug compatible with C99!?)
- *	
  */
 
 #include <float.h>
@@ -204,28 +201,28 @@ sin_pi(double x)
 double
 __ieee754_lgamma_r(double x, int *signgamp)
 {
-	double t,y,z,nadj,p,p1,p2,p3,q,r,w;
+	double nadj,p,p1,p2,p3,q,r,t,w,y,z;
 	int32_t hx;
 	int i,ix,lx;
 
 	EXTRACT_WORDS(hx,lx,x);
 
-    /* purge off +-inf, NaN, +-0, tiny and negative arguments */
+    /* purge +-Inf and NaNs */
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
 	if(ix>=0x7ff00000) return x*x;
-	if((ix|lx)==0) {
+
-	   if(hx<0)
+    /* purge +-0 and tiny arguments */
-		*signgamp = -1;
+	*signgamp = 1-2*((uint32_t)hx>>31);
-	    return one/vzero;
+	if(ix<0x3c700000) {	/* |x|<2**-56, return -log(|x|) */
-	}
+	    if((ix|lx)==0)
-	if(ix<0x3b900000) {	/* |x|<2**-70, return -log(|x|) */
+	        return one/vzero;
-	    if(hx<0) {
+	    return -__ieee754_log(fabs(x));
-	        *signgamp = -1;
-	        return -__ieee754_log(-x);
-	    } else return -__ieee754_log(x);
 	}
+
+    /* purge negative integers and start evaluation for other x < 0 */
 	if(hx<0) {
+	    *signgamp = 1;
 	    if(ix>=0x43300000) 	/* |x|>=2**52, must be -integer */
 		return one/vzero;
 	    t = sin_pi(x);
@@ -235,7 +232,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
 	    x = -x;
 	}
 
-    /* purge off 1 and 2 */
+    /* purge 1 and 2 */
 	if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0;
     /* for x < 2.0 */
 	else if(ix<0x40000000) {
@@ -256,7 +253,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
 		p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
 		p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
 		p  = y*p1+p2;
-		r  += (p-y/2); break;
+		r  += p-y/2; break;
 	      case 1:
 		z = y*y;
 		w = z*y;
@@ -264,19 +261,20 @@ __ieee754_lgamma_r(double x, int *signgamp)
 		p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
 		p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
 		p  = z*p1-(tt-w*(p2+y*p3));
-		r += (tf + p); break;
+		r += tf + p; break;
 	      case 2:
 		p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
 		p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
-		r += (-0.5*y + p1/p2);
+		r += p1/p2-y/2;
 	    }
 	}
-	else if(ix<0x40200000) { 			/* x < 8.0 */
+    /* x < 8.0 */
-	    i = (int)x;
+	else if(ix<0x40200000) {
-	    y = x-(double)i;
+	    i = x;
+	    y = x-i;
 	    p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
 	    q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
-	    r = half*y+p/q;
+	    r = y/2+p/q;
 	    z = one;	/* lgamma(1+s) = log(s) + lgamma(s) */
 	    switch(i) {
 	    case 7: z *= (y+6);		/* FALLTHRU */
@@ -286,15 +284,15 @@ __ieee754_lgamma_r(double x, int *signgamp)
 	    case 3: z *= (y+2);		/* FALLTHRU */
 		    r += __ieee754_log(z); break;
 	    }
-    /* 8.0 <= x < 2**58 */
+    /* 8.0 <= x < 2**56 */
-	} else if (ix < 0x43900000) {
+	} else if (ix < 0x43700000) {
 	    t = __ieee754_log(x);
 	    z = one/x;
 	    y = z*z;
 	    w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
 	    r = (x-half)*(t-one)+w;
 	} else
-    /* 2**58 <= x <= inf */
+    /* 2**56 <= x <= inf */
 	    r =  x*(__ieee754_log(x)-one);
 	if(hx<0) r = nadj - r;
 	return r;
@@ -303,4 +301,3 @@ __ieee754_lgamma_r(double x, int *signgamp)
 #if (LDBL_MANT_DIG == 53)
 __weak_reference(lgamma_r, lgammal_r);
 #endif
-

libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c

@@ -122,29 +122,29 @@ sin_pif(float x)
 float
 __ieee754_lgammaf_r(float x, int *signgamp)
 {
-	float t,y,z,nadj,p,p1,p2,p3,q,r,w;
+	float nadj,p,p1,p2,p3,q,r,t,w,y,z;
 	int32_t hx;
 	int i,ix;
 
 	GET_FLOAT_WORD(hx,x);
 
-    /* purge off +-inf, NaN, +-0, tiny and negative arguments */
+    /* purge +-Inf and NaNs */
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
 	if(ix>=0x7f800000) return x*x;
-	if(ix==0) {
+
-	    if(hx<0)
+    /* purge +-0 and tiny arguments */
-		*signgamp = -1;
+	*signgamp = 1-2*((uint32_t)hx>>31);
-	    return one/vzero;
+	if(ix<0x32000000) {		/* |x|<2**-27, return -log(|x|) */
-	}
+	    if(ix==0)
-	if(ix<0x35000000) {	/* |x|<2**-21, return -log(|x|) */
+	        return one/vzero;
-	    if(hx<0) {
+	    return -__ieee754_logf(fabsf(x));
-	        *signgamp = -1;
-	        return -__ieee754_logf(-x);
-	    } else return -__ieee754_logf(x);
 	}
+
+    /* purge negative integers and start evaluation for other x < 0 */
 	if(hx<0) {
-	    if(ix>=0x4b000000) 	/* |x|>=2**23, must be -integer */
+	    *signgamp = 1;
+	    if(ix>=0x4b000000) 		/* |x|>=2**23, must be -integer */
 		return one/vzero;
 	    t = sin_pif(x);
 	    if(t==zero) return one/vzero; /* -integer */
@@ -153,7 +153,7 @@ __ieee754_lgammaf_r(float x, int *signgamp)
 	    x = -x;
 	}
 
-    /* purge off 1 and 2 */
+    /* purge 1 and 2 */
 	if (ix==0x3f800000||ix==0x40000000) r = 0;
     /* for x < 2.0 */
 	else if(ix<0x40000000) {
@@ -174,17 +174,18 @@ __ieee754_lgammaf_r(float x, int *signgamp)
 		p1 = a0+z*(a2+z*a4);
 		p2 = z*(a1+z*(a3+z*a5));
 		p  = y*p1+p2;
-		r  += (p-y/2); break;
+		r  += p-y/2; break;
 	      case 1:
 		p = t0+y*t1+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*t7)))));
-		r += (tf + p); break;
+		r += tf + p; break;
 	      case 2:
 		p1 = y*(u0+y*(u1+y*u2));
 		p2 = one+y*(v1+y*(v2+y*v3));
-		r += (p1/p2-y/2);
+		r += p1/p2-y/2;
 	    }
 	}
-	else if(ix<0x41000000) { 			/* x < 8.0 */
+    /* x < 8.0 */
+	else if(ix<0x41000000) {
 	    i = x;
 	    y = x-i;
 	    p = y*(s0+y*(s1+y*(s2+y*s3)));
@@ -199,15 +200,15 @@ __ieee754_lgammaf_r(float x, int *signgamp)
 	    case 3: z *= (y+2);		/* FALLTHRU */
 		    r += __ieee754_logf(z); break;
 	    }
-    /* 8.0 <= x < 2**24 */
+    /* 8.0 <= x < 2**27 */
-	} else if (ix < 0x4b800000) {
+	} else if (ix < 0x4d000000) {
 	    t = __ieee754_logf(x);
 	    z = one/x;
 	    y = z*z;
 	    w = w0+z*(w1+y*w2);
 	    r = (x-half)*(t-one)+w;
 	} else
-    /* 2**24 <= x <= inf */
+    /* 2**27 <= x <= inf */
 	    r =  x*(__ieee754_logf(x)-one);
 	if(hx<0) r = nadj - r;
 	return r;

libm/upstream-freebsd/lib/msun/src/fenv-softfloat.h

@@ -1,184 +0,0 @@
-/*-
- * Copyright (c) 2004-2011 David Schultz <das@FreeBSD.ORG>
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- *
- * $FreeBSD$
- */
-
-#ifndef	_FENV_H_
-#error "This file is meant to be included only by <fenv.h>."
-#endif
-
-/*
- * This file implements the functionality of <fenv.h> on platforms that
- * lack an FPU and use softfloat in libc for floating point.  To use it,
- * you must write an <fenv.h> that provides the following:
- *
- *   - a typedef for fenv_t, which may be an integer or struct type
- *   - a typedef for fexcept_t (XXX This file assumes fexcept_t is a
- *     simple integer type containing the exception mask.)
- *   - definitions of FE_* constants for the five exceptions and four
- *     rounding modes in IEEE 754, as described in fenv(3)
- *   - a definition, and the corresponding external symbol, for FE_DFL_ENV
- *   - a macro __set_env(env, flags, mask, rnd), which sets the given fenv_t
- *     from the exception flags, mask, and rounding mode
- *   - macros __env_flags(env), __env_mask(env), and __env_round(env), which
- *     extract fields from an fenv_t
- *   - a definition of __fenv_static
- *
- * If the architecture supports an optional FPU, it's recommended that you
- * define fenv_t and fexcept_t to match the hardware ABI.  Otherwise, it
- * doesn't matter how you define them.
- */
-
-extern int __softfloat_float_exception_flags;
-extern int __softfloat_float_exception_mask;
-extern int __softfloat_float_rounding_mode;
-void __softfloat_float_raise(int);
-
-__fenv_static inline int
-feclearexcept(int __excepts)
-{
-
-	__softfloat_float_exception_flags &= ~__excepts;
-	return (0);
-}
-
-__fenv_static inline int
-fegetexceptflag(fexcept_t *__flagp, int __excepts)
-{
-
-	*__flagp = __softfloat_float_exception_flags & __excepts;
-	return (0);
-}
-
-__fenv_static inline int
-fesetexceptflag(const fexcept_t *__flagp, int __excepts)
-{
-
-	__softfloat_float_exception_flags &= ~__excepts;
-	__softfloat_float_exception_flags |= *__flagp & __excepts;
-	return (0);
-}
-
-__fenv_static inline int
-feraiseexcept(int __excepts)
-{
-
-	__softfloat_float_raise(__excepts);
-	return (0);
-}
-
-__fenv_static inline int
-fetestexcept(int __excepts)
-{
-
-	return (__softfloat_float_exception_flags & __excepts);
-}
-
-__fenv_static inline int
-fegetround(void)
-{
-
-	return (__softfloat_float_rounding_mode);
-}
-
-__fenv_static inline int
-fesetround(int __round)
-{
-
-	__softfloat_float_rounding_mode = __round;
-	return (0);
-}
-
-__fenv_static inline int
-fegetenv(fenv_t *__envp)
-{
-
-	__set_env(*__envp, __softfloat_float_exception_flags,
-	    __softfloat_float_exception_mask, __softfloat_float_rounding_mode);
-	return (0);
-}
-
-__fenv_static inline int
-feholdexcept(fenv_t *__envp)
-{
-	fenv_t __env;
-
-	fegetenv(__envp);
-	__softfloat_float_exception_flags = 0;
-	__softfloat_float_exception_mask = 0;
-	return (0);
-}
-
-__fenv_static inline int
-fesetenv(const fenv_t *__envp)
-{
-
-	__softfloat_float_exception_flags = __env_flags(*__envp);
-	__softfloat_float_exception_mask = __env_mask(*__envp);
-	__softfloat_float_rounding_mode = __env_round(*__envp);
-	return (0);
-}
-
-__fenv_static inline int
-feupdateenv(const fenv_t *__envp)
-{
-	int __oflags = __softfloat_float_exception_flags;
-
-	fesetenv(__envp);
-	feraiseexcept(__oflags);
-	return (0);
-}
-
-#if __BSD_VISIBLE
-
-/* We currently provide no external definitions of the functions below. */
-
-static inline int
-feenableexcept(int __mask)
-{
-	int __omask = __softfloat_float_exception_mask;
-
-	__softfloat_float_exception_mask |= __mask;
-	return (__omask);
-}
-
-static inline int
-fedisableexcept(int __mask)
-{
-	int __omask = __softfloat_float_exception_mask;
-
-	__softfloat_float_exception_mask &= ~__mask;
-	return (__omask);
-}
-
-static inline int
-fegetexcept(void)
-{
-
-	return (__softfloat_float_exception_mask);
-}
-
-#endif /* __BSD_VISIBLE */