Merge commit 'refs/top-bases/t/bigmem' into t/bigmem

37 files changed, 480 insertions, 364 deletions

diff --git a/sysdeps/ieee754/ldbl-96/e_acoshl.c b/sysdeps/ieee754/ldbl-96/e_acoshl.c
index bbaef680bd..cf9a6db0ef 100644
--- a/sysdeps/ieee754/ldbl-96/e_acoshl.c
+++ b/sysdeps/ieee754/ldbl-96/e_acoshl.c
@@ -48,7 +48,7 @@ __ieee754_acoshl(long double x)
 		return x+x;
 	    } else
 		return __ieee754_logl(x)+ln2;	/* acoshl(huge)=logl(2x) */
-	} else if(((se-0x3fff)|i0|i1)==0) {
+	} else if(((se-0x3fff)|(i0^0x80000000)|i1)==0) {
 	    return 0.0;			/* acosh(1) = 0 */
 	} else if (se > 0x4000) {	/* 2**28 > x > 2 */
 	    t=x*x;
diff --git a/sysdeps/ieee754/ldbl-96/e_asinl.c b/sysdeps/ieee754/ldbl-96/e_asinl.c
index c1ffa3e0d9..2973bf071b 100644
--- a/sysdeps/ieee754/ldbl-96/e_asinl.c
+++ b/sysdeps/ieee754/ldbl-96/e_asinl.c
@@ -58,6 +58,7 @@
  */
 
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -111,6 +112,11 @@ __ieee754_asinl (long double x)
     {				/* |x|<0.5 */
       if (ix < 0x3fde8000)
 	{			/* if |x| < 2**-33 */
+	  if (fabsl (x) < LDBL_MIN)
+	    {
+	      long double force_underflow = x * x;
+	      math_force_eval (force_underflow);
+	    }
 	  if (huge + x > one)
 	    return x;		/* return x with inexact if x!=0 */
 	}
diff --git a/sysdeps/ieee754/ldbl-96/e_atanhl.c b/sysdeps/ieee754/ldbl-96/e_atanhl.c
index 4c29243695..9a957c9065 100644
--- a/sysdeps/ieee754/ldbl-96/e_atanhl.c
+++ b/sysdeps/ieee754/ldbl-96/e_atanhl.c
@@ -32,6 +32,7 @@
  *
  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -52,9 +53,14 @@ __ieee754_atanhl(long double x)
 	    return (x-x)/(x-x);
 	if(ix==0x3fff)
 	    return x/zero;
-	if(ix<0x3fe3) {
+	if(ix<0x3fdf) {
 	    math_force_eval(huge+x);
-	    return x;	/* x<2**-28 */
+	    if (fabsl (x) < LDBL_MIN)
+	      {
+		long double force_underflow = x * x;
+		math_force_eval (force_underflow);
+	      }
+	    return x;	/* x<2**-32 */
 	}
 	SET_LDOUBLE_EXP(x,ix);
 	if(ix<0x3ffe) {		/* x < 0.5 */
diff --git a/sysdeps/ieee754/ldbl-96/e_coshl.c b/sysdeps/ieee754/ldbl-96/e_coshl.c
index 995c9922cb..dd22cae363 100644
--- a/sysdeps/ieee754/ldbl-96/e_coshl.c
+++ b/sysdeps/ieee754/ldbl-96/e_coshl.c
@@ -54,9 +54,9 @@ __ieee754_coshl (long double x)
 	if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) {
 	    /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
 		if(ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) {
+		    if (ex<0x3fbc) return one;	/* cosh(tiny) = 1 */
 		    t = __expm1l(fabsl(x));
 		    w = one+t;
-		    if (ex<0x3fbc) return w;	/* cosh(tiny) = 1 */
 		    return one+(t*t)/(w+w);
 		}
 
diff --git a/sysdeps/ieee754/ldbl-96/e_gammal_r.c b/sysdeps/ieee754/ldbl-96/e_gammal_r.c
index 477c3a61d3..9da5db33f0 100644
--- a/sysdeps/ieee754/ldbl-96/e_gammal_r.c
+++ b/sysdeps/ieee754/ldbl-96/e_gammal_r.c
@@ -1,5 +1,5 @@
 /* Implementation of gamma function according to ISO C.
-   Copyright (C) 1997-2014 Free Software Foundation, Inc.
+   Copyright (C) 1997-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -102,7 +102,7 @@ gammal_positive (long double x, int *exp2_adj)
 			 * __ieee754_expl (-x_adj)
 			 * __ieee754_sqrtl (2 * M_PIl / x_adj)
 			 / prod);
-      exp_adj += x_eps * __ieee754_logl (x);
+      exp_adj += x_eps * __ieee754_logl (x_adj);
       long double bsum = gamma_coeff[NCOEFF - 1];
       long double x_adj2 = x_adj * x_adj;
       for (size_t i = 1; i <= NCOEFF - 1; i++)
@@ -116,22 +116,23 @@ long double
 __ieee754_gammal_r (long double x, int *signgamp)
 {
   u_int32_t es, hx, lx;
+  long double ret;
 
   GET_LDOUBLE_WORDS (es, hx, lx, x);
 
-  if (__builtin_expect (((es & 0x7fff) | hx | lx) == 0, 0))
+  if (__glibc_unlikely (((es & 0x7fff) | hx | lx) == 0))
     {
       /* Return value for x == 0 is Inf with divide by zero exception.  */
       *signgamp = 0;
       return 1.0 / x;
     }
-  if (__builtin_expect (es == 0xffffffff && ((hx & 0x7fffffff) | lx) == 0, 0))
+  if (__glibc_unlikely (es == 0xffffffff && ((hx & 0x7fffffff) | lx) == 0))
     {
       /* x == -Inf.  According to ISO this is NaN.  */
       *signgamp = 0;
       return x - x;
     }
-  if (__builtin_expect ((es & 0x7fff) == 0x7fff, 0))
+  if (__glibc_unlikely ((es & 0x7fff) == 0x7fff))
     {
       /* Positive infinity (return positive infinity) or NaN (return
 	 NaN).  */
@@ -151,35 +152,58 @@ __ieee754_gammal_r (long double x, int *signgamp)
       *signgamp = 0;
       return LDBL_MAX * LDBL_MAX;
     }
-  else if (x > 0.0L)
+  else
     {
-      *signgamp = 0;
-      int exp2_adj;
-      long double ret = gammal_positive (x, &exp2_adj);
-      return __scalbnl (ret, exp2_adj);
+      SET_RESTORE_ROUNDL (FE_TONEAREST);
+      if (x > 0.0L)
+	{
+	  *signgamp = 0;
+	  int exp2_adj;
+	  ret = gammal_positive (x, &exp2_adj);
+	  ret = __scalbnl (ret, exp2_adj);
+	}
+      else if (x >= -LDBL_EPSILON / 4.0L)
+	{
+	  *signgamp = 0;
+	  ret = 1.0L / x;
+	}
+      else
+	{
+	  long double tx = __truncl (x);
+	  *signgamp = (tx == 2.0L * __truncl (tx / 2.0L)) ? -1 : 1;
+	  if (x <= -1766.0L)
+	    /* Underflow.  */
+	    ret = LDBL_MIN * LDBL_MIN;
+	  else
+	    {
+	      long double frac = tx - x;
+	      if (frac > 0.5L)
+		frac = 1.0L - frac;
+	      long double sinpix = (frac <= 0.25L
+				    ? __sinl (M_PIl * frac)
+				    : __cosl (M_PIl * (0.5L - frac)));
+	      int exp2_adj;
+	      ret = M_PIl / (-x * sinpix
+			     * gammal_positive (-x, &exp2_adj));
+	      ret = __scalbnl (ret, -exp2_adj);
+	    }
+	}
     }
-  else if (x >= -LDBL_EPSILON / 4.0L)
+  if (isinf (ret) && x != 0)
     {
-      *signgamp = 0;
-      return 1.0f / x;
+      if (*signgamp < 0)
+	return -(-__copysignl (LDBL_MAX, ret) * LDBL_MAX);
+      else
+	return __copysignl (LDBL_MAX, ret) * LDBL_MAX;
     }
-  else
+  else if (ret == 0)
     {
-      long double tx = __truncl (x);
-      *signgamp = (tx == 2.0L * __truncl (tx / 2.0L)) ? -1 : 1;
-      if (x <= -1766.0L)
-	/* Underflow.  */
-	return LDBL_MIN * LDBL_MIN;
-      long double frac = tx - x;
-      if (frac > 0.5L)
-	frac = 1.0L - frac;
-      long double sinpix = (frac <= 0.25L
-			    ? __sinl (M_PIl * frac)
-			    : __cosl (M_PIl * (0.5L - frac)));
-      int exp2_adj;
-      long double ret = M_PIl / (-x * sinpix
-				 * gammal_positive (-x, &exp2_adj));
-      return __scalbnl (ret, -exp2_adj);
+      if (*signgamp < 0)
+	return -(-__copysignl (LDBL_MIN, ret) * LDBL_MIN);
+      else
+	return __copysignl (LDBL_MIN, ret) * LDBL_MIN;
     }
+  else
+    return ret;
 }
 strong_alias (__ieee754_gammal_r, __gammal_r_finite)
diff --git a/sysdeps/ieee754/ldbl-96/e_j0l.c b/sysdeps/ieee754/ldbl-96/e_j0l.c
index 3939259d81..a536054cde 100644
--- a/sysdeps/ieee754/ldbl-96/e_j0l.c
+++ b/sysdeps/ieee754/ldbl-96/e_j0l.c
@@ -112,7 +112,7 @@ __ieee754_j0l (long double x)
 
   GET_LDOUBLE_EXP (se, x);
   ix = se & 0x7fff;
-  if (__builtin_expect (ix >= 0x7fff, 0))
+  if (__glibc_unlikely (ix >= 0x7fff))
     return one / (x * x);
   x = fabsl (x);
   if (ix >= 0x4000)		/* |x| >= 2.0 */
@@ -132,7 +132,7 @@ __ieee754_j0l (long double x)
        * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
        * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
        */
-      if (__builtin_expect (ix > 0x4080, 0))	/* 2^129 */
+      if (__glibc_unlikely (ix > 0x4080))      	/* 2^129 */
 	z = (invsqrtpi * cc) / __ieee754_sqrtl (x);
       else
 	{
@@ -142,7 +142,7 @@ __ieee754_j0l (long double x)
 	}
       return z;
     }
-  if (__builtin_expect (ix < 0x3fef, 0)) /* |x| < 2**-16 */
+  if (__glibc_unlikely (ix < 0x3fef))       /* |x| < 2**-16 */
     {
       /* raise inexact if x != 0 */
       math_force_eval (huge + x);
@@ -199,11 +199,11 @@ __ieee754_y0l (long double x)
   GET_LDOUBLE_WORDS (se, i0, i1, x);
   ix = se & 0x7fff;
   /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
-  if (__builtin_expect (se & 0x8000, 0))
+  if (__glibc_unlikely (se & 0x8000))
     return zero / (zero * x);
-  if (__builtin_expect (ix >= 0x7fff, 0))
+  if (__glibc_unlikely (ix >= 0x7fff))
     return one / (x + x * x);
-  if (__builtin_expect ((i0 | i1) == 0, 0))
+  if (__glibc_unlikely ((i0 | i1) == 0))
     return -HUGE_VALL + x;  /* -inf and overflow exception.  */
   if (ix >= 0x4000)
     {				/* |x| >= 2.0 */
@@ -234,7 +234,7 @@ __ieee754_y0l (long double x)
 	  else
 	    ss = z / cc;
 	}
-      if (__builtin_expect (ix > 0x4080, 0))	/* 1e39 */
+      if (__glibc_unlikely (ix > 0x4080))      	/* 1e39 */
 	z = (invsqrtpi * ss) / __ieee754_sqrtl (x);
       else
 	{
@@ -244,7 +244,7 @@ __ieee754_y0l (long double x)
 	}
       return z;
     }
-  if (__builtin_expect (ix <= 0x3fde, 0)) /* x < 2^-33 */
+  if (__glibc_unlikely (ix <= 0x3fde))       /* x < 2^-33 */
     {
       z = -7.380429510868722527629822444004602747322E-2L
 	+ tpi * __ieee754_logl (x);
@@ -356,6 +356,7 @@ pzero (long double x)
 
   GET_LDOUBLE_WORDS (se, i0, i1, x);
   ix = se & 0x7fff;
+  /* ix >= 0x4000 for all calls to this function.  */
   if (ix >= 0x4002)
     {
       p = pR8;
@@ -374,7 +375,7 @@ pzero (long double x)
 	  p = pR3;
 	  q = pS3;
 	}
-      else if (ix >= 0x4000)	/* x better be >= 2 */
+      else	/* x >= 2 */
 	{
 	  p = pR2;
 	  q = pS2;
@@ -493,6 +494,7 @@ qzero (long double x)
 
   GET_LDOUBLE_WORDS (se, i0, i1, x);
   ix = se & 0x7fff;
+  /* ix >= 0x4000 for all calls to this function.  */
   if (ix >= 0x4002)		/* x >= 8 */
     {
       p = qR8;
@@ -511,7 +513,7 @@ qzero (long double x)
 	  p = qR3;
 	  q = qS3;
 	}
-      else if (ix >= 0x4000)	/* x better be >= 2 */
+      else	/* x >= 2 */
 	{
 	  p = qR2;
 	  q = qS2;
diff --git a/sysdeps/ieee754/ldbl-96/e_j1l.c b/sysdeps/ieee754/ldbl-96/e_j1l.c
index 4c13018aea..46e28dfe11 100644
--- a/sysdeps/ieee754/ldbl-96/e_j1l.c
+++ b/sysdeps/ieee754/ldbl-96/e_j1l.c
@@ -71,6 +71,8 @@
  *	   by method mentioned above.
  */
 
+#include <errno.h>
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -114,7 +116,7 @@ __ieee754_j1l (long double x)
 
   GET_LDOUBLE_EXP (se, x);
   ix = se & 0x7fff;
-  if (__builtin_expect (ix >= 0x7fff, 0))
+  if (__glibc_unlikely (ix >= 0x7fff))
     return one / x;
   y = fabsl (x);
   if (ix >= 0x4000)
@@ -134,7 +136,7 @@ __ieee754_j1l (long double x)
        * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
        * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
        */
-      if (__builtin_expect (ix > 0x4080, 0))
+      if (__glibc_unlikely (ix > 0x4080))
 	z = (invsqrtpi * cc) / __ieee754_sqrtl (y);
       else
 	{
@@ -147,10 +149,18 @@ __ieee754_j1l (long double x)
       else
 	return z;
     }
-  if (__builtin_expect (ix < 0x3fde, 0)) /* |x| < 2^-33 */
+  if (__glibc_unlikely (ix < 0x3fde))       /* |x| < 2^-33 */
     {
-      if (huge + x > one)
-	return 0.5 * x;		/* inexact if x!=0 necessary */
+      if (huge + x > one)		/* inexact if x!=0 necessary */
+	{
+	  long double ret = 0.5 * x;
+	  if (fabsl (ret) < LDBL_MIN)
+	    {
+	      long double force_underflow = ret * ret;
+	      math_force_eval (force_underflow);
+	    }
+	  return ret;
+	}
     }
   z = x * x;
   r = z * (R[0] + z * (R[1]+ z * (R[2] + z * (R[3] + z * R[4]))));
@@ -192,11 +202,11 @@ __ieee754_y1l (long double x)
   GET_LDOUBLE_WORDS (se, i0, i1, x);
   ix = se & 0x7fff;
   /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
-  if (__builtin_expect (se & 0x8000, 0))
+  if (__glibc_unlikely (se & 0x8000))
     return zero / (zero * x);
-  if (__builtin_expect (ix >= 0x7fff, 0))
+  if (__glibc_unlikely (ix >= 0x7fff))
     return one / (x + x * x);
-  if (__builtin_expect ((i0 | i1) == 0, 0))
+  if (__glibc_unlikely ((i0 | i1) == 0))
     return -HUGE_VALL + x;  /* -inf and overflow exception.  */
   if (ix >= 0x4000)
     {				/* |x| >= 2.0 */
@@ -222,7 +232,7 @@ __ieee754_y1l (long double x)
        *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
        * to compute the worse one.
        */
-      if (__builtin_expect (ix > 0x4080, 0))
+      if (__glibc_unlikely (ix > 0x4080))
 	z = (invsqrtpi * ss) / __ieee754_sqrtl (x);
       else
 	{
@@ -232,9 +242,12 @@ __ieee754_y1l (long double x)
 	}
       return z;
     }
-  if (__builtin_expect (ix <= 0x3fbe, 0))
+  if (__glibc_unlikely (ix <= 0x3fbe))
     {				/* x < 2**-65 */
-      return (-tpi / x);
+      z = -tpi / x;
+      if (isinf (z))
+	__set_errno (ERANGE);
+      return z;
     }
   z = x * x;
  u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * (U0[4] + z * U0[5]))));
@@ -355,6 +368,7 @@ pone (long double x)
 
   GET_LDOUBLE_WORDS (se, i0, i1, x);
   ix = se & 0x7fff;
+  /* ix >= 0x4000 for all calls to this function.  */
   if (ix >= 0x4002) /* x >= 8 */
     {
       p = pr8;
@@ -373,7 +387,7 @@ pone (long double x)
 	  p = pr3;
 	  q = ps3;
 	}
-      else if (ix >= 0x4000)	/* x better be >= 2 */
+      else	/* x >= 2 */
 	{
 	  p = pr2;
 	  q = ps2;
@@ -501,6 +515,7 @@ qone (long double x)
 
   GET_LDOUBLE_WORDS (se, i0, i1, x);
   ix = se & 0x7fff;
+  /* ix >= 0x4000 for all calls to this function.  */
   if (ix >= 0x4002)		/* x >= 8 */
     {
       p = qr8;
@@ -519,7 +534,7 @@ qone (long double x)
 	  p = qr3;
 	  q = qs3;
 	}
-      else if (ix >= 0x4000)	/* x better be >= 2 */
+      else	/* x >= 2 */
 	{
 	  p = qr2;
 	  q = qs2;
diff --git a/sysdeps/ieee754/ldbl-96/e_jnl.c b/sysdeps/ieee754/ldbl-96/e_jnl.c
index fa8e27efec..2f3a452f55 100644
--- a/sysdeps/ieee754/ldbl-96/e_jnl.c
+++ b/sysdeps/ieee754/ldbl-96/e_jnl.c
@@ -57,6 +57,7 @@
  */
 
 #include <errno.h>
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -70,7 +71,7 @@ __ieee754_jnl (int n, long double x)
 {
   u_int32_t se, i0, i1;
   int32_t i, ix, sgn;
-  long double a, b, temp, di;
+  long double a, b, temp, di, ret;
   long double z, w;
 
   /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
@@ -81,7 +82,7 @@ __ieee754_jnl (int n, long double x)
   ix = se & 0x7fff;
 
   /* if J(n,NaN) is NaN */
-  if (__builtin_expect ((ix == 0x7fff) && ((i0 & 0x7fffffff) != 0), 0))
+  if (__glibc_unlikely ((ix == 0x7fff) && ((i0 & 0x7fffffff) != 0)))
     return x + x;
   if (n < 0)
     {
@@ -95,195 +96,206 @@ __ieee754_jnl (int n, long double x)
     return (__ieee754_j1l (x));
   sgn = (n & 1) & (se >> 15);	/* even n -- 0, odd n -- sign(x) */
   x = fabsl (x);
-  if (__builtin_expect ((ix | i0 | i1) == 0 || ix >= 0x7fff, 0))
-    /* if x is 0 or inf */
-    b = zero;
-  else if ((long double) n <= x)
-    {
-      /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
-      if (ix >= 0x412D)
-	{			/* x > 2**302 */
+  {
+    SET_RESTORE_ROUNDL (FE_TONEAREST);
+    if (__glibc_unlikely ((ix | i0 | i1) == 0 || ix >= 0x7fff))
+      /* if x is 0 or inf */
+      return sgn == 1 ? -zero : zero;
+    else if ((long double) n <= x)
+      {
+	/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
+	if (ix >= 0x412D)
+	  {			/* x > 2**302 */
 
-	  /* ??? This might be a futile gesture.
-	     If x exceeds X_TLOSS anyway, the wrapper function
-	     will set the result to zero. */
+	    /* ??? This might be a futile gesture.
+	       If x exceeds X_TLOSS anyway, the wrapper function
+	       will set the result to zero. */
 
-	  /* (x >> n**2)
-	   *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-	   *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-	   *      Let s=sin(x), c=cos(x),
-	   *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
-	   *
-	   *             n    sin(xn)*sqt2    cos(xn)*sqt2
-	   *          ----------------------------------
-	   *             0     s-c             c+s
-	   *             1    -s-c            -c+s
-	   *             2    -s+c            -c-s
-	   *             3     s+c             c-s
-	   */
-	  long double s;
-	  long double c;
-	  __sincosl (x, &s, &c);
-	  switch (n & 3)
-	    {
-	    case 0:
-	      temp = c + s;
-	      break;
-	    case 1:
-	      temp = -c + s;
-	      break;
-	    case 2:
-	      temp = -c - s;
-	      break;
-	    case 3:
-	      temp = c - s;
-	      break;
-	    }
-	  b = invsqrtpi * temp / __ieee754_sqrtl (x);
-	}
-      else
-	{
-	  a = __ieee754_j0l (x);
-	  b = __ieee754_j1l (x);
-	  for (i = 1; i < n; i++)
-	    {
-	      temp = b;
-	      b = b * ((long double) (i + i) / x) - a;	/* avoid underflow */
-	      a = temp;
-	    }
-	}
-    }
-  else
-    {
-      if (ix < 0x3fde)
-	{			/* x < 2**-33 */
-	  /* x is tiny, return the first Taylor expansion of J(n,x)
-	   * J(n,x) = 1/n!*(x/2)^n  - ...
-	   */
-	  if (n >= 400)		/* underflow, result < 10^-4952 */
-	    b = zero;
-	  else
-	    {
-	      temp = x * 0.5;
-	      b = temp;
-	      for (a = one, i = 2; i <= n; i++)
-		{
-		  a *= (long double) i;	/* a = n! */
-		  b *= temp;	/* b = (x/2)^n */
-		}
-	      b = b / a;
-	    }
-	}
-      else
-	{
-	  /* use backward recurrence */
-	  /*                      x      x^2      x^2
-	   *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
-	   *                      2n  - 2(n+1) - 2(n+2)
-	   *
-	   *                      1      1        1
-	   *  (for large x)   =  ----  ------   ------   .....
-	   *                      2n   2(n+1)   2(n+2)
-	   *                      -- - ------ - ------ -
-	   *                       x     x         x
-	   *
-	   * Let w = 2n/x and h=2/x, then the above quotient
-	   * is equal to the continued fraction:
-	   *                  1
-	   *      = -----------------------
-	   *                     1
-	   *         w - -----------------
-	   *                        1
-	   *              w+h - ---------
-	   *                     w+2h - ...
-	   *
-	   * To determine how many terms needed, let
-	   * Q(0) = w, Q(1) = w(w+h) - 1,
-	   * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
-	   * When Q(k) > 1e4      good for single
-	   * When Q(k) > 1e9      good for double
-	   * When Q(k) > 1e17     good for quadruple
-	   */
-	  /* determine k */
-	  long double t, v;
-	  long double q0, q1, h, tmp;
-	  int32_t k, m;
-	  w = (n + n) / (long double) x;
-	  h = 2.0L / (long double) x;
-	  q0 = w;
-	  z = w + h;
-	  q1 = w * z - 1.0L;
-	  k = 1;
-	  while (q1 < 1.0e11L)
-	    {
-	      k += 1;
-	      z += h;
-	      tmp = z * q1 - q0;
-	      q0 = q1;
-	      q1 = tmp;
-	    }
-	  m = n + n;
-	  for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
-	    t = one / (i / x - t);
-	  a = t;
-	  b = one;
-	  /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
-	   *  Hence, if n*(log(2n/x)) > ...
-	   *  single 8.8722839355e+01
-	   *  double 7.09782712893383973096e+02
-	   *  long double 1.1356523406294143949491931077970765006170e+04
-	   *  then recurrent value may overflow and the result is
-	   *  likely underflow to zero
-	   */
-	  tmp = n;
-	  v = two / x;
-	  tmp = tmp * __ieee754_logl (fabsl (v * tmp));
+	    /* (x >> n**2)
+	     *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
+	     *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
+	     *      Let s=sin(x), c=cos(x),
+	     *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
+	     *
+	     *             n    sin(xn)*sqt2    cos(xn)*sqt2
+	     *          ----------------------------------
+	     *             0     s-c             c+s
+	     *             1    -s-c            -c+s
+	     *             2    -s+c            -c-s
+	     *             3     s+c             c-s
+	     */
+	    long double s;
+	    long double c;
+	    __sincosl (x, &s, &c);
+	    switch (n & 3)
+	      {
+	      case 0:
+		temp = c + s;
+		break;
+	      case 1:
+		temp = -c + s;
+		break;
+	      case 2:
+		temp = -c - s;
+		break;
+	      case 3:
+		temp = c - s;
+		break;
+	      }
+	    b = invsqrtpi * temp / __ieee754_sqrtl (x);
+	  }
+	else
+	  {
+	    a = __ieee754_j0l (x);
+	    b = __ieee754_j1l (x);
+	    for (i = 1; i < n; i++)
+	      {
+		temp = b;
+		b = b * ((long double) (i + i) / x) - a;	/* avoid underflow */
+		a = temp;
+	      }
+	  }
+      }
+    else
+      {
+	if (ix < 0x3fde)
+	  {			/* x < 2**-33 */
+	    /* x is tiny, return the first Taylor expansion of J(n,x)
+	     * J(n,x) = 1/n!*(x/2)^n  - ...
+	     */
+	    if (n >= 400)		/* underflow, result < 10^-4952 */
+	      b = zero;
+	    else
+	      {
+		temp = x * 0.5;
+		b = temp;
+		for (a = one, i = 2; i <= n; i++)
+		  {
+		    a *= (long double) i;	/* a = n! */
+		    b *= temp;	/* b = (x/2)^n */
+		  }
+		b = b / a;
+	      }
+	  }
+	else
+	  {
+	    /* use backward recurrence */
+	    /*                      x      x^2      x^2
+	     *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
+	     *                      2n  - 2(n+1) - 2(n+2)
+	     *
+	     *                      1      1        1
+	     *  (for large x)   =  ----  ------   ------   .....
+	     *                      2n   2(n+1)   2(n+2)
+	     *                      -- - ------ - ------ -
+	     *                       x     x         x
+	     *
+	     * Let w = 2n/x and h=2/x, then the above quotient
+	     * is equal to the continued fraction:
+	     *                  1
+	     *      = -----------------------
+	     *                     1
+	     *         w - -----------------
+	     *                        1
+	     *              w+h - ---------
+	     *                     w+2h - ...
+	     *
+	     * To determine how many terms needed, let
+	     * Q(0) = w, Q(1) = w(w+h) - 1,
+	     * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
+	     * When Q(k) > 1e4      good for single
+	     * When Q(k) > 1e9      good for double
+	     * When Q(k) > 1e17     good for quadruple
+	     */
+	    /* determine k */
+	    long double t, v;
+	    long double q0, q1, h, tmp;
+	    int32_t k, m;
+	    w = (n + n) / (long double) x;
+	    h = 2.0L / (long double) x;
+	    q0 = w;
+	    z = w + h;
+	    q1 = w * z - 1.0L;
+	    k = 1;
+	    while (q1 < 1.0e11L)
+	      {
+		k += 1;
+		z += h;
+		tmp = z * q1 - q0;
+		q0 = q1;
+		q1 = tmp;
+	      }
+	    m = n + n;
+	    for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
+	      t = one / (i / x - t);
+	    a = t;
+	    b = one;
+	    /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
+	     *  Hence, if n*(log(2n/x)) > ...
+	     *  single 8.8722839355e+01
+	     *  double 7.09782712893383973096e+02
+	     *  long double 1.1356523406294143949491931077970765006170e+04
+	     *  then recurrent value may overflow and the result is
+	     *  likely underflow to zero
+	     */
+	    tmp = n;
+	    v = two / x;
+	    tmp = tmp * __ieee754_logl (fabsl (v * tmp));
 
-	  if (tmp < 1.1356523406294143949491931077970765006170e+04L)
-	    {
-	      for (i = n - 1, di = (long double) (i + i); i > 0; i--)
-		{
-		  temp = b;
-		  b *= di;
-		  b = b / x - a;
-		  a = temp;
-		  di -= two;
-		}
-	    }
-	  else
-	    {
-	      for (i = n - 1, di = (long double) (i + i); i > 0; i--)
-		{
-		  temp = b;
-		  b *= di;
-		  b = b / x - a;
-		  a = temp;
-		  di -= two;
-		  /* scale b to avoid spurious overflow */
-		  if (b > 1e100L)
-		    {
-		      a /= b;
-		      t /= b;
-		      b = one;
-		    }
-		}
-	    }
-	  /* j0() and j1() suffer enormous loss of precision at and
-	   * near zero; however, we know that their zero points never
-	   * coincide, so just choose the one further away from zero.
-	   */
-	  z = __ieee754_j0l (x);
-	  w = __ieee754_j1l (x);
-	  if (fabsl (z) >= fabsl (w))
-	    b = (t * z / b);
-	  else
-	    b = (t * w / a);
-	}
+	    if (tmp < 1.1356523406294143949491931077970765006170e+04L)
+	      {
+		for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+		  {
+		    temp = b;
+		    b *= di;
+		    b = b / x - a;
+		    a = temp;
+		    di -= two;
+		  }
+	      }
+	    else
+	      {
+		for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+		  {
+		    temp = b;
+		    b *= di;
+		    b = b / x - a;
+		    a = temp;
+		    di -= two;
+		    /* scale b to avoid spurious overflow */
+		    if (b > 1e100L)
+		      {
+			a /= b;
+			t /= b;
+			b = one;
+		      }
+		  }
+	      }
+	    /* j0() and j1() suffer enormous loss of precision at and
+	     * near zero; however, we know that their zero points never
+	     * coincide, so just choose the one further away from zero.
+	     */
+	    z = __ieee754_j0l (x);
+	    w = __ieee754_j1l (x);
+	    if (fabsl (z) >= fabsl (w))
+	      b = (t * z / b);
+	    else
+	      b = (t * w / a);
+	  }
+      }
+    if (sgn == 1)
+      ret = -b;
+    else
+      ret = b;
+  }
+  if (ret == 0)
+    ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN;
+  else if (fabsl (ret) < LDBL_MIN)
+    {
+      long double force_underflow = ret * ret;
+      math_force_eval (force_underflow);
     }
-  if (sgn == 1)
-    return -b;
-  else
-    return b;
+  return ret;
 }
 strong_alias (__ieee754_jnl, __jnl_finite)
 
@@ -293,7 +305,7 @@ __ieee754_ynl (int n, long double x)
   u_int32_t se, i0, i1;
   int32_t i, ix;
   int32_t sign;
-  long double a, b, temp;
+  long double a, b, temp, ret;
 
 
   GET_LDOUBLE_WORDS (se, i0, i1, x);
@@ -314,69 +326,79 @@ __ieee754_ynl (int n, long double x)
     }
   if (n == 0)
     return (__ieee754_y0l (x));
-  if (n == 1)
-    return (sign * __ieee754_y1l (x));
-  if (__builtin_expect (ix == 0x7fff, 0))
-    return zero;
-  if (ix >= 0x412D)
-    {				/* x > 2**302 */
+  {
+    SET_RESTORE_ROUNDL (FE_TONEAREST);
+    if (n == 1)
+      {
+	ret = sign * __ieee754_y1l (x);
+	goto out;
+      }
+    if (__glibc_unlikely (ix == 0x7fff))
+      return zero;
+    if (ix >= 0x412D)
+      {				/* x > 2**302 */
 
-      /* ??? See comment above on the possible futility of this.  */
+	/* ??? See comment above on the possible futility of this.  */
 
-      /* (x >> n**2)
-       *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-       *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-       *      Let s=sin(x), c=cos(x),
-       *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
-       *
-       *             n    sin(xn)*sqt2    cos(xn)*sqt2
-       *          ----------------------------------
-       *             0     s-c             c+s
-       *             1    -s-c            -c+s
-       *             2    -s+c            -c-s
-       *             3     s+c             c-s
-       */
-      long double s;
-      long double c;
-      __sincosl (x, &s, &c);
-      switch (n & 3)
-	{
-	case 0:
-	  temp = s - c;
-	  break;
-	case 1:
-	  temp = -s - c;
-	  break;
-	case 2:
-	  temp = -s + c;
-	  break;
-	case 3:
-	  temp = s + c;
-	  break;
-	}
-      b = invsqrtpi * temp / __ieee754_sqrtl (x);
-    }
-  else
-    {
-      a = __ieee754_y0l (x);
-      b = __ieee754_y1l (x);
-      /* quit if b is -inf */
-      GET_LDOUBLE_WORDS (se, i0, i1, b);
-      /* Use 0xffffffff since GET_LDOUBLE_WORDS sign-extends SE.  */
-      for (i = 1; i < n && se != 0xffffffff; i++)
-	{
-	  temp = b;
-	  b = ((long double) (i + i) / x) * b - a;
-	  GET_LDOUBLE_WORDS (se, i0, i1, b);
-	  a = temp;
-	}
-    }
-  /* If B is +-Inf, set up errno accordingly.  */
-  if (! __finitel (b))
-    __set_errno (ERANGE);
-  if (sign > 0)
-    return b;
-  else
-    return -b;
+	/* (x >> n**2)
+	 *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
+	 *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
+	 *      Let s=sin(x), c=cos(x),
+	 *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
+	 *
+	 *             n    sin(xn)*sqt2    cos(xn)*sqt2
+	 *          ----------------------------------
+	 *             0     s-c             c+s
+	 *             1    -s-c            -c+s
+	 *             2    -s+c            -c-s
+	 *             3     s+c             c-s
+	 */
+	long double s;
+	long double c;
+	__sincosl (x, &s, &c);
+	switch (n & 3)
+	  {
+	  case 0:
+	    temp = s - c;
+	    break;
+	  case 1:
+	    temp = -s - c;
+	    break;
+	  case 2:
+	    temp = -s + c;
+	    break;
+	  case 3:
+	    temp = s + c;
+	    break;
+	  }
+	b = invsqrtpi * temp / __ieee754_sqrtl (x);
+      }
+    else
+      {
+	a = __ieee754_y0l (x);
+	b = __ieee754_y1l (x);
+	/* quit if b is -inf */
+	GET_LDOUBLE_WORDS (se, i0, i1, b);
+	/* Use 0xffffffff since GET_LDOUBLE_WORDS sign-extends SE.  */
+	for (i = 1; i < n && se != 0xffffffff; i++)
+	  {
+	    temp = b;
+	    b = ((long double) (i + i) / x) * b - a;
+	    GET_LDOUBLE_WORDS (se, i0, i1, b);
+	    a = temp;
+	  }
+      }
+    /* If B is +-Inf, set up errno accordingly.  */
+    if (! isfinite (b))
+      __set_errno (ERANGE);
+    if (sign > 0)
+      ret = b;
+    else
+      ret = -b;
+  }
+ out:
+  if (isinf (ret))
+    ret = __copysignl (LDBL_MAX, ret) * LDBL_MAX;
+  return ret;
 }
 strong_alias (__ieee754_ynl, __ynl_finite)
diff --git a/sysdeps/ieee754/ldbl-96/e_lgammal_r.c b/sysdeps/ieee754/ldbl-96/e_lgammal_r.c
index 1a95f846ae..0cc35f9252 100644
--- a/sysdeps/ieee754/ldbl-96/e_lgammal_r.c
+++ b/sysdeps/ieee754/ldbl-96/e_lgammal_r.c
@@ -91,6 +91,7 @@
  *
  */
 
+#include <libc-internal.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -423,8 +424,18 @@ __ieee754_lgammal_r (long double x, int *signgamp)
   else
     /* 2**66 <= x <= inf */
     r = x * (__ieee754_logl (x) - one);
+  /* NADJ is set for negative arguments but not otherwise, resulting
+     in warnings that it may be used uninitialized although in the
+     cases where it is used it has always been set.  */
+  DIAG_PUSH_NEEDS_COMMENT;
+#if __GNUC_PREREQ (4, 7)
+  DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wmaybe-uninitialized");
+#else
+  DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wuninitialized");
+#endif
   if (se & 0x8000)
     r = nadj - r;
+  DIAG_POP_NEEDS_COMMENT;
   return r;
 }
 strong_alias (__ieee754_lgammal_r, __lgammal_r_finite)
diff --git a/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c b/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c
index e18be6ee0c..d05afce623 100644
--- a/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c
+++ b/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c
@@ -1,5 +1,5 @@
 /* Extended-precision floating point argument reduction.
-   Copyright (C) 1999-2014 Free Software Foundation, Inc.
+   Copyright (C) 1999-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Based on quad-precision code by Jakub Jelinek <jj@ultra.linux.cz>
 
diff --git a/sysdeps/ieee754/ldbl-96/gamma_product.c b/sysdeps/ieee754/ldbl-96/gamma_product.c
index 214745624f..e113e23723 100644
--- a/sysdeps/ieee754/ldbl-96/gamma_product.c
+++ b/sysdeps/ieee754/ldbl-96/gamma_product.c
@@ -1,5 +1,5 @@
 /* Compute a product of X, X+1, ..., with an error estimate.
-   Copyright (C) 2013-2014 Free Software Foundation, Inc.
+   Copyright (C) 2013-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
diff --git a/sysdeps/ieee754/ldbl-96/gamma_productl.c b/sysdeps/ieee754/ldbl-96/gamma_productl.c
index 2c87649677..32990b8f1c 100644
--- a/sysdeps/ieee754/ldbl-96/gamma_productl.c
+++ b/sysdeps/ieee754/ldbl-96/gamma_productl.c
@@ -1,5 +1,5 @@
 /* Compute a product of X, X+1, ..., with an error estimate.
-   Copyright (C) 2013-2014 Free Software Foundation, Inc.
+   Copyright (C) 2013-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
diff --git a/sysdeps/ieee754/ldbl-96/k_cosl.c b/sysdeps/ieee754/ldbl-96/k_cosl.c
index daf7c060d2..95d64bd014 100644
--- a/sysdeps/ieee754/ldbl-96/k_cosl.c
+++ b/sysdeps/ieee754/ldbl-96/k_cosl.c
@@ -1,5 +1,5 @@
 /* Extended-precision floating point cosine on <-pi/4,pi/4>.
-   Copyright (C) 1999-2014 Free Software Foundation, Inc.
+   Copyright (C) 1999-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Based on quad-precision cosine by Jakub Jelinek <jj@ultra.linux.cz>
 
diff --git a/sysdeps/ieee754/ldbl-96/k_sinl.c b/sysdeps/ieee754/ldbl-96/k_sinl.c
index f2d1c860e8..b7b5ae359d 100644
--- a/sysdeps/ieee754/ldbl-96/k_sinl.c
+++ b/sysdeps/ieee754/ldbl-96/k_sinl.c
@@ -1,5 +1,5 @@
 /* Quad-precision floating point sine on <-pi/4,pi/4>.
-   Copyright (C) 1999-2014 Free Software Foundation, Inc.
+   Copyright (C) 1999-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Based on quad-precision sine by Jakub Jelinek <jj@ultra.linux.cz>
 
@@ -20,6 +20,7 @@
 /* The polynomials have not been optimized for extended-precision and
    may contain more terms than needed.  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -94,7 +95,14 @@ __kernel_sinl(long double x, long double y, int iy)
       /* Argument is small enough to approximate it by a Chebyshev
 	 polynomial of degree 17.  */
       if (absx < 0x1p-33L)
-	if (!((int)x)) return x;	/* generate inexact */
+	{
+	  if (fabsl (x) < LDBL_MIN)
+	    {
+	      long double force_underflow = x * x;
+	      math_force_eval (force_underflow);
+	    }
+	  if (!((int)x)) return x;	/* generate inexact */
+	}
       z = x * x;
       return x + (x * (z*(SIN1+z*(SIN2+z*(SIN3+z*(SIN4+
 		       z*(SIN5+z*(SIN6+z*(SIN7+z*SIN8)))))))));
diff --git a/sysdeps/ieee754/ldbl-96/ldbl2mpn.c b/sysdeps/ieee754/ldbl-96/ldbl2mpn.c
index be8b93ffa5..26efea18f6 100644
--- a/sysdeps/ieee754/ldbl-96/ldbl2mpn.c
+++ b/sysdeps/ieee754/ldbl-96/ldbl2mpn.c
@@ -1,4 +1,4 @@
-/* Copyright (C) 1995-2014 Free Software Foundation, Inc.
+/* Copyright (C) 1995-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
diff --git a/sysdeps/ieee754/ldbl-96/mpn2ldbl.c b/sysdeps/ieee754/ldbl-96/mpn2ldbl.c
index 8fdb335fe4..b219de0fca 100644
--- a/sysdeps/ieee754/ldbl-96/mpn2ldbl.c
+++ b/sysdeps/ieee754/ldbl-96/mpn2ldbl.c
@@ -1,4 +1,4 @@
-/* Copyright (C) 1995-2014 Free Software Foundation, Inc.
+/* Copyright (C) 1995-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
diff --git a/sysdeps/ieee754/ldbl-96/printf_fphex.c b/sysdeps/ieee754/ldbl-96/printf_fphex.c
index 798c788485..564f062727 100644
--- a/sysdeps/ieee754/ldbl-96/printf_fphex.c
+++ b/sysdeps/ieee754/ldbl-96/printf_fphex.c
@@ -1,5 +1,5 @@
 /* Print floating point number in hexadecimal notation according to ISO C99.
-   Copyright (C) 1997-2014 Free Software Foundation, Inc.
+   Copyright (C) 1997-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
diff --git a/sysdeps/ieee754/ldbl-96/s_asinhl.c b/sysdeps/ieee754/ldbl-96/s_asinhl.c
index 5e679bd719..75e47c71ea 100644
--- a/sysdeps/ieee754/ldbl-96/s_asinhl.c
+++ b/sysdeps/ieee754/ldbl-96/s_asinhl.c
@@ -29,6 +29,7 @@ static char rcsid[] = "$NetBSD: $";
  *		  := signl(x)*log1pl(|x| + x^2/(1 + sqrtl(1+x^2)))
  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -44,6 +45,11 @@ long double __asinhl(long double x)
 	GET_LDOUBLE_EXP(hx,x);
 	ix = hx&0x7fff;
 	if(__builtin_expect(ix< 0x3fde, 0)) {	/* |x|<2**-34 */
+	    if (fabsl (x) < LDBL_MIN)
+	      {
+		long double force_underflow = x * x;
+		math_force_eval (force_underflow);
+	      }
 	    if(huge+x>one) return x;	/* return x inexact except 0 */
 	}
 	if(__builtin_expect(ix>0x4020,0)) {		/* |x| > 2**34 */
diff --git a/sysdeps/ieee754/ldbl-96/s_cbrtl.c b/sysdeps/ieee754/ldbl-96/s_cbrtl.c
index 49a6891642..bc3c2abfeb 100644
--- a/sysdeps/ieee754/ldbl-96/s_cbrtl.c
+++ b/sysdeps/ieee754/ldbl-96/s_cbrtl.c
@@ -1,5 +1,5 @@
 /* Compute cubic root of double value.
-   Copyright (C) 1997-2014 Free Software Foundation, Inc.
+   Copyright (C) 1997-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Dirk Alboth <dirka@uni-paderborn.de> and
    Ulrich Drepper <drepper@cygnus.com>, 1997.
diff --git a/sysdeps/ieee754/ldbl-96/s_erfl.c b/sysdeps/ieee754/ldbl-96/s_erfl.c
index 47e4b9e909..c27de812ca 100644
--- a/sysdeps/ieee754/ldbl-96/s_erfl.c
+++ b/sysdeps/ieee754/ldbl-96/s_erfl.c
@@ -105,6 +105,7 @@
 
 
 #include <errno.h>
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -120,8 +121,6 @@ tiny = 1e-4931L,
  */
   /* 2/sqrt(pi) - 1 */
   efx = 1.2837916709551257389615890312154517168810E-1L,
-  /* 8 * (2/sqrt(pi) - 1) */
-  efx8 = 1.0270333367641005911692712249723613735048E0L,
 
   pp[6] = {
     1.122751350964552113068262337278335028553E6L,
@@ -272,7 +271,16 @@ __erfl (long double x)
       if (ix < 0x3fde8000) /* |x|<2**-33 */
 	{
 	  if (ix < 0x00080000)
-	    return 0.125 * (8.0 * x + efx8 * x);	/*avoid underflow */
+	    {
+	      /* Avoid spurious underflow.  */
+	      long double ret = 0.0625 * (16.0 * x + (16.0 * efx) * x);
+	      if (fabsl (ret) < LDBL_MIN)
+		{
+		  long double force_underflow = ret * ret;
+		  math_force_eval (force_underflow);
+		}
+	      return ret;
+	    }
 	  return x + efx * x;
 	}
       z = x * x;
diff --git a/sysdeps/ieee754/ldbl-96/s_fma.c b/sysdeps/ieee754/ldbl-96/s_fma.c
index fde2811040..8fd238cfe4 100644
--- a/sysdeps/ieee754/ldbl-96/s_fma.c
+++ b/sysdeps/ieee754/ldbl-96/s_fma.c
@@ -1,5 +1,5 @@
 /* Compute x * y + z as ternary operation.
-   Copyright (C) 2010-2014 Free Software Foundation, Inc.
+   Copyright (C) 2010-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Jakub Jelinek <jakub@redhat.com>, 2010.
 
@@ -21,6 +21,7 @@
 #include <math.h>
 #include <fenv.h>
 #include <ieee754.h>
+#include <math_private.h>
 
 /* This implementation uses rounding to odd to avoid problems with
    double rounding.  See a paper by Boldo and Melquiond:
@@ -29,17 +30,17 @@
 double
 __fma (double x, double y, double z)
 {
-  if (__builtin_expect (isinf (z), 0))
+  if (__glibc_unlikely (isinf (z)))
     {
       /* If z is Inf, but x and y are finite, the result should be
 	 z rather than NaN.  */
-      if (finite (x) && finite (y))
+      if (isfinite (x) && isfinite (y))
 	return (z + x) + y;
       return (x * y) + z;
     }
 
   /* Ensure correct sign of exact 0 + 0.  */
-  if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0))
+  if (__glibc_unlikely ((x == 0 || y == 0) && z == 0))
     return x * y + z;
 
   fenv_t env;
@@ -64,16 +65,17 @@ __fma (double x, double y, double z)
   t1 = m1 - t1;
   t2 = z - t2;
   long double a2 = t1 + t2;
+  /* Ensure the arithmetic is not scheduled after feclearexcept call.  */
+  math_force_eval (m2);
+  math_force_eval (a2);
   feclearexcept (FE_INEXACT);
 
-  /* If the result is an exact zero, ensure it has the correct
-     sign.  */
+  /* If the result is an exact zero, ensure it has the correct sign.  */
   if (a1 == 0 && m2 == 0)
     {
       feupdateenv (&env);
-      /* Ensure that round-to-nearest value of z + m1 is not
-	 reused.  */
-      asm volatile ("" : "=m" (z) : "m" (z));
+      /* Ensure that round-to-nearest value of z + m1 is not reused.  */
+      z = math_opt_barrier (z);
       return z + m1;
     }
 
diff --git a/sysdeps/ieee754/ldbl-96/s_fmal.c b/sysdeps/ieee754/ldbl-96/s_fmal.c
index 0564321354..eec5a02e0b 100644
--- a/sysdeps/ieee754/ldbl-96/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-96/s_fmal.c
@@ -1,5 +1,5 @@
 /* Compute x * y + z as ternary operation.
-   Copyright (C) 2010-2014 Free Software Foundation, Inc.
+   Copyright (C) 2010-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Jakub Jelinek <jakub@redhat.com>, 2010.
 
@@ -176,7 +176,7 @@ __fmal (long double x, long double y, long double z)
     }
 
   /* Ensure correct sign of exact 0 + 0.  */
-  if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0))
+  if (__glibc_unlikely ((x == 0 || y == 0) && z == 0))
     return x * y + z;
 
   fenv_t env;
@@ -201,16 +201,17 @@ __fmal (long double x, long double y, long double z)
   t1 = m1 - t1;
   t2 = z - t2;
   long double a2 = t1 + t2;
+  /* Ensure the arithmetic is not scheduled after feclearexcept call.  */
+  math_force_eval (m2);
+  math_force_eval (a2);
   feclearexcept (FE_INEXACT);
 
-  /* If the result is an exact zero, ensure it has the correct
-     sign.  */
+  /* If the result is an exact zero, ensure it has the correct sign.  */
   if (a1 == 0 && m2 == 0)
     {
       feupdateenv (&env);
-      /* Ensure that round-to-nearest value of z + m1 is not
-	 reused.  */
-      asm volatile ("" : "=m" (z) : "m" (z));
+      /* Ensure that round-to-nearest value of z + m1 is not reused.  */
+      z = math_opt_barrier (z);
       return z + m1;
     }
 
@@ -218,7 +219,7 @@ __fmal (long double x, long double y, long double z)
   /* Perform m2 + a2 addition with round to odd.  */
   u.d = a2 + m2;
 
-  if (__builtin_expect (adjust == 0, 1))
+  if (__glibc_likely (adjust == 0))
     {
       if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7fff)
 	u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
@@ -226,7 +227,7 @@ __fmal (long double x, long double y, long double z)
       /* Result is a1 + u.d.  */
       return a1 + u.d;
     }
-  else if (__builtin_expect (adjust > 0, 1))
+  else if (__glibc_likely (adjust > 0))
     {
       if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7fff)
 	u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
diff --git a/sysdeps/ieee754/ldbl-96/s_issignalingl.c b/sysdeps/ieee754/ldbl-96/s_issignalingl.c
index bbb98735f5..adbee9823c 100644
--- a/sysdeps/ieee754/ldbl-96/s_issignalingl.c
+++ b/sysdeps/ieee754/ldbl-96/s_issignalingl.c
@@ -1,5 +1,5 @@
 /* Test for signaling NaN.
-   Copyright (C) 2013-2014 Free Software Foundation, Inc.
+   Copyright (C) 2013-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
diff --git a/sysdeps/ieee754/ldbl-96/s_llrintl.c b/sysdeps/ieee754/ldbl-96/s_llrintl.c
index 0102f18242..070246cb6a 100644
--- a/sysdeps/ieee754/ldbl-96/s_llrintl.c
+++ b/sysdeps/ieee754/ldbl-96/s_llrintl.c
@@ -1,6 +1,6 @@
 /* Round argument to nearest integral value according to current rounding
    direction.
-   Copyright (C) 1997-2014 Free Software Foundation, Inc.
+   Copyright (C) 1997-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
diff --git a/sysdeps/ieee754/ldbl-96/s_llroundl.c b/sysdeps/ieee754/ldbl-96/s_llroundl.c
index 1559540a7a..8381649269 100644
--- a/sysdeps/ieee754/ldbl-96/s_llroundl.c
+++ b/sysdeps/ieee754/ldbl-96/s_llroundl.c
@@ -1,5 +1,5 @@
 /* Round long double value to long long int.
-   Copyright (C) 1997-2014 Free Software Foundation, Inc.
+   Copyright (C) 1997-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
diff --git a/sysdeps/ieee754/ldbl-96/s_lrintl.c b/sysdeps/ieee754/ldbl-96/s_lrintl.c
index a668934468..96c6a5ea7d 100644
--- a/sysdeps/ieee754/ldbl-96/s_lrintl.c
+++ b/sysdeps/ieee754/ldbl-96/s_lrintl.c
@@ -1,6 +1,6 @@
 /* Round argument to nearest integral value according to current rounding
    direction.
-   Copyright (C) 1997-2014 Free Software Foundation, Inc.
+   Copyright (C) 1997-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
diff --git a/sysdeps/ieee754/ldbl-96/s_lroundl.c b/sysdeps/ieee754/ldbl-96/s_lroundl.c
index a5c08d6391..e34edd307d 100644
--- a/sysdeps/ieee754/ldbl-96/s_lroundl.c
+++ b/sysdeps/ieee754/ldbl-96/s_lroundl.c
@@ -1,5 +1,5 @@
 /* Round long double value to long int.
-   Copyright (C) 1997-2014 Free Software Foundation, Inc.
+   Copyright (C) 1997-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
diff --git a/sysdeps/ieee754/ldbl-96/s_remquol.c b/sysdeps/ieee754/ldbl-96/s_remquol.c
index 866935c6e7..230d0fa57c 100644
--- a/sysdeps/ieee754/ldbl-96/s_remquol.c
+++ b/sysdeps/ieee754/ldbl-96/s_remquol.c
@@ -1,5 +1,5 @@
 /* Compute remainder and a congruent to the quotient.
-   Copyright (C) 1997-2014 Free Software Foundation, Inc.
+   Copyright (C) 1997-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -44,7 +44,7 @@ __remquol (long double x, long double p, int *quo)
     return (x * p) / (x * p); 			/* p = 0 */
   if ((ex == 0x7fff)				/* x not finite */
       || ((ep == 0x7fff)			/* p is NaN */
-	  && ((hp | lp) != 0)))
+	  && (((hp & 0x7fffffff) | lp) != 0)))
     return (x * p) / (x * p);
 
   if (ep <= 0x7ffb)
@@ -60,12 +60,12 @@ __remquol (long double x, long double p, int *quo)
   p  = fabsl (p);
   cquo = 0;
 
-  if (x >= 4 * p)
+  if (ep <= 0x7ffc && x >= 4 * p)
     {
       x -= 4 * p;
       cquo += 4;
     }
-  if (x >= 2 * p)
+  if (ep <= 0x7ffd && x >= 2 * p)
     {
       x -= 2 * p;
       cquo += 2;
@@ -101,6 +101,9 @@ __remquol (long double x, long double p, int *quo)
 
   *quo = qs ? -cquo : cquo;
 
+  /* Ensure correct sign of zero result in round-downward mode.  */
+  if (x == 0.0L)
+    x = 0.0L;
   if (sx)
     x = -x;
   return x;
diff --git a/sysdeps/ieee754/ldbl-96/s_roundl.c b/sysdeps/ieee754/ldbl-96/s_roundl.c
index e16e9c6f65..560f7afb99 100644
--- a/sysdeps/ieee754/ldbl-96/s_roundl.c
+++ b/sysdeps/ieee754/ldbl-96/s_roundl.c
@@ -1,5 +1,5 @@
 /* Round long double to integer away from zero.
-   Copyright (C) 1997-2014 Free Software Foundation, Inc.
+   Copyright (C) 1997-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
diff --git a/sysdeps/ieee754/ldbl-96/s_scalblnl.c b/sysdeps/ieee754/ldbl-96/s_scalblnl.c
index 76a4c59e6b..457e999c6c 100644
--- a/sysdeps/ieee754/ldbl-96/s_scalblnl.c
+++ b/sysdeps/ieee754/ldbl-96/s_scalblnl.c
@@ -25,8 +25,8 @@
 #include <math_private.h>
 
 static const long double
-two63   =  4.50359962737049600000e+15,
-twom63  =  1.08420217248550443400e-19,
+two63   =  0x1p63L,
+twom64  =  0x1p-64L,
 huge   = 1.0e+4900L,
 tiny   = 1.0e-4900L;
 
@@ -40,7 +40,7 @@ __scalblnl (long double x, long int n)
 	    if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
 	    x *= two63;
 	    GET_LDOUBLE_EXP(es,x);
-	    k = (hx&0x7fff) - 63;
+	    k = (es&0x7fff) - 63;
 	    }
 	if (__builtin_expect(k==0x7fff, 0)) return x+x;	/* NaN or Inf */
 	if (__builtin_expect(n< -50000, 0))
@@ -52,10 +52,9 @@ __scalblnl (long double x, long int n)
 	k = k+n;
 	if (__builtin_expect(k > 0, 1))		/* normal result */
 	    {SET_LDOUBLE_EXP(x,(es&0x8000)|k); return x;}
-	if (k <= -63)
+	if (k <= -64)
 	    return tiny*__copysignl(tiny,x); 	/*underflow*/
-	k += 63;				/* subnormal result */
+	k += 64;				/* subnormal result */
 	SET_LDOUBLE_EXP(x,(es&0x8000)|k);
-	return x*twom63;
+	return x*twom64;
 }
-weak_alias (__scalblnl, scalblnl)
diff --git a/sysdeps/ieee754/ldbl-96/s_signbitl.c b/sysdeps/ieee754/ldbl-96/s_signbitl.c
index 8256d9ddb0..bbe72a62ea 100644
--- a/sysdeps/ieee754/ldbl-96/s_signbitl.c
+++ b/sysdeps/ieee754/ldbl-96/s_signbitl.c
@@ -1,5 +1,5 @@
 /* Return nonzero value if number is negative.
-   Copyright (C) 1997-2014 Free Software Foundation, Inc.
+   Copyright (C) 1997-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
diff --git a/sysdeps/ieee754/ldbl-96/s_sincosl.c b/sysdeps/ieee754/ldbl-96/s_sincosl.c
index af8f899c24..37067bdaeb 100644
--- a/sysdeps/ieee754/ldbl-96/s_sincosl.c
+++ b/sysdeps/ieee754/ldbl-96/s_sincosl.c
@@ -1,5 +1,5 @@
 /* Compute sine and cosine of argument.
-   Copyright (C) 1997-2014 Free Software Foundation, Inc.
+   Copyright (C) 1997-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -17,6 +17,7 @@
    License along with the GNU C Library; if not, see
    <http://www.gnu.org/licenses/>.  */
 
+#include <errno.h>
 #include <math.h>
 
 #include <math_private.h>
@@ -41,6 +42,8 @@ __sincosl (long double x, long double *sinx, long double *cosx)
     {
       /* sin(Inf or NaN) is NaN */
       *sinx = *cosx = x - x;
+      if (__isinf_nsl (x))
+	__set_errno (EDOM);
     }
   else
     {
diff --git a/sysdeps/ieee754/ldbl-96/strtold_l.c b/sysdeps/ieee754/ldbl-96/strtold_l.c
index 82c33e52d6..c082e7472e 100644
--- a/sysdeps/ieee754/ldbl-96/strtold_l.c
+++ b/sysdeps/ieee754/ldbl-96/strtold_l.c
@@ -1,4 +1,4 @@
-/* Copyright (C) 1999-2014 Free Software Foundation, Inc.
+/* Copyright (C) 1999-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
diff --git a/sysdeps/ieee754/ldbl-96/t_sincosl.c b/sysdeps/ieee754/ldbl-96/t_sincosl.c
index da8fcfa112..4efc4d9a02 100644
--- a/sysdeps/ieee754/ldbl-96/t_sincosl.c
+++ b/sysdeps/ieee754/ldbl-96/t_sincosl.c
@@ -1,5 +1,5 @@
 /* Extended-precision floating point sine and cosine tables.
-   Copyright (C) 1999-2014 Free Software Foundation, Inc.
+   Copyright (C) 1999-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Based on quad-precision tables by Jakub Jelinek <jj@ultra.linux.cz>
 
diff --git a/sysdeps/ieee754/ldbl-96/w_expl.c b/sysdeps/ieee754/ldbl-96/w_expl.c
index af221184d1..3cc165b39d 100644
--- a/sysdeps/ieee754/ldbl-96/w_expl.c
+++ b/sysdeps/ieee754/ldbl-96/w_expl.c
@@ -1,4 +1,4 @@
-/* Copyright (C) 2011-2014 Free Software Foundation, Inc.
+/* Copyright (C) 2011-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
 
@@ -24,9 +24,9 @@ long double
 __expl (long double x)
 {
   long double z = __ieee754_expl (x);
-  if (__builtin_expect (!__finitel (z) || z == 0, 0)
-      && __finitel (x) && _LIB_VERSION != _IEEE_)
-    return __kernel_standard_l (x, x, 206 + !!__signbitl (x));
+  if (__builtin_expect (!isfinite (z) || z == 0, 0)
+      && isfinite (x) && _LIB_VERSION != _IEEE_)
+    return __kernel_standard_l (x, x, 206 + !!signbit (x));
 
   return z;
 }
diff --git a/sysdeps/ieee754/ldbl-96/x2y2m1.c b/sysdeps/ieee754/ldbl-96/x2y2m1.c
index 6ae2a63cbe..a6cc82cb29 100644
--- a/sysdeps/ieee754/ldbl-96/x2y2m1.c
+++ b/sysdeps/ieee754/ldbl-96/x2y2m1.c
@@ -1,5 +1,5 @@
 /* Compute x^2 + y^2 - 1, without large cancellation error.
-   Copyright (C) 2012-2014 Free Software Foundation, Inc.
+   Copyright (C) 2012-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
diff --git a/sysdeps/ieee754/ldbl-96/x2y2m1l.c b/sysdeps/ieee754/ldbl-96/x2y2m1l.c
index 575bcd812d..11757c6af8 100644
--- a/sysdeps/ieee754/ldbl-96/x2y2m1l.c
+++ b/sysdeps/ieee754/ldbl-96/x2y2m1l.c
@@ -1,5 +1,5 @@
 /* Compute x^2 + y^2 - 1, without large cancellation error.
-   Copyright (C) 2012-2014 Free Software Foundation, Inc.
+   Copyright (C) 2012-2015 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or