Accuracy of Mathematical Functions

[Paul Zimmermann]

Hi,

I hope this is not off-topic for this list (Technical discussions relating
to FreeBSD).

We have updated our comparison:

https://members.loria.fr/PZimmermann/papers/accuracy.pdf

This new update includes for the first time the FreeBSD math library,
whose accuracy is quite good, except:

* single precision: the Bessel functions, lgammaf, cospif, sinpif, tanpif, powf
* double precision: the Bessel functions, lgammaf, tgammaf, cospi, sinpi,
tanpi, pow
* double-extended precision: erfcl, lgammal, tgammal, cospil, sinpil, tanpil,
powl

Some issues have already been fixed in the development version by Steve
Kargl (we used FreeBSD 13.2).

Best regards,
Paul


--
Posted automagically by a mail2news gateway at muc.de e.V.
Please direct questions, flames, donations, etc. to news-***@muc.de

[Dimitry Andric]

Post by Paul Zimmermann
I hope this is not off-topic for this list (Technical discussions relating
to FreeBSD).

The freebsd-numerics@ list might have been a better match, but it
receives very low traffic, and the audience on this list will be larger.

Post by Paul Zimmermann
https://members.loria.fr/PZimmermann/papers/accuracy.pdf
This new update includes for the first time the FreeBSD math library,
* single precision: the Bessel functions, lgammaf, cospif, sinpif, tanpif, powf
* double precision: the Bessel functions, lgammaf, tgammaf, cospi, sinpi,
tanpi, pow
* double-extended precision: erfcl, lgammal, tgammal, cospil, sinpil, tanpil,
powl
Some issues have already been fixed in the development version by Steve
Kargl (we used FreeBSD 13.2).

Very interesting paper! Of course we are always interested in
improvements for libm, and Steve semi-regularly posts patches in our bug
tracker. (Steve's no longer a committer, but usually these get committed
by others quickly enough.)

At the moment FreeBSD 14.0 is in beta phase, and there have been a
number of updates to libm, so it would be interesting to see if it makes
the ulp situation a bit better.

-Dimitry

[Steve Kargl]

Post by Dimitry Andric

Post by Paul Zimmermann
I hope this is not off-topic for this list (Technical discussions relating
to FreeBSD).

receives very low traffic, and the audience on this list will be larger.

I suggested Paul post here as numerics@ has had 32 post since
2021 and many of those are due to someone assigning a bug
report to numerics.

Post by Dimitry Andric

Post by Paul Zimmermann
https://members.loria.fr/PZimmermann/papers/accuracy.pdf
This new update includes for the first time the FreeBSD math library,
* single precision: the Bessel functions, lgammaf, cospif, sinpif, tanpif, powf
* double precision: the Bessel functions, lgammaf, tgammaf, cospi, sinpi,
tanpi, pow

lgammaf and tgammaf are single precision, but it is certainly possible
that there are issues. I'll take a look when I have time.

Post by Dimitry Andric

Post by Paul Zimmermann
* double-extended precision: erfcl, lgammal, tgammal, cospil, sinpil, tanpil,
powl
Some issues have already been fixed in the development version by Steve
Kargl (we used FreeBSD 13.2).

Very interesting paper! Of course we are always interested in
improvements for libm, and Steve semi-regularly posts patches in our bug
tracker. (Steve's no longer a committer, but usually these get committed
by others quickly enough.)

The issues with half-cycle trig function (i.e., cospi and friends)
should be fixed with git 2d3b0a687 on the main branch. It looks
like kib merged the patch into stable/13 with git 6fe5d4d8.

https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=272742

[Alexander Leidinger]

Post by Paul Zimmermann
https://members.loria.fr/PZimmermann/papers/accuracy.pdf
This new update includes for the first time the FreeBSD math library,

I wonder how those functions/libs you tested compare in terms of
speed...
It would allow to provide a hint to the question
"Which lib is the fastest and fulfills the needs in terms of accuracy
for the intended use-case?"

I agree that the best way to do this requires to run all libs on the
same hardware and OS, which is not feasible in your approach. What may
be feasible is to compare the relative performance of those subsets,
which you run on the same hardware.

Bye,
Alexander.

[Paul Zimmermann]

Dear Alexander,

since you ask, you will find some graphs in the last slides of my talk
about the correctly-rounded power function (x^y) at Arith23:

https://members.loria.fr/PZimmermann/talks/pow-arith23.pdf

Once CORE-MATH provides all binary64 functions, we plan to redo these graphs
with all existing libraries, including freebsd of course.

Best regards,
Paul

Date: Tue, 26 Sep 2023 15:26:16 +0200
Organization: No organization, this is a private message.
[1:text/plain Hide]

Post by Paul Zimmermann
https://members.loria.fr/PZimmermann/papers/accuracy.pdf
This new update includes for the first time the FreeBSD math library,

I wonder how those functions/libs you tested compare in terms of
speed...
It would allow to provide a hint to the question
"Which lib is the fastest and fulfills the needs in terms of accuracy
for the intended use-case?"
I agree that the best way to do this requires to run all libs on the
same hardware and OS, which is not feasible in your approach. What may
be feasible is to compare the relative performance of those subsets,
which you run on the same hardware.
Bye,
Alexander.
--
[2:application/pgp-signature Show Save:signature.asc (833B)]

--
Posted automagically by a mail2news gateway at muc.de e.V.
Please direct questions, flames, donations, etc. to news-***@muc.de

[Paul Zimmermann]

Post by Steve Kargl

Post by Paul Zimmermann
* double precision: the Bessel functions, lgammaf, tgammaf, cospi, sinpi,
tanpi, pow

lgammaf and tgammaf are single precision, but it is certainly possible
that there are issues. I'll take a look when I have time.

sorry, one should read:

* double precision: the Bessel functions, lgamma, tgamma, cospi, sinpi,
tanpi, pow

Paul


--
Posted automagically by a mail2news gateway at muc.de e.V.
Please direct questions, flames, donations, etc. to news-***@muc.de

[Steve Kargl]

Post by Paul Zimmermann
https://members.loria.fr/PZimmermann/papers/accuracy.pdf
This new update includes for the first time the FreeBSD math library,

I wonder how those functions/libs you tested compare in terms of speed...
It would allow to provide a hint to the question
"Which lib is the fastest and fulfills the needs in terms of accuracy for
the intended use-case?"
I agree that the best way to do this requires to run all libs on the same
hardware and OS, which is not feasible in your approach. What may be
feasible is to compare the relative performance of those subsets, which you
run on the same hardware.

Speed vs accuracy is always a trade-off. Consider

% ./tlibm j0 -f -x 0x1p-9 -X 2 -s 1 -N 4
Interval tested for j0f: [0.00195312,2]
4000000 calls, 0.040901 secs, 0.01023 usecs/call

% ./tlibm j0 -f -x 0x1p-9 -X 2 -s 1 -N 4
Interval tested for j0f: [0.00195312,2]
4000000 calls, 0.092471 secs, 0.02312 usecs/call

The former timing is for FreeBSD libm on an AMD FX-8350 using
only -O2 optimization. The latter is a patched libm, which uses
-O2 -funroll-loops -march=bdver2, and is twice as slow! The
difference lies in accuracy. The former gives

% ./tlibm j0 -f -x 0x1p-9 -X 2 -PD -N 4
Interval tested for j0f: [0.00195312,2]
ulp <= 0.5: 85.04% 3401545 | 85.039% 3401545
0.5 < ulp < 0.6: 5.376% 215028 | 90.414% 3616573
0.6 < ulp < 0.7: 3.107% 124266 | 93.521% 3740839
0.7 < ulp < 0.8: 2.432% 97284 | 95.953% 3838123
0.8 < ulp < 0.9: 1.740% 69612 | 97.693% 3907735
0.9 < ulp < 1.0: 0.941% 37646 | 98.635% 3945381
1.0 < ulp < 1.5: 1.108% 44312 | 99.742% 3989693
1.5 < ulp < 2.0: 0.195% 7791 | 99.937% 3997484
2.0 < ulp < 3.0: 0.062% 2491 | 99.999% 3999975
3.0 < ulp < 0.0: 0.001% 25 | 100.000% 4000000
Max ulp: 3.259556 at 1.9667229652404785e+00

while the latter has

% ./tlibm j0 -f -x 0x1p-9 -X 2 -PD -N 4
Interval tested for j0f: [0.00195312,2]
ulp <= 0.5: 86.76% 3470362 | 86.759% 3470362
0.5 < ulp < 0.6: 5.531% 221257 | 92.290% 3691619
0.6 < ulp < 0.7: 2.761% 110437 | 95.051% 3802056
0.7 < ulp < 0.8: 1.705% 68195 | 96.756% 3870251
0.8 < ulp < 0.9: 1.228% 49134 | 97.985% 3919385
0.9 < ulp < 1.0: 0.841% 33628 | 98.825% 3953013
1.0 < ulp < 1.5: 1.087% 43475 | 99.912% 3996488
1.5 < ulp < 2.0: 0.087% 3473 | 99.999% 3999961
2.0 < ulp < 3.0: 0.001% 39 | 100.000% 4000000
Max ulp: 2.157274 at 1.9673234224319458e+00

The latter is more accurate, but its underlying algorithm
uses summation-and-carry of the ascending series. This
algorithm is sensitive to compiler options, so I haven't
pushed it FreeBSD (, yet).

[Alexander Leidinger]

Post by Steve Kargl

Post by Paul Zimmermann
https://members.loria.fr/PZimmermann/papers/accuracy.pdf
This new update includes for the first time the FreeBSD math library,

I wonder how those functions/libs you tested compare in terms of speed...
It would allow to provide a hint to the question
"Which lib is the fastest and fulfills the needs in terms of accuracy for
the intended use-case?"
I agree that the best way to do this requires to run all libs on the same
hardware and OS, which is not feasible in your approach. What may be
feasible is to compare the relative performance of those subsets, which you
run on the same hardware.

Speed vs accuracy is always a trade-off. Consider

Yes.

[examples]

Post by Steve Kargl
The latter is more accurate, but its underlying algorithm
uses summation-and-carry of the ascending series. This
algorithm is sensitive to compiler options, so I haven't
pushed it FreeBSD (, yet).

A thought just crossed my mind... should we consider to provide two ABI
compatible math libs, one fast (and "acceptable accurate"... whatever
this means), and one accurate (and this one being the default to link
against)? Users then could use libmap.conf(5) to use the one according
to their needs.

Bye,
Alexander.

[Steve Kargl]

Post by Alexander Leidinger

Post by Steve Kargl

Post by Paul Zimmermann
https://members.loria.fr/PZimmermann/papers/accuracy.pdf
This new update includes for the first time the FreeBSD math library,

I wonder how those functions/libs you tested compare in terms of speed...
It would allow to provide a hint to the question
"Which lib is the fastest and fulfills the needs in terms of accuracy for
the intended use-case?"
I agree that the best way to do this requires to run all libs on the same
hardware and OS, which is not feasible in your approach. What may be
feasible is to compare the relative performance of those subsets, which you
run on the same hardware.

Speed vs accuracy is always a trade-off. Consider

Yes.
[examples]

Post by Steve Kargl
The latter is more accurate, but its underlying algorithm
uses summation-and-carry of the ascending series. This
algorithm is sensitive to compiler options, so I haven't
pushed it FreeBSD (, yet).

A thought just crossed my mind... should we consider to provide two ABI
compatible math libs, one fast (and "acceptable accurate"... whatever this
means), and one accurate (and this one being the default to link against)?
Users then could use libmap.conf(5) to use the one according to their needs.

This is certainly possible, but implementing it across all supported
architecture might be problematic. For float and double, FreeBSD
may be able to leverage CoreMath that Paul and his colleagues at
Inria are developing (https://core-math.gitlabpages.inria.fr/). I
have not tried to build and test it on FreeBSD, yet.

Unfortunately, compiler options like -ffast-math and -Ofast tend
to lead users down a potentially hazardous path. Everyone wants
their codes to run "fast". They just don't know that "fast" might
also mean wrong.

[Paul Zimmermann]

Hi Alexander,

Post by Alexander Leidinger
A thought just crossed my mind... should we consider to provide two ABI
compatible math libs, one fast (and "acceptable accurate"... whatever
this means), and one accurate (and this one being the default to link
against)? Users then could use libmap.conf(5) to use the one according
to their needs.

this is not the same idea, but the new C standard allows to have two
functions exp and cr_exp:

See https://www.open-std.org/jtc1/sc22/wg14/www/docs/n3096.pdf, page 454,
item 4:

Function names that begin with cr_ are potentially reserved identifiers and
may be added to the <math.h> header. The cr_ prefix is intended to indicate
a correctly rounded version of the function.

Having both under the same name, and having a runtime dispatcher that chooses
either the fast version or the correctly rounded version, is a nice idea!

Paul


--
Posted automagically by a mail2news gateway at muc.de e.V.
Please direct questions, flames, donations, etc. to news-***@muc.de