3 ms·
I had once surveyed algorithms and libraries that implement accurate and reasonably efficient mathematical operations. The status quo is a bit of mess. - For m
by lifthrasiir 7y ago
I had once surveyed algorithms and libraries that implement accurate and reasonably efficient mathematical operations. The status quo is a bit of mess.
- For many transcendental functions, 1 ulp error bound is reasonable but not necessarily reproducible, and 0.5 ulp (optimal) error bound algorithm exists but is significantly slower, to the point that Kahan [1] has once pointed out that some versions of MATLAB hard-coded log10(10^n) cases to give a false illusion of accuracy.
- The optimal algorithm is, fortunately, computationally tractable. Many functions have been thoroughly analyzed for the closest outputs to IEEE 754 binary64 values and therefore the number of bits required to accurately round them, and that number of bits is small enough (~150 bits) that all trig/log/exp/pow functions can be computed without bignums.
- For a long time CRlibm from INRIA was the only library implementing optimal algorithms, and it was licensed in GPL. A successor to CRlibm, Metalibm, was in development with no licensing statements available, but nowadays a version of Metalibm [2] is licensed in MIT so you can probably make use of them.
- If I recall correctly, CRlibm had way too big code footprint (several hundreds of KBs when compiled) that I couldn't really use that anyway. The main culprit seems to be five sets of functions for each possible rounding mode. I haven't checked yet if a new Metalibm is still big.
- There are many approximation algorithms and libraries if you don't want the full accuracy. One notable example is Sleef [3] that has multiple functions for desired accuracy (from 0.5001 ulp to 3.5 ulp). Assuming that we will eventually have a set of small and portable libms necessary, I believe a variation of Sleef's approach is desirable for future programming languages: a (compile-time) argument for accuracy, all the way from the optimal (0.5 ulp) to the crude (2^20 ulp?).
[1] https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT
[2] https://github.com/kalray/metalibm https://github.com/kalray/metalibm
[3] https://sleef.org/purec.xhtml https://sleef.org/purec.xhtml