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.
This commit is contained in:
Szabolcs Nagy 2018-07-03 13:05:31 +01:00 committed by Corinna Vinschen
parent 6a85e1a4e5
commit 2805b07fa1
1 changed files with 4 additions and 2 deletions

View File

@ -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;