Fix erfcf spurious underflows (bug 18217).

The flt-32 implementation of erfcf produces spurious underflow
exceptions for some arguments close to 0, because of calculations
squaring the argument and then multiplying by small constants.  This
patch fixes this by adjusting the threshold for arguments for which
the result is so close to 1 that 1 - x will give the right result from
2**-56 to 2**-26.  (If 1 - x * 2/sqrt(pi) were used, the errors would be
on the order of x^3 and a much larger threshold could be used.)

Tested for x86_64 and x86.

	[BZ #18217]
	* sysdeps/ieee754/flt-32/s_erff.c (__erfcf): Use 2**-26 not 2**-56
	as threshold for returning 1 - x.
	* math/auto-libm-test-in: Add more tests of erfc.
	* math/auto-libm-test-out: Regenerated.

diff --git a/ChangeLog b/ChangeLog
index 1179142c92334c35d3dcb4c91965c91c82d5663c..454a8d26e172ecf09f1a7e24c7f2fa6e7db48ea7 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,11 @@
+2015-05-15  Joseph Myers  <joseph@codesourcery.com>
+
+       [BZ #18217]
+       * sysdeps/ieee754/flt-32/s_erff.c (__erfcf): Use 2**-26 not 2**-56
+       as threshold for returning 1 - x.
+       * math/auto-libm-test-in: Add more tests of erfc.
+       * math/auto-libm-test-out: Regenerated.
+
 2015-05-14  Joseph Myers  <joseph@codesourcery.com>
 
        [BZ #18196]

diff --git a/NEWS b/NEWS
index c46ee35103570e1efe79e3fd6bfb5703568d30a2..9251333195e54be29a506815c00ab0664ee1d25d 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -16,8 +16,8 @@ Version 2.22
   17964, 17965, 17967, 17969, 17978, 17987, 17991, 17996, 17998, 17999,
   18007, 18019, 18020, 18029, 18030, 18032, 18036, 18038, 18039, 18042,
   18043, 18046, 18047, 18068, 18080, 18093, 18100, 18104, 18110, 18111,
-  18125, 18128, 18138, 18185, 18196, 18197, 18206, 18210, 18211, 18247,
-  18287, 18319, 18333, 18346, 18397, 18409.
+  18125, 18128, 18138, 18185, 18196, 18197, 18206, 18210, 18211, 18217,
+  18247, 18287, 18319, 18333, 18346, 18397, 18409.
 
 * Cache information can be queried via sysconf() function on s390 e.g. with
   _SC_LEVEL1_ICACHE_SIZE as argument.

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 1deb99ad0dea73d621e4cd5b3629a819b86e4b59..c4bfe7439c1fe030052af29717535a80c4035963 100644 (file)
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -1129,6 +1129,8 @@ erf -0x1.3a0d48p+0
 
 erfc 0.0
 erfc -0
+erfc 0x1p-55
+erfc -0x1p-55
 erfc 0.125
 erfc 0.75
 erfc 1.25

diff --git a/math/auto-libm-test-out b/math/auto-libm-test-out
index 9f6bbbbf6f00f2a79af7f156d205f550d34e0680..d717ea4431dcd8290abbc585beba334c39ec083c 100644 (file)
--- a/math/auto-libm-test-out
+++ b/math/auto-libm-test-out
@@ -101555,6 +101555,56 @@ erfc -0
 = erfc tonearest ldbl-128ibm -0x0p+0L : 0x1p+0L : inexact-ok
 = erfc towardzero ldbl-128ibm -0x0p+0L : 0x1p+0L : inexact-ok
 = erfc upward ldbl-128ibm -0x0p+0L : 0x1p+0L : inexact-ok
+erfc 0x1p-55
+= erfc downward flt-32 0x2p-56f : 0xf.fffffp-4f : inexact-ok
+= erfc tonearest flt-32 0x2p-56f : 0x1p+0f : inexact-ok
+= erfc towardzero flt-32 0x2p-56f : 0xf.fffffp-4f : inexact-ok
+= erfc upward flt-32 0x2p-56f : 0x1p+0f : inexact-ok
+= erfc downward dbl-64 0x2p-56 : 0xf.ffffffffffff8p-4 : inexact-ok
+= erfc tonearest dbl-64 0x2p-56 : 0x1p+0 : inexact-ok
+= erfc towardzero dbl-64 0x2p-56 : 0xf.ffffffffffff8p-4 : inexact-ok
+= erfc upward dbl-64 0x2p-56 : 0x1p+0 : inexact-ok
+= erfc downward ldbl-96-intel 0x2p-56L : 0xf.ffffffffffffdbep-4L : inexact-ok
+= erfc tonearest ldbl-96-intel 0x2p-56L : 0xf.ffffffffffffdbep-4L : inexact-ok
+= erfc towardzero ldbl-96-intel 0x2p-56L : 0xf.ffffffffffffdbep-4L : inexact-ok
+= erfc upward ldbl-96-intel 0x2p-56L : 0xf.ffffffffffffdbfp-4L : inexact-ok
+= erfc downward ldbl-96-m68k 0x2p-56L : 0xf.ffffffffffffdbep-4L : inexact-ok
+= erfc tonearest ldbl-96-m68k 0x2p-56L : 0xf.ffffffffffffdbep-4L : inexact-ok
+= erfc towardzero ldbl-96-m68k 0x2p-56L : 0xf.ffffffffffffdbep-4L : inexact-ok
+= erfc upward ldbl-96-m68k 0x2p-56L : 0xf.ffffffffffffdbfp-4L : inexact-ok
+= erfc downward ldbl-128 0x2p-56L : 0xf.ffffffffffffdbe4515f7ac925d8p-4L : inexact-ok
+= erfc tonearest ldbl-128 0x2p-56L : 0xf.ffffffffffffdbe4515f7ac925ep-4L : inexact-ok
+= erfc towardzero ldbl-128 0x2p-56L : 0xf.ffffffffffffdbe4515f7ac925d8p-4L : inexact-ok
+= erfc upward ldbl-128 0x2p-56L : 0xf.ffffffffffffdbe4515f7ac925ep-4L : inexact-ok
+= erfc downward ldbl-128ibm 0x2p-56L : 0xf.ffffffffffffdbe4515f7ac924p-4L : inexact-ok
+= erfc tonearest ldbl-128ibm 0x2p-56L : 0xf.ffffffffffffdbe4515f7ac924p-4L : inexact-ok
+= erfc towardzero ldbl-128ibm 0x2p-56L : 0xf.ffffffffffffdbe4515f7ac924p-4L : inexact-ok
+= erfc upward ldbl-128ibm 0x2p-56L : 0xf.ffffffffffffdbe4515f7ac928p-4L : inexact-ok
+erfc -0x1p-55
+= erfc downward flt-32 -0x2p-56f : 0x1p+0f : inexact-ok
+= erfc tonearest flt-32 -0x2p-56f : 0x1p+0f : inexact-ok
+= erfc towardzero flt-32 -0x2p-56f : 0x1p+0f : inexact-ok
+= erfc upward flt-32 -0x2p-56f : 0x1.000002p+0f : inexact-ok
+= erfc downward dbl-64 -0x2p-56 : 0x1p+0 : inexact-ok
+= erfc tonearest dbl-64 -0x2p-56 : 0x1p+0 : inexact-ok
+= erfc towardzero dbl-64 -0x2p-56 : 0x1p+0 : inexact-ok
+= erfc upward dbl-64 -0x2p-56 : 0x1.0000000000001p+0 : inexact-ok
+= erfc downward ldbl-96-intel -0x2p-56L : 0x1.000000000000024p+0L : inexact-ok
+= erfc tonearest ldbl-96-intel -0x2p-56L : 0x1.0000000000000242p+0L : inexact-ok
+= erfc towardzero ldbl-96-intel -0x2p-56L : 0x1.000000000000024p+0L : inexact-ok
+= erfc upward ldbl-96-intel -0x2p-56L : 0x1.0000000000000242p+0L : inexact-ok
+= erfc downward ldbl-96-m68k -0x2p-56L : 0x1.000000000000024p+0L : inexact-ok
+= erfc tonearest ldbl-96-m68k -0x2p-56L : 0x1.0000000000000242p+0L : inexact-ok
+= erfc towardzero ldbl-96-m68k -0x2p-56L : 0x1.000000000000024p+0L : inexact-ok
+= erfc upward ldbl-96-m68k -0x2p-56L : 0x1.0000000000000242p+0L : inexact-ok
+= erfc downward ldbl-128 -0x2p-56L : 0x1.0000000000000241baea08536da2p+0L : inexact-ok
+= erfc tonearest ldbl-128 -0x2p-56L : 0x1.0000000000000241baea08536da2p+0L : inexact-ok
+= erfc towardzero ldbl-128 -0x2p-56L : 0x1.0000000000000241baea08536da2p+0L : inexact-ok
+= erfc upward ldbl-128 -0x2p-56L : 0x1.0000000000000241baea08536da3p+0L : inexact-ok
+= erfc downward ldbl-128ibm -0x2p-56L : 0x1.0000000000000241baea08536d8p+0L : inexact-ok
+= erfc tonearest ldbl-128ibm -0x2p-56L : 0x1.0000000000000241baea08536d8p+0L : inexact-ok
+= erfc towardzero ldbl-128ibm -0x2p-56L : 0x1.0000000000000241baea08536d8p+0L : inexact-ok
+= erfc upward ldbl-128ibm -0x2p-56L : 0x1.0000000000000241baea08536ep+0L : inexact-ok
 erfc 0.125
 = erfc downward flt-32 0x2p-4f : 0xd.c143cp-4f : inexact-ok
 = erfc tonearest flt-32 0x2p-4f : 0xd.c143dp-4f : inexact-ok

diff --git a/sysdeps/ieee754/flt-32/s_erff.c b/sysdeps/ieee754/flt-32/s_erff.c
index d23cff9293894a5ebe51a7925a5eb34b579041f7..2be44cc40cbdc84c8d4f09535959399741e718dd 100644 (file)
--- a/sysdeps/ieee754/flt-32/s_erff.c
+++ b/sysdeps/ieee754/flt-32/s_erff.c
@@ -169,7 +169,7 @@ float __erfcf(float x)
        }
 
        if(ix < 0x3f580000) {           /* |x|<0.84375 */
-           if(ix < 0x23800000)         /* |x|<2**-56 */
+           if(ix < 0x32800000)         /* |x|<2**-26 */
                return one-x;
            z = x*x;
            r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));