diff options
| author | Joseph Myers <joseph@codesourcery.com> | 2020-02-14 14:16:25 +0000 |
|---|---|---|
| committer | Joseph Myers <joseph@codesourcery.com> | 2020-02-14 14:16:25 +0000 |
| commit | ad180676b83dc1782d407dbff57dabbaab0c1f71 (patch) | |
| tree | 4cc4cb4f8a8addccc9df54163408fa52096c48cd /sysdeps/ieee754/ldbl-96 | |
| parent | fa00db0a6eb755837ae5d413515e0da582b304f3 (diff) | |
| download | glibc-ad180676b83dc1782d407dbff57dabbaab0c1f71.tar.gz glibc-ad180676b83dc1782d407dbff57dabbaab0c1f71.tar.xz glibc-ad180676b83dc1782d407dbff57dabbaab0c1f71.zip | |
Adjust thresholds in Bessel function implementations (bug 14469).
A recent discussion in bug 14469 notes that a threshold in float Bessel function implementations, used to determine when to use a simpler implementation approach, results in substantially inaccurate results. As I discussed in <https://sourceware.org/ml/libc-alpha/2013-03/msg00345.html>, a heuristic argument suggests 2^(S+P) as the right order of magnitude for a suitable threshold, where S is the number of significand bits in the floating-point type and P is the number of significant bits in the representation of the floating-point type, and the float and ldbl-96 implementations use thresholds that are too small. Some threshold does need using, there or elsewhere in the implementation, to avoid spurious underflow and overflow for large arguments. This patch sets the thresholds in the affected implementations to more heuristically justifiable values. Results will still be inaccurate close to zeroes of the functions (thus this patch does *not* fix any of the bugs for Bessel function inaccuracy); fixing that would require a different implementation approach, likely along the lines described in <http://www.cl.cam.ac.uk/~jrh13/papers/bessel.ps.gz>. So the justification for a change such as this would be statistical rather than based on particular tests that had excessive errors and no longer do so (no doubt such tests could be found, but would probably be too fragile to add to the testsuite, as liable to give large errors again from very small implementation changes or even from compiler changes). See <https://sourceware.org/ml/libc-alpha/2020-02/msg00638.html> for such statistics of the resulting improvements for float functions. Tested (glibc testsuite) for x86_64.
Diffstat (limited to 'sysdeps/ieee754/ldbl-96')
| -rw-r--r-- | sysdeps/ieee754/ldbl-96/e_j0l.c | 4 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-96/e_j1l.c | 4 |
2 files changed, 4 insertions, 4 deletions
diff --git a/sysdeps/ieee754/ldbl-96/e_j0l.c b/sysdeps/ieee754/ldbl-96/e_j0l.c index 715f56fb0b..d1f06c78e8 100644 --- a/sysdeps/ieee754/ldbl-96/e_j0l.c +++ b/sysdeps/ieee754/ldbl-96/e_j0l.c @@ -134,7 +134,7 @@ __ieee754_j0l (long double x) * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) */ - if (__glibc_unlikely (ix > 0x4080)) /* 2^129 */ + if (__glibc_unlikely (ix > 0x408e)) /* 2^143 */ z = (invsqrtpi * cc) / sqrtl (x); else { @@ -236,7 +236,7 @@ __ieee754_y0l (long double x) else ss = z / cc; } - if (__glibc_unlikely (ix > 0x4080)) /* 1e39 */ + if (__glibc_unlikely (ix > 0x408e)) /* 2^143 */ z = (invsqrtpi * ss) / sqrtl (x); else { diff --git a/sysdeps/ieee754/ldbl-96/e_j1l.c b/sysdeps/ieee754/ldbl-96/e_j1l.c index 2c967a6e56..b8ace5afd1 100644 --- a/sysdeps/ieee754/ldbl-96/e_j1l.c +++ b/sysdeps/ieee754/ldbl-96/e_j1l.c @@ -138,7 +138,7 @@ __ieee754_j1l (long double x) * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x) * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x) */ - if (__glibc_unlikely (ix > 0x4080)) + if (__glibc_unlikely (ix > 0x408e)) z = (invsqrtpi * cc) / sqrtl (y); else { @@ -232,7 +232,7 @@ __ieee754_y1l (long double x) * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) * to compute the worse one. */ - if (__glibc_unlikely (ix > 0x4080)) + if (__glibc_unlikely (ix > 0x408e)) z = (invsqrtpi * ss) / sqrtl (x); else { |