diff options
Diffstat (limited to 'libquadmath/math/erfq.c')
| -rw-r--r-- | libquadmath/math/erfq.c | 33 |
1 files changed, 23 insertions, 10 deletions
diff --git a/libquadmath/math/erfq.c b/libquadmath/math/erfq.c index 8d383e9..45a8c01 100644 --- a/libquadmath/math/erfq.c +++ b/libquadmath/math/erfq.c @@ -11,9 +11,9 @@ /* Modifications and expansions for 128-bit long double are Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> - and are incorporated herein by permission of the author. The author + and are incorporated herein by permission of the author. The author reserves the right to distribute this material elsewhere under different - copying permissions. These modifications are distributed here under + copying permissions. These modifications are distributed here under the following terms: This library is free software; you can redistribute it and/or @@ -96,6 +96,7 @@ * erfc/erf(NaN) is NaN */ +#include <errno.h> #include "quadmath-imp.h" @@ -142,13 +143,10 @@ deval (__float128 x, const __float128 *p, int n) static const __float128 tiny = 1e-4931Q, - half = 0.5Q, one = 1.0Q, two = 2.0Q, /* 2/sqrt(pi) - 1 */ - efx = 1.2837916709551257389615890312154517168810E-1Q, - /* 8 * (2/sqrt(pi) - 1) */ - efx8 = 1.0270333367641005911692712249723613735048E0Q; + efx = 1.2837916709551257389615890312154517168810E-1Q; /* erf(x) = x + x R(x^2) @@ -773,6 +771,8 @@ erfq (__float128 x) if (ix >= 0x3fff0000) /* |x| >= 1.0 */ { + if (ix >= 0x40030000 && sign > 0) + return one; /* x >= 16, avoid spurious underflow from erfc. */ y = erfcq (x); return (one - y); /* return (one - erfcq (x)); */ @@ -785,7 +785,12 @@ erfq (__float128 x) if (ix < 0x3fc60000) /* |x|<2**-57 */ { if (ix < 0x00080000) - return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */ + { + /* Avoid spurious underflow. */ + __float128 ret = 0.0625 * (16.0 * x + (16.0 * efx) * x); + math_check_force_underflow (ret); + return ret; + } return x + efx * x; } y = a + a * neval (z, TN1, NTN1) / deval (z, TD1, NTD1); @@ -867,7 +872,7 @@ erfcq (__float128 x) y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19); y += C19a; break; - case 9: + default: /* i == 9. */ z = x - 1.125Q; y = C20b + z * neval (z, RNr20, NRNr20) / deval (z, RDr20, NRDr20); y += C20a; @@ -921,14 +926,22 @@ erfcq (__float128 x) z = u.value; r = expq (-z * z - 0.5625) * expq ((z - x) * (z + x) + p); if ((sign & 0x80000000) == 0) - return r / x; + { + __float128 ret = r / x; + if (ret == 0) + errno = ERANGE; + return ret; + } else return two - r / x; } else { if ((sign & 0x80000000) == 0) - return tiny * tiny; + { + errno = ERANGE; + return tiny * tiny; + } else return two - tiny; } |