From 6043738b3692bd542fc47596628aefc8dba75523 Mon Sep 17 00:00:00 2001 From: Adhemerval Zanella Date: Tue, 5 Jun 2012 10:42:49 -0300 Subject: 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. --- ChangeLog | 5 +++++ sysdeps/ieee754/ldbl-128ibm/s_erfl.c | 21 ++++++++++++++++++--- 2 files changed, 23 insertions(+), 3 deletions(-) diff --git a/ChangeLog b/ChangeLog index 925e986f1f..dc61b0dd89 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,8 @@ +2012-06-05 Adhemerval Zanella + + * sysdeps/ieee754/ldbl-128ibm/s_erfl.c: Fix spurious underflow for + values higher than 25.6283. + 2012-06-04 Adhemerval Zanella * sysdeps/ieee754/ldbl-128ibm/e_fmodl.c (__ieee754_fmodl): Fix 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; -- cgit 1.4.1