Fix erfc errno setting on underflow (bug 6786).

5 files changed, 65 insertions, 10 deletions

diff --git a/sysdeps/ieee754/dbl-64/s_erf.c b/sysdeps/ieee754/dbl-64/s_erf.c
index aab4ed5..3f37397 100644
--- a/sysdeps/ieee754/dbl-64/s_erf.c
+++ b/sysdeps/ieee754/dbl-64/s_erf.c
@@ -112,6 +112,8 @@ static char rcsid[] = "$NetBSD: s_erf.c,v 1.8 1995/05/10 20:47:05 jtc Exp $";
  */
 
 
+#include <errno.h>
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -391,14 +393,25 @@ __erfc (double x)
       r = __ieee754_exp (-z * z - 0.5625) *
 	  __ieee754_exp ((z - x) * (z + x) + R / S);
       if (hx > 0)
-	return r / x;
+	{
+#if FLT_EVAL_METHOD != 0
+	  volatile
+#endif
+	  double ret = r / x;
+	  if (ret == 0)
+	    __set_errno (ERANGE);
+	  return ret;
+	}
       else
 	return two - r / x;
     }
   else
     {
       if (hx > 0)
-	return tiny * tiny;
+	{
+	  __set_errno (ERANGE);
+	  return tiny * tiny;
+	}
       else
 	return two - tiny;
     }
diff --git a/sysdeps/ieee754/flt-32/s_erff.c b/sysdeps/ieee754/flt-32/s_erff.c
index 7d17f42..7c09589 100644
--- a/sysdeps/ieee754/flt-32/s_erff.c
+++ b/sysdeps/ieee754/flt-32/s_erff.c
@@ -17,6 +17,8 @@
 static char rcsid[] = "$NetBSD: s_erff.c,v 1.4 1995/05/10 20:47:07 jtc Exp $";
 #endif
 
+#include <errno.h>
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -203,9 +205,22 @@ float __erfcf(float x)
 	    SET_FLOAT_WORD(z,ix&0xffffe000);
 	    r  =  __ieee754_expf(-z*z-(float)0.5625)*
 			__ieee754_expf((z-x)*(z+x)+R/S);
-	    if(hx>0) return r/x; else return two-r/x;
+	    if(hx>0) {
+#if FLT_EVAL_METHOD != 0
+		volatile
+#endif
+		float ret = r/x;
+		if (ret == 0)
+		    __set_errno (ERANGE);
+		return ret;
+	    } else
+		return two-r/x;
 	} else {
-	    if(hx>0) return tiny*tiny; else return two-tiny;
+	    if(hx>0) {
+		__set_errno (ERANGE);
+		return tiny*tiny;
+	    } else
+		return two-tiny;
 	}
 }
 weak_alias (__erfcf, erfcf)
diff --git a/sysdeps/ieee754/ldbl-128/s_erfl.c b/sysdeps/ieee754/ldbl-128/s_erfl.c
index d9d4195..ef65ed8 100644
--- a/sysdeps/ieee754/ldbl-128/s_erfl.c
+++ b/sysdeps/ieee754/ldbl-128/s_erfl.c
@@ -96,6 +96,7 @@
  *		erfc/erf(NaN) is NaN
  */
 
+#include <errno.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -918,14 +919,22 @@ __erfcl (long double x)
       r = __ieee754_expl (-z * z - 0.5625) *
 	__ieee754_expl ((z - x) * (z + x) + p);
       if ((sign & 0x80000000) == 0)
-	return r / x;
+	{
+	  long double ret = r / x;
+	  if (ret == 0)
+	    __set_errno (ERANGE);
+	  return ret;
+	}
       else
 	return two - r / x;
     }
   else
     {
       if ((sign & 0x80000000) == 0)
-	return tiny * tiny;
+	{
+	  __set_errno (ERANGE);
+	  return tiny * tiny;
+	}
       else
 	return two - tiny;
     }
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c
index c861c65..95dc415 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c
@@ -101,6 +101,7 @@
  *		erfc/erf(NaN) is NaN
  */
 
+#include <errno.h>
 #include <math.h>
 #include <math_private.h>
 #include <math_ldbl_opt.h>
@@ -940,14 +941,22 @@ __erfcl (long double x)
       r = __ieee754_expl (-z * z - 0.5625) *
 	__ieee754_expl ((z - x) * (z + x) + p);
       if ((hx & 0x80000000) == 0)
-	return r / x;
+	{
+	  long double ret = r / x;
+	  if (ret == 0)
+	    __set_errno (ERANGE);
+	  return ret;
+	}
       else
 	return two - r / x;
     }
   else
     {
       if ((hx & 0x80000000) == 0)
-	return tiny * tiny;
+	{
+	  __set_errno (ERANGE);
+	  return tiny * tiny;
+	}
       else
 	return two - tiny;
     }
diff --git a/sysdeps/ieee754/ldbl-96/s_erfl.c b/sysdeps/ieee754/ldbl-96/s_erfl.c
index 17d2278..47e4b9e 100644
--- a/sysdeps/ieee754/ldbl-96/s_erfl.c
+++ b/sysdeps/ieee754/ldbl-96/s_erfl.c
@@ -104,6 +104,7 @@
  */
 
 
+#include <errno.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -422,14 +423,22 @@ __erfcl (long double x)
       r = __ieee754_expl (-z * z - 0.5625) *
 	__ieee754_expl ((z - x) * (z + x) + R / S);
       if ((se & 0x8000) == 0)
-	return r / x;
+	{
+	  long double ret = r / x;
+	  if (ret == 0)
+	    __set_errno (ERANGE);
+	  return ret;
+	}
       else
 	return two - r / x;
     }
   else
     {
       if ((se & 0x8000) == 0)
-	return tiny * tiny;
+	{
+	  __set_errno (ERANGE);
+	  return tiny * tiny;
+	}
       else
 	return two - tiny;
     }