[FFmpeg-devel] [PATCH 00/15] replace pow(10, x) by exp10(x) across FFmpeg

[Ganesh Ajjanagadde]

On Thu, Dec 24, 2015 at 3:55 PM, Ganesh Ajjanagadde <gajjanag at mit.edu> wrote:
[...]
> 2. accuracy - yes, I am the only one who seems to care about it enough
> to bring it up everytime. On the other hand, I have documented the
> caveat and will transfer relevant information to avpriv_exp10 if we go
> that route, so I am fine with it.

My long standing faith in GNU libm has been shattered, and I am
perfectly alright with this accuracy wise. BTW, I can reduce the error
by ~ 30% with 2 extra multiplications and an addition (a negligible
cost in front of the exp) in a very easy to understand way (no "magic"
numbers). Belongs in separate patch IMHO.
For those curious, here is the sequence:
1. GNU libm makes a huge fuss about correct rounding (even 0.5 ulp),
refusing to take in slightly less accurate, but much faster functions:
https://news.ycombinator.com/item?id=8828936, particularly
https://news.ycombinator.com/item?id=8830486. Ok, I respect that
sentiment as long as they actually live by that. Experiments with sin,
cos, and other relatively simple libm functions confirmed that their
implementations are very accurate.
2. Beginning of suspicion: while working on swr/resample (and merging
in Boost's code for bessel), I noticed GNU libm actually implements j0
and other Bessel functions (man j0). They have a nice BUGS section
detailing errors up to 2e-16 on -8 to 8.
3. Work on erf - I noticed that even here, GNU's implementation is not
correctly rounded in all cases, and Boost's is ~30% faster at similar
levels of accuracy: Boost's math function implementers seem to be
pragmatists wrt such rounding,
http://www.boost.org/doc/libs/1_48_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_erf/error_function.html,
and come clean on how/to what degree things are correct. I do a man
erf, no BUGS section, nothing telling me anything regarding its
quality. I have to dig into the source to see that the claim is 1ulp,
which seems correct from some simple testing. BTW, this increased
speed, up front discussions of accuracy, readable and clean
implementations, and licensing issues are why I pull stuff from Boost
in case some of you wondered.
4. Work on exp10 - turns out their initial implementation was an
exp(log(10)*x), which suffers from accuracy loss at large/small
numbers. Old bug report:
https://sourceware.org/bugzilla/show_bug.cgi?id=13884, and apparently
"fixed" by computing 2 exps (one being a small correction term, the
other the main term),
https://github.com/andikleen/glibc/blob/rtm-devel9/sysdeps/ieee754/dbl-64/e_exp10.c.
I assumed with all that effort and "magic" constants log10_high,
log10_low (what are they?), this would actually solve the rounding
issue: there is essentially no excuse for slowing down clients 2x
unless it actually achieves GNU libm's goal of correct rounding.
The beauty is, it does not. Illustration:
arg   : -303.137207600000010643270798027515
exp10 : 7.2910890073523505e-304, 2 ulp
exp10l: 7.2910890073523489e-304, 0 ulp
simple: 7.2910890073526541e-304, 377 ulp
corr  : 7.2910890073524274e-304, 97 ulp
real  : 7.2910890073523489e-304, 0 ulp

where correction denotes the approach with a trivial correction
factor, real an arbitrary precision result via Mathematica. Looks like
exp10l is correctly rounded here, exp10 itself is 2 ulp off.

[...]