On Thu, Dec 24, 2015 at 3:55 PM, Ganesh Ajjanagadde <gajja...@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. [...] _______________________________________________ ffmpeg-devel mailing list ffmpeg-devel@ffmpeg.org http://ffmpeg.org/mailman/listinfo/ffmpeg-devel
