On Fri, Dec 25, 2015 at 12:42 AM, Michael Niedermayer <mich...@niedermayer.cc> wrote: > On Thu, Dec 24, 2015 at 06:29:09PM -0800, Ganesh Ajjanagadde wrote: >> On Thu, Dec 24, 2015 at 6:18 PM, Michael Niedermayer >> <mich...@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 <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 >> > >> > 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)) > > 377 ulp looks rather good to me if the closest representable arguments > are 1178 ulp apart
patch set posted introducing avpriv_exp10, and getting rid of exp10 detection in configure. > > [...] > -- > Michael GnuPG fingerprint: 9FF2128B147EF6730BADF133611EC787040B0FAB > > Concerning the gods, I have no means of knowing whether they exist or not > or of what sort they may be, because of the obscurity of the subject, and > the brevity of human life -- Protagoras > > _______________________________________________ > ffmpeg-devel mailing list > ffmpeg-devel@ffmpeg.org > http://ffmpeg.org/mailman/listinfo/ffmpeg-devel > _______________________________________________ ffmpeg-devel mailing list ffmpeg-devel@ffmpeg.org http://ffmpeg.org/mailman/listinfo/ffmpeg-devel