diff options
| author | Samuel Thibault <samuel.thibault@ens-lyon.org> | 2016-08-20 19:50:45 +0200 |
|---|---|---|
| committer | Samuel Thibault <samuel.thibault@ens-lyon.org> | 2016-08-20 19:50:45 +0200 |
| commit | 4dd9e35bfd35d3138bc44169baba098005bad51e (patch) | |
| tree | a4939c43a9c3fe00eb27f023e14acc5e1fe8808c /sysdeps/ieee754/ldbl-128ibm/s_erfl.c | |
| parent | bd42a4599d1b6f77bcfe1e4f67b7cbd9e1cb2dfd (diff) | |
| parent | f76453c31593957fec1a99b986bfa5506618b79c (diff) | |
Merge commit 'refs/top-bases/t/bigmem' into t/bigmem
Diffstat (limited to 'sysdeps/ieee754/ldbl-128ibm/s_erfl.c')
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_erfl.c | 18 |
1 files changed, 12 insertions, 6 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c index 95dc415845..f6fcf48cfa 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c @@ -102,6 +102,7 @@ */ #include <errno.h> +#include <float.h> #include <math.h> #include <math_private.h> #include <math_ldbl_opt.h> @@ -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; } @@ -888,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; |