Accuracy of single precision functions

[ViralBShah]

@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.

[kargl]

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+00

The 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 zero z0 split 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+00

A 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 < 45 and an asymptotic expansion about the zeros for x > 45. I haven't had time to pursue this.

[1] (https://www.cl.cam.ac.uk/~jrh13/papers/bessel.html)

[zimmermann6]

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]

@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.

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]

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]

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]

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]

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]

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]

This patch fixes the issue reported by Paul Zimmerman.

--- /usr/src/lib/msun/src/e_powf.c	2021-02-21 03:29:00.956878000 -0800
+++ src/e_powf.c	2021-09-06 08:17:09.800008000 -0700
@@ -136,7 +136,7 @@
     /* |y| is huge */
 	if(iy>0x4d000000) { /* if |y| > 2**27 */
 	/* over/underflow if x is not close to one */
-	    if(ix<0x3f7ffff7) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
+	    if(ix<0x3f7ffff6) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
 	    if(ix>0x3f800007) return (hy>0)? sn*huge*huge:sn*tiny*tiny;
 	/* now |1-x| is tiny <= 2**-20, suffice to compute
 	   log(x) by x-x^2/2+x^3/3-x^4/4 */

[ViralBShah]

Fix in a090d3e

[zimmermann6]

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.