about summary refs log tree commit diff
path: root/sysdeps
diff options
context:
space:
mode:
authorAdhemerval Zanella <azanella@linux.vnet.ibm.com>2012-06-05 10:42:49 -0300
committerAdhemerval Zanella <azanella@linux.vnet.ibm.com>2012-06-05 10:42:49 -0300
commit6043738b3692bd542fc47596628aefc8dba75523 (patch)
tree5a2b8b4bfe0f0d9e34021e408f9d04a33c8549a6 /sysdeps
parent34ae0b3270c67cae0c54ac98b693fdf7d010a206 (diff)
downloadglibc-6043738b3692bd542fc47596628aefc8dba75523.tar.gz
glibc-6043738b3692bd542fc47596628aefc8dba75523.tar.xz
glibc-6043738b3692bd542fc47596628aefc8dba75523.zip
Fix spurious undeflow for ldbl-128ibm erfl
For values higher than 25.6283 erflc underflow, so adjust erfl
to return a constant value based argument sign.
Diffstat (limited to 'sysdeps')
-rw-r--r--sysdeps/ieee754/ldbl-128ibm/s_erfl.c21
1 files changed, 18 insertions, 3 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c
index f91a00ff5d..6a4475ed6b 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c
@@ -72,6 +72,8 @@
  *              z=1/x^2
  *         The interval is partitioned into several segments
  *         of width 1/8 in 1/x.
+ *	erf(x) = 1.0 - erfc(x) if x < 25.6283 else
+ *	erf(x) = sign(x)*(1.0 - tiny)
  *
  *      Note1:
  *	   To compute exp(-x*x-0.5625+R/S), let s be a single
@@ -85,6 +87,9 @@
  *		erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
  *			  x*sqrt(pi)
  *
+ *      Note3:
+ * 	   For x higher than 25.6283, erf(x) underflows.
+ *
  *      5. For inf > x >= 107
  *	erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
  *	erfc(x) = tiny*tiny (raise underflow) if x > 0
@@ -770,9 +775,19 @@ __erfl (long double x)
 
   if (ix >= 0x3ff00000) /* |x| >= 1.0 */
     {
-      y = __erfcl (x);
-      return (one - y);
-      /*    return (one - __erfcl (x)); */
+      if (ix >= 0x4039A0DE)
+	{
+	/* __erfcl (x) underflows if x > 25.6283 */
+	  if (sign)
+	    return one-tiny;
+	  else
+	    return tiny-one;
+	}
+      else
+	{
+	  y = __erfcl (x);
+	  return (one - y);
+	}
     }
   u.parts32.w0 = ix;
   a = u.value;