[Bug math/28472] pow(10, i) accuracy

["wdijkstr at arm dot com" <sourceware-bugzilla@sourceware.org>]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #19 from Wilco <wdijkstr at arm dot com> ---
(In reply to Vincent Lefèvre from comment #18)
> (In reply to Wilco from comment #17)
> > (In reply to Vincent Lefèvre from comment #16)
> > > (In reply to Wilco from comment #15)
> > > > GLIBC double precision pow is the most accurate of all libraries tested at
> > > > 0.523 ULP [1].
> > > 
> > > What you forget is that this is the accuracy *tested* on arbitrary values.
> > > The actual accuracy may be worse. And this is the case here, with an
> > > accuracy larger than 1 ulp, according to the results in Comment #6!
> > 
> > Please see the implementation - it documents the accuracy across the full
> > input ranges. The worst-case reported by random testing is slightly lower
> > due to not being able to test all input values.
> 
> OK, this is something that is new to me. So perhaps [1] should also document
> the proved error bounds.

That might be too much work to figure out for every math library (AFAIK only
GLIBC documents these error bounds in the implementations). Note the final
result vary depending on FMA support, evaluation order etc.

> > Ie. if we have an algorithm that does < 0.55ULP before rounding, we can't ever > get a 2 ULP error.
> 
> OK, but in addition to the source, such proved error bounds (assuming no
> mistakes in the analysis) should be given in the glibc manual.
> 
> Also, note that an error < 0.55 ULP before rounding does not guarantee that
> pow(2,i) will be correct for the integers i. So, perhaps more should be said.

As it happens, it does work out perfectly for pow(2, i) and exp2(i). Only exp10
and pow(10, i) have a few cases that are not:

exp10(126): 0x1.7a2ecc414a04p+418 0x1.7a2ecc414a03f7ff6ca1cb527788p+418

pow(10, 23): 0x1.52d02c7e14af7p+76 0x1.52d02c7e14af68p+76 
pow(10, 210): 0x1.8557f31326bbcp+697 0x1.8557f31326bbb7fcd59ff127a01dp+697 

These are very close to the halfway rounding point and would need an extra 16
bits internal precision to get right. pow(10, 23) is exactly on the mid-point,
so is even harder.

Note GLIBC previously had some correctly rounded math functions, but they were
a total nightmare. Whenever you got close to a halfway point, performance
dropped by a factor of a million times, even for cases that were relatively
easy.

That's a lesson learnt that we don't ever want to repeat again.

-- 
You are receiving this mail because:
You are on the CC list for the bug.