5 files changed, 54 insertions, 17 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_erf.c b/sysdeps/ieee754/dbl-64/s_erf.c
index 3f37397224..ea0a73424e 100644
--- a/sysdeps/ieee754/dbl-64/s_erf.c
+++ b/sysdeps/ieee754/dbl-64/s_erf.c
@@ -128,7 +128,6 @@ static const double
  * Coefficients for approximation to  erf on [0,0.84375]
  */
   efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
-  efx8 = 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
   pp[] = { 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
 	   -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
 	   -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
@@ -211,7 +210,16 @@ __erf (double x)
       if (ix < 0x3e300000)              /* |x|<2**-28 */
 	{
 	  if (ix < 0x00800000)
-	    return 0.125 * (8.0 * x + efx8 * x);  /*avoid underflow */
+	    {
+	      /* Avoid spurious underflow.  */
+	      double ret = 0.0625 * (16.0 * x + (16.0 * efx) * x);
+	      if (fabs (ret) < DBL_MIN)
+		{
+		  double force_underflow = ret * ret;
+		  math_force_eval (force_underflow);
+		}
+	      return ret;
+	    }
 	  return x + efx * x;
 	}
       z = x * x;
diff --git a/sysdeps/ieee754/flt-32/s_erff.c b/sysdeps/ieee754/flt-32/s_erff.c
index 7c09589648..d23cff9293 100644
--- a/sysdeps/ieee754/flt-32/s_erff.c
+++ b/sysdeps/ieee754/flt-32/s_erff.c
@@ -33,7 +33,6 @@ erx =  8.4506291151e-01, /* 0x3f58560b */
  * Coefficients for approximation to  erf on [0,0.84375]
  */
 efx =  1.2837916613e-01, /* 0x3e0375d4 */
-efx8=  1.0270333290e+00, /* 0x3f8375d4 */
 pp0  =  1.2837916613e-01, /* 0x3e0375d4 */
 pp1  = -3.2504209876e-01, /* 0xbea66beb */
 pp2  = -2.8481749818e-02, /* 0xbce9528f */
@@ -111,8 +110,16 @@ float __erff(float x)
 	if(ix < 0x3f580000) {		/* |x|<0.84375 */
 	    if(ix < 0x31800000) { 	/* |x|<2**-28 */
 	        if (ix < 0x04000000)
-		    /*avoid underflow */
-		    return (float)0.125*((float)8.0*x+efx8*x);
+		  {
+		    /* Avoid spurious underflow.  */
+		    float ret = 0.0625f * (16.0f * x + (16.0f * efx) * x);
+		    if (fabsf (ret) < FLT_MIN)
+		      {
+			float force_underflow = ret * ret;
+			math_force_eval (force_underflow);
+		      }
+		    return ret;
+		  }
 		return x + efx*x;
 	    }
 	    z = x*x;
diff --git a/sysdeps/ieee754/ldbl-128/s_erfl.c b/sysdeps/ieee754/ldbl-128/s_erfl.c
index ef65ed8922..35ca8c1360 100644
--- a/sysdeps/ieee754/ldbl-128/s_erfl.c
+++ b/sysdeps/ieee754/ldbl-128/s_erfl.c
@@ -97,6 +97,7 @@
  */
 
 #include <errno.h>
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -143,9 +144,7 @@ tiny = 1e-4931L,
   one = 1.0L,
   two = 2.0L,
   /* 2/sqrt(pi) - 1 */
-  efx = 1.2837916709551257389615890312154517168810E-1L,
-  /* 8 * (2/sqrt(pi) - 1) */
-  efx8 = 1.0270333367641005911692712249723613735048E0L;
+  efx = 1.2837916709551257389615890312154517168810E-1L;
 
 
 /* erf(x)  = x  + x R(x^2)
@@ -782,7 +781,16 @@ __erfl (long double x)
       if (ix < 0x3fc60000) /* |x|<2**-57 */
 	{
 	  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;
 	}
       y = a + a * neval (z, TN1, NTN1) / deval (z, TD1, NTD1);
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c
index 95dc415845..f55e8b7879 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c
@@ -102,6 +102,7 @@
  */
 
 #include <errno.h>
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 #include <math_ldbl_opt.h>
@@ -149,9 +150,7 @@ tiny = 1e-300L,
   one = 1.0L,
   two = 2.0L,
   /* 2/sqrt(pi) - 1 */
-  efx = 1.2837916709551257389615890312154517168810E-1L,
-  /* 8 * (2/sqrt(pi) - 1) */
-  efx8 = 1.0270333367641005911692712249723613735048E0L;
+  efx = 1.2837916709551257389615890312154517168810E-1L;
 
 
 /* erf(x)  = x  + x R(x^2)
@@ -801,10 +800,17 @@ __erfl (long double x)
 	  if (ix < 0x00800000)
 	    {
 	      /* erf (-0) = -0.  Unfortunately, for IBM extended double
-		 0.125 * (8.0 * x + efx8 * x) for x = -0 evaluates to 0.  */
+		 0.0625 * (16.0 * x + (16.0 * efx) * x) for x = -0
+		 evaluates to 0.  */
 	      if (x == 0)
 		return x;
-	      return 0.125 * (8.0 * x + efx8 * x);	/*avoid 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;
 	}
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;