This is the mail archive of the
newlib-cvs@sourceware.org
mailing list for the newlib project.
[newlib-cygwin] Fix large ulp error in pow without fma very near 1.0
- From: Corinna Vinschen <corinna at sourceware dot org>
- To: newlib-cvs at sourceware dot org
- Date: 6 Jul 2018 08:45:00 -0000
- Subject: [newlib-cygwin] Fix large ulp error in pow without fma very near 1.0
https://sourceware.org/git/gitweb.cgi?p=newlib-cygwin.git;h=2805b07fa1f54ef001c2ce78618e8efc87155232
commit 2805b07fa1f54ef001c2ce78618e8efc87155232
Author: Szabolcs Nagy <szabolcs.nagy@arm.com>
Date: Tue Jul 3 13:05:31 2018 +0100
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.
Diff:
---
newlib/libm/common/pow.c | 6 ++++--
1 file 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;