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

Ganesh Ajjanagadde gajjanag at mit.edu
Fri Dec 25 03:29:09 CET 2015


On Thu, Dec 24, 2015 at 6:18 PM, Michael Niedermayer
<michael at niedermayer.cc> wrote:
> On Thu, Dec 24, 2015 at 06:07:17PM -0800, Ganesh Ajjanagadde wrote:
>> 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
>
> how many ulps apart are exp10(x) and exp10(x + epsilon)
> that is the double and immedeatly next representable double arguments?

More precisely I think you mean exp10(nextafter(x, INFINITY)). Here
are the answers (with incorrectly rounded exp):
next  : 7.2910890073533049e-304, 1179
prev  : 7.2910890073513962e-304, 1179

i.e
exp10(nextafter(x, INFINITY))
exp10(nextafter(x, -INFINITY))

or with the correct exp10l:
next  : 7.2910890073533033e-304, 1178
prev  : 7.2910890073513954e-304, 1178

i.e
exp10l(nextafter(x, INFINITY))
exp10l(nextafter(x, -INFINITY))

>
> [...]
>
> --
> Michael     GnuPG fingerprint: 9FF2128B147EF6730BADF133611EC787040B0FAB
>
> The educated differ from the uneducated as much as the living from the
> dead. -- Aristotle
>
> _______________________________________________
> ffmpeg-devel mailing list
> ffmpeg-devel at ffmpeg.org
> http://ffmpeg.org/mailman/listinfo/ffmpeg-devel
>


More information about the ffmpeg-devel mailing list