{{ message }}
Repository navigation
Accuracy of single precision functions #212
Description
@zimmermann6 has a publication comparing the accuracy of single precision functions across several libms.
https://members.loria.fr/PZimmermann/papers/accuracy.pdf
The tests are with openlibm 0.4.1, and it would be nice to be able to try the latest release. In general musl fares better.
Activity
kargl commentedon Aug 30, 2020
While the max ULP is an important metric, two additional metrics may be more important: speed and distribution of ULP (if max ULP > 1). Consider
j0f(x). It is well-known that polynomial approximations break down near zeros of a function. Exhaustive testing in the domain [1:4], which contains 1 zero, shows% tlibm -DfEsP -x 1 -X 4 j0 Interval tested for j0f: [1,4] 1000000 calls, 0.090046 secs, 0.09005 usecs/call ulp <= 0.5: 49.102% 8237966 | 49.102% 8237966 0.5 < ulp < 0.6: 7.534% 1263940 | 56.636% 9501906 0.6 < ulp < 0.7: 6.538% 1096871 | 63.174% 10598777 0.7 < ulp < 0.8: 5.690% 954682 | 68.864% 11553459 0.8 < ulp < 0.9: 4.883% 819227 | 73.747% 12372686 0.9 < ulp < 1.0: 4.015% 673617 | 77.762% 13046303 1.0 < ulp < 1.5: 12.171% 2042032 | 89.933% 15088335 1.5 < ulp < 2.0: 5.488% 920737 | 95.422% 16009072 2.0 < ulp < 3.0: 3.391% 568873 | 98.812% 16577945 3.0 < ulp < : 1.188% 199271 | 100.000% 16777216 Max ulp: 900690.687500 at 2.40482545e+00The reference function was
j0(x). The zeros of the regular Bessel functions are simple zeros, so the polynomial approximation can take this into account. For example,j0f(x) = (z - z0) * P(x)with the zeroz0split into high and low parts with a total of 48 bits of precision. This gives% ./tlibm -DfEsP -x 1 -X 4 j0 Interval tested for j0f: [1,4] 1000000 calls, 0.045775 secs, 0.04578 usecs/call ulp <= 0.5: 57.112% 9581798 | 57.112% 9581798 0.5 < ulp < 0.6: 7.607% 1276254 | 64.719% 10858052 0.6 < ulp < 0.7: 5.920% 993139 | 70.639% 11851191 0.7 < ulp < 0.8: 4.656% 781177 | 75.295% 12632368 0.8 < ulp < 0.9: 3.957% 663922 | 79.252% 13296290 0.9 < ulp < 1.0: 3.321% 557242 | 82.573% 13853532 1.0 < ulp < 1.5: 10.080% 1691123 | 92.653% 15544655 1.5 < ulp < 2.0: 4.563% 765544 | 97.216% 16310199 2.0 < ulp < 3.0: 2.529% 424342 | 99.746% 16734541 3.0 < ulp < : 0.254% 42675 | 100.000% 16777216 Max ulp: 5.246863 at 3.98081970e+00A much more accurate approximation that is twice as fast!
Unfortunately, Bessel functions are not periodic so argument reduction cannot be used. Thus, one needs multiple polynomial approximations to cover a significant domain. J. Harrison [1] proposed an algorithm that uses multiple polynomial for
x < 45and an asymptotic expansion about the zeros forx > 45. I haven't had time to pursue this.[1] (https://www.cl.cam.ac.uk/~jrh13/papers/bessel.html)
zimmermann6 commentedon Aug 31, 2020
Dear Steve, on my side I am more interested by the max ULP. I know other people are more interested by speed. I guess the Intel Math Library implements the algorithm from John Harrison.
zimmermann6 commentedon Sep 16, 2020
I have updated my results with OpenLibm 0.7.0. The results are almost the same, except for y1f whose largest error has improved. The url is still the same. Previous versions are available from https://members.loria.fr/PZimmermann/papers/
zimmermann6 commentedon Dec 4, 2020
I have revised my note with double and quadruple precision, and bivariate functions like atan2, hypot, pow:
https://homepages.loria.fr/PZimmermann/papers/#accuracy
For double precision, the largest error I found for OpenLibm 0.7.0 is less than 4 ulps, except for the Bessel functions, and the lgamma and tgamma functions.
If something has changed in newer versions of OpenLibm, please tell me and I will update my note.
zimmermann6 commentedon Feb 5, 2021
a new version of my note is now available, including the "long double" format:
https://homepages.loria.fr/PZimmermann/papers/#accuracy
I found a few issues with Openlibm, which I reported in separate new issues here.
ViralBShah commentedon Feb 6, 2021
I wonder what the easiest way is to fix some of these. Maybe just pick the good implementations from other libm libraries and bring them in here?
zimmermann6 commentedon Aug 3, 2021
I am about to release of new version of the accuracy comparison. It will be with Openlibm 0.7.5, the last release. For single precision, the large error for powf has been fixed. It remains large errors for the Bessel functions (j0f, j1f, y0f, y1f) and lgammaf. All other functions have a maximal error less than 3.17 ulps.
zimmermann6 commentedon Sep 6, 2021
I said above that in OpenLibm 0.7.5, for single precision, the large error for powf has been fixed. But I have found a new pair of inputs giving a huge error: for x=0x1.ffffecp-1 and y=-0x1.000002p+27 OpenLibm 0.7.5 yields +Inf, whereas the correct result is 0x1.557a86p115 (rounding to nearest).
kargl commentedon Sep 6, 2021
This patch fixes the issue reported by Paul Zimmerman.
Fix powf().
Fix powf: #212 (comment)
ViralBShah commentedon Sep 7, 2021
Fix in a090d3e
Fix powf().
zimmermann6 commentedon Mar 24, 2022
within the CORE-MATH project, we now have a tool to measure the efficiency of math library functions. Here is a comparison between the CORE-MATH routines and OpenLibm: https://sympa.inria.fr/sympa/arc/core-math/2022-03/msg00011.html. The CORE-MATH rountines are under MIT license, and should be easy to integrate in OpenLibm.
zimmermann6 commentedon May 18, 2022
here are updated timings for the single precision functions, with revision d14da3d of CORE-MATH, on an Intel i5-4590:
First column is CORE-MATH's reciprocal throughput, 2nd is for GNU libc, and 3rd for OpenLibm.
Squashed 'openlibm/' changes from a844d58..ae2d916
Fix powf().
awxkee commentedon Jul 18, 2025
If someone interested here is j1f with ULP 0.5 across inputs include regions near zeros.
zimmermann6 commentedon Jul 21, 2025
great ! Did you perform an exhaustive check with all binary32 values? Is that only for rounding to nearest?
awxkee commentedon Jul 21, 2025
Yes. Rust works only in rounding to nearest, thus no need to deal with other rounding modes.
I'm using libFuzzer for ULP checking and fuzzing harness. It may not be completely exhaustive, but fuzzers are quite good at finding interesting inputs.
zimmermann6 commentedon Jul 21, 2025
I strongly recommend to perform an exhaustive test. It should not take long (a few hours on a multi-core machine).
Fix powf: JuliaMath/openlibm#212 (comment)
Fix powf: JuliaMath/openlibm#212 (comment)
Squashed 'openlibm/' changes from a844d58..9fbeafc