From 0bf061d3e37e0f72c2378cbded67c58df4f58a4b Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Wed, 14 May 2014 12:34:03 +0000 Subject: Fix erf underflow handling near 0 (bug 16516). Bug 16516 reports spurious underflows from erf (for all floating-point types), when the result is close to underflowing but does not actually underflow. erf (x) is about (2/sqrt(pi))*x for x close to 0, so there are subnormal arguments for which it does not underflow. The various implementations do (x + efx*x) (for efx = 2/sqrt(pi) - 1), for greater accuracy than if just using a single multiplication by an approximation to 2/sqrt(pi) (effectively, this way there are a few more bits in the approximation to 2/sqrt(pi)). This can introduce underflows when efx*x underflows even though the final result does not, so a scaled calculation with 8*efx is done in these cases - but 8 is not a big enough scale factor to avoid all such underflows. 16 is (any underflows with a scale factor of 16 would only occur when the final result underflows), so this patch changes the code to use that factor. Rather than recomputing all the values of the efx8 variable, it is removed, leaving it to the compiler's constant folding to compute 16*efx. As such scaling can also lose underflows when the final scaling down happens to be exact, appropriate checks are added to ensure underflow exceptions occur when required in such cases. Tested x86_64 and x86; no ulps updates needed. Also spot-checked for powerpc32 and mips64 to verify the changes to the ldbl-128ibm and ldbl-128 implementations. [BZ #16516] * sysdeps/ieee754/dbl-64/s_erf.c (efx8): Remove variable. (__erf): Scale by 16 instead of 8 in potentially underflowing case. Ensure exception if result actually underflows. * sysdeps/ieee754/flt-32/s_erff.c (efx8): Remove variable. (__erff): Scale by 16 instead of 8 in potentially underflowing case. Ensure exception if result actually underflows. * sysdeps/ieee754/ldbl-128/s_erfl.c: Include . (efx8): Remove variable. (__erfl): Scale by 16 instead of 8 in potentially underflowing case. Ensure exception if result actually underflows. * sysdeps/ieee754/ldbl-128ibm/s_erfl.c: Include . (efx8): Remove variable. (__erfl): Scale by 16 instead of 8 in potentially underflowing case. Ensure exception if result actually underflows. * sysdeps/ieee754/ldbl-96/s_erfl.c: Include . (efx8): Remove variable. (__erfl): Scale by 16 instead of 8 in potentially underflowing case. Ensure exception if result actually underflows. * math/auto-libm-test-in: Add more tests of erf. * math/auto-libm-test-out: Regenerated. --- sysdeps/ieee754/ldbl-128ibm/s_erfl.c | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) (limited to 'sysdeps/ieee754/ldbl-128ibm/s_erfl.c') 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 +#include #include #include #include @@ -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; } -- cgit v1.2.3 From 31a8780d0b41934c6ab5160f834f56a2d5d97f4a Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Fri, 22 May 2015 17:48:45 +0000 Subject: Fix ldbl-128 / ldbl-128ibm erfcl for -Wuninitialized The ldbl-128 and ldbl-128ibm implementations of erfcl produce uninitialized variable warnings with -Wuninitialized because of switch statements where in fact one of the cases will always be executed, but the compiler does not see that these cases cover all possibilities (and because the reasoning that it does involves inequalities on the representation of a floating point value leading to a set of possible values for 8.0 times that value, converted to int, it's highly nontrivial for the compiler to see that). This patch fixes those warnings by converting the last case in those switch statements to a "default" case. Tested for powerpc and mips64. * sysdeps/ieee754/ldbl-128/s_erfl.c (__erfcl): Make case 9 in switch statement into default case. * sysdeps/ieee754/ldbl-128ibm/s_erfl.c (__erfcl): Likewise. --- ChangeLog | 4 ++++ sysdeps/ieee754/ldbl-128/s_erfl.c | 2 +- sysdeps/ieee754/ldbl-128ibm/s_erfl.c | 2 +- 3 files changed, 6 insertions(+), 2 deletions(-) (limited to 'sysdeps/ieee754/ldbl-128ibm/s_erfl.c') diff --git a/ChangeLog b/ChangeLog index 21c4fa68e0..5b31961d17 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,9 @@ 2015-05-22 Joseph Myers + * sysdeps/ieee754/ldbl-128/s_erfl.c (__erfcl): Make case 9 in + switch statement into default case. + * sysdeps/ieee754/ldbl-128ibm/s_erfl.c (__erfcl): Likewise. + * sysdeps/ieee754/ldbl-128/e_asinl.c (__ieee754_asinl): Don't use a conditional in forcing "inexact". * sysdeps/ieee754/ldbl-128ibm/e_asinl.c (__ieee754_asinl): diff --git a/sysdeps/ieee754/ldbl-128/s_erfl.c b/sysdeps/ieee754/ldbl-128/s_erfl.c index 9bfc9c580b..fa4609f136 100644 --- a/sysdeps/ieee754/ldbl-128/s_erfl.c +++ b/sysdeps/ieee754/ldbl-128/s_erfl.c @@ -874,7 +874,7 @@ __erfcl (long double x) y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19); y += C19a; break; - case 9: + default: /* i == 9. */ z = x - 1.125L; y = C20b + z * neval (z, RNr20, NRNr20) / deval (z, RDr20, NRDr20); y += C20a; diff --git a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c index f55e8b7879..f6fcf48cfa 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c @@ -894,7 +894,7 @@ __erfcl (long double x) y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19); y += C19a; break; - case 9: + default: /* i == 9. */ z = x - 1.125L; y = C20b + z * neval (z, RNr20, NRNr20) / deval (z, RDr20, NRDr20); y += C20a; -- cgit v1.2.3