diff options
author | Szabolcs Nagy <szabolcs.nagy@arm.com> | 2018-07-03 13:05:31 +0100 |
---|---|---|
committer | Corinna Vinschen <corinna@vinschen.de> | 2018-07-06 10:29:01 +0200 |
commit | 2805b07fa1f54ef001c2ce78618e8efc87155232 (patch) | |
tree | 1d54b22c61463070bd3ee28fa9014c19962685e4 | |
parent | 6a85e1a4e5a01dec3b423eb172b3ad9155487e75 (diff) | |
download | newlib-2805b07fa1f54ef001c2ce78618e8efc87155232.zip newlib-2805b07fa1f54ef001c2ce78618e8efc87155232.tar.gz newlib-2805b07fa1f54ef001c2ce78618e8efc87155232.tar.bz2 |
Fix large ulp error in pow without fma very near 1.0
The !HAVE_FAST_FMA code path split r = z/c - 1 into r = rhi + rlo such
that when z = 1-tiny and c = 1 then rlo and rhi could have much larger
magnitude than r which later caused large rounding errors.
So do a nearest rounding instead of truncation at the split.
In newlib with default settings this was observable on some arm targets
that enable the new math code but has no fma.
-rw-r--r-- | newlib/libm/common/pow.c | 6 |
1 files changed, 4 insertions, 2 deletions
diff --git a/newlib/libm/common/pow.c b/newlib/libm/common/pow.c index 11964e3..a5504e5 100644 --- a/newlib/libm/common/pow.c +++ b/newlib/libm/common/pow.c @@ -79,11 +79,13 @@ log_inline (uint64_t ix, double_t *tail) logc = T[i].logc; logctail = T[i].logctail; - /* r = z/c - 1, arranged to be exact. */ + /* Note: 1/c is j/N or j/N/2 where j is an integer in [N,2N) and + |z/c - 1| < 1/N, so r = z/c - 1 is exactly representible. */ #if HAVE_FAST_FMA r = fma (z, invc, -1.0); #else - double_t zhi = asdouble (iz & (-1ULL << 32)); + /* Split z such that rhi, rlo and rhi*rhi are exact and |rlo| <= |r|. */ + double_t zhi = asdouble ((iz + (1ULL << 31)) & (-1ULL << 32)); double_t zlo = z - zhi; double_t rhi = zhi * invc - 1.0; double_t rlo = zlo * invc; |