diff options
Diffstat (limited to 'sysdeps')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_erf.c | 12 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/s_erff.c | 13 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128/s_erfl.c | 16 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_erfl.c | 16 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-96/s_erfl.c | 14 |
5 files changed, 54 insertions, 17 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_erf.c b/sysdeps/ieee754/dbl-64/s_erf.c index 3f37397224..ea0a73424e 100644 --- a/sysdeps/ieee754/dbl-64/s_erf.c +++ b/sysdeps/ieee754/dbl-64/s_erf.c @@ -128,7 +128,6 @@ static const double * Coefficients for approximation to erf on [0,0.84375] */ efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */ - efx8 = 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */ pp[] = { 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */ -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */ -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */ @@ -211,7 +210,16 @@ __erf (double x) if (ix < 0x3e300000) /* |x|<2**-28 */ { if (ix < 0x00800000) - return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */ + { + /* Avoid spurious underflow. */ + double ret = 0.0625 * (16.0 * x + (16.0 * efx) * x); + if (fabs (ret) < DBL_MIN) + { + double force_underflow = ret * ret; + math_force_eval (force_underflow); + } + return ret; + } return x + efx * x; } z = x * x; diff --git a/sysdeps/ieee754/flt-32/s_erff.c b/sysdeps/ieee754/flt-32/s_erff.c index 7c09589648..d23cff9293 100644 --- a/sysdeps/ieee754/flt-32/s_erff.c +++ b/sysdeps/ieee754/flt-32/s_erff.c @@ -33,7 +33,6 @@ erx = 8.4506291151e-01, /* 0x3f58560b */ * Coefficients for approximation to erf on [0,0.84375] */ efx = 1.2837916613e-01, /* 0x3e0375d4 */ -efx8= 1.0270333290e+00, /* 0x3f8375d4 */ pp0 = 1.2837916613e-01, /* 0x3e0375d4 */ pp1 = -3.2504209876e-01, /* 0xbea66beb */ pp2 = -2.8481749818e-02, /* 0xbce9528f */ @@ -111,8 +110,16 @@ float __erff(float x) if(ix < 0x3f580000) { /* |x|<0.84375 */ if(ix < 0x31800000) { /* |x|<2**-28 */ if (ix < 0x04000000) - /*avoid underflow */ - return (float)0.125*((float)8.0*x+efx8*x); + { + /* Avoid spurious underflow. */ + float ret = 0.0625f * (16.0f * x + (16.0f * efx) * x); + if (fabsf (ret) < FLT_MIN) + { + float force_underflow = ret * ret; + math_force_eval (force_underflow); + } + return ret; + } return x + efx*x; } z = x*x; diff --git a/sysdeps/ieee754/ldbl-128/s_erfl.c b/sysdeps/ieee754/ldbl-128/s_erfl.c index ef65ed8922..35ca8c1360 100644 --- a/sysdeps/ieee754/ldbl-128/s_erfl.c +++ b/sysdeps/ieee754/ldbl-128/s_erfl.c @@ -97,6 +97,7 @@ */ #include <errno.h> +#include <float.h> #include <math.h> #include <math_private.h> @@ -143,9 +144,7 @@ tiny = 1e-4931L, 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) @@ -782,7 +781,16 @@ __erfl (long double x) if (ix < 0x3fc60000) /* |x|<2**-57 */ { if (ix < 0x00080000) - return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */ + { + /* Avoid spurious 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; } y = a + a * neval (z, TN1, NTN1) / deval (z, TD1, NTD1); 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 <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; } diff --git a/sysdeps/ieee754/ldbl-96/s_erfl.c b/sysdeps/ieee754/ldbl-96/s_erfl.c index 47e4b9e909..c27de812ca 100644 --- a/sysdeps/ieee754/ldbl-96/s_erfl.c +++ b/sysdeps/ieee754/ldbl-96/s_erfl.c @@ -105,6 +105,7 @@ #include <errno.h> +#include <float.h> #include <math.h> #include <math_private.h> @@ -120,8 +121,6 @@ tiny = 1e-4931L, */ /* 2/sqrt(pi) - 1 */ efx = 1.2837916709551257389615890312154517168810E-1L, - /* 8 * (2/sqrt(pi) - 1) */ - efx8 = 1.0270333367641005911692712249723613735048E0L, pp[6] = { 1.122751350964552113068262337278335028553E6L, @@ -272,7 +271,16 @@ __erfl (long double x) if (ix < 0x3fde8000) /* |x|<2**-33 */ { if (ix < 0x00080000) - return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */ + { + /* Avoid spurious 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; } z = x * x; |