* Fix the inaccuracy of j0f/y0f/j1f/y1f
@ 2021-02-19 9:02 Paul Zimmermann
2021-02-19 16:44 ` Michael Morrell
2021-02-24 18:59 ` Adhemerval Zanella
0 siblings, 2 replies; 10+ messages in thread
From: Paul Zimmermann @ 2021-02-19 9:02 UTC (permalink / raw)
To: libc-alpha
Hi,
I just sent a new version of patches to fix #14469, #14471, #14470, #14472.
Following Wilco's advice, this new version uses the "reduce_large" function
for argument reduction modulo pi/2.
If anyone is willing to perform exhaustive tests on some exotic platforms,
some code is available on
<https://members.loria.fr/PZimmermann/software/check_exhaustive-555.c>.
Warning: it takes about one hour on a 4-core computer for each function
and each rounding mode.
Paul
^ permalink raw reply [flat|nested] 10+ messages in thread
* RE: Fix the inaccuracy of j0f/y0f/j1f/y1f
2021-02-19 9:02 Fix the inaccuracy of j0f/y0f/j1f/y1f Paul Zimmermann
@ 2021-02-19 16:44 ` Michael Morrell
2021-02-19 21:15 ` Joseph Myers
2021-02-21 7:37 ` Paul Zimmermann
2021-02-24 18:59 ` Adhemerval Zanella
1 sibling, 2 replies; 10+ messages in thread
From: Michael Morrell @ 2021-02-19 16:44 UTC (permalink / raw)
To: libc-alpha
I'm curious about something. Where is it specified what the required accuracy is for a given math function? Before I discovered the libm-test-ulps file, I would have expected a C library to be required to computed the exact result for any defined function (possibly allowing some inaccuracy if -ffast-math were specified). I look through the POSIX spec, but couldn't find any place that mentioned an accuracy requirement.
Michael
-----Original Message-----
From: Libc-alpha <libc-alpha-bounces@sourceware.org> On Behalf Of Paul Zimmermann
Sent: Friday, February 19, 2021 1:03 AM
To: libc-alpha@sourceware.org
Subject: Fix the inaccuracy of j0f/y0f/j1f/y1f
Hi,
I just sent a new version of patches to fix #14469, #14471, #14470, #14472.
Following Wilco's advice, this new version uses the "reduce_large" function for argument reduction modulo pi/2.
If anyone is willing to perform exhaustive tests on some exotic platforms, some code is available on <https://members.loria.fr/PZimmermann/software/check_exhaustive-555.c>.
Warning: it takes about one hour on a 4-core computer for each function and each rounding mode.
Paul
^ permalink raw reply [flat|nested] 10+ messages in thread
* RE: Fix the inaccuracy of j0f/y0f/j1f/y1f
2021-02-19 16:44 ` Michael Morrell
@ 2021-02-19 21:15 ` Joseph Myers
2021-02-21 7:37 ` Paul Zimmermann
1 sibling, 0 replies; 10+ messages in thread
From: Joseph Myers @ 2021-02-19 21:15 UTC (permalink / raw)
To: Michael Morrell; +Cc: libc-alpha
On Fri, 19 Feb 2021, Michael Morrell wrote:
> I'm curious about something. Where is it specified what the required
> accuracy is for a given math function? Before I discovered the
If the relevant standards bind a function to an IEEE 754 operation (e.g.
sqrt and fma) then, for IEEE formats, it needs to be correctly rounded
(including correct exceptions, including being correct about whether
"inexact" is raised or not, including being correct about the
architecture-specific choice of whether tininess is detected before or
after rounding; it may or may not always be correct regarding whether the
underflow exception is signaled in the case of exact underflow, which does
not raise the underflow flag and so can only be detected with traps
enabled).
If the relevant standards do not bind a function to an IEEE 754 operation
(most functions), glibc's accuracy goals are described in the manual
("Errors in Math Functions"). The specific numerical limits on errors in
ulps that are accepted in libm-test-ulps files - such that errors above
those bounds are considered to be bugs in glibc - may be found in
libm-test-support.c:init_max_error. (Those are empirical bounds in that
they were based on errors actually observed with glibc implementations and
the glibc testsuite some time ago. It would be reasonable to reduce them
if we improve the implementations in glibc so that the largest error is
smaller. I expect a reduction of the bounds to 1 or 2 ulps for IEEE
formats could be achieved without harming performance, although with
significant effort involved to reduce error accumulation in various
functions.)
--
Joseph S. Myers
joseph@codesourcery.com
^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: Fix the inaccuracy of j0f/y0f/j1f/y1f
2021-02-19 16:44 ` Michael Morrell
2021-02-19 21:15 ` Joseph Myers
@ 2021-02-21 7:37 ` Paul Zimmermann
2021-02-23 19:07 ` Michael Morrell
1 sibling, 1 reply; 10+ messages in thread
From: Paul Zimmermann @ 2021-02-21 7:37 UTC (permalink / raw)
To: Michael Morrell; +Cc: libc-alpha
Dear Michael,
> I'm curious about something. Where is it specified what the required accuracy is for a given math function? Before I discovered the libm-test-ulps file, I would have expected a C library to be required to computed the exact result for any defined function (possibly allowing some inaccuracy if -ffast-math were specified). I look through the POSIX spec, but couldn't find any place that mentioned an accuracy requirement.
to complete Joseph's answer:
* the IEEE 754 standard (even in its latest revision from 2019) does not
require "correct rounding" for the mathematical functions (except for
sqrt). More precisely it defines "Additional mathematical operations"
which shall return correctly rounded results, but these operations are
not mandatory. This means that mathematical libraries can do whatever
they want, and in practice they do [1].
* for GNU libc, the current requirement is an error of at most 9 ulps
(units in last place) for all rounding modes. (In some cases this is
16 ulps instead.) Any error larger than that is considered as a bug.
For j0f/y0f/j1f/y1f, the maximal error is millions of ulps for those
functions in GNU libc 2.33 [1]. My patch reduces that maximal error
to 9 ulps.
* there is some discussion in the C standard to reserve names like cr_sin
for correctly rounded functions. I have a project to implement such
functions, and contribute them to any interested mathematical library.
If you are interested in joining that project, please tell me.
Best regards,
Paul
[1] https://members.loria.fr/PZimmermann/papers/accuracy.pdf
[2] http://www.open-std.org/jtc1/sc22/wg14/www/docs/n2596.pdf (7.31.8)
[3] https://members.loria.fr/PZimmermann/papers/coremath.pdf
^ permalink raw reply [flat|nested] 10+ messages in thread
* RE: Fix the inaccuracy of j0f/y0f/j1f/y1f
2021-02-21 7:37 ` Paul Zimmermann
@ 2021-02-23 19:07 ` Michael Morrell
2021-02-23 19:34 ` Joseph Myers
0 siblings, 1 reply; 10+ messages in thread
From: Michael Morrell @ 2021-02-23 19:07 UTC (permalink / raw)
To: Paul Zimmermann; +Cc: libc-alpha
Paul,
You state that the IEEE 754 standard doesn't require correct rounding, but I don't see where it says that. I see Section 9.2 ("Additional mathematical operations") and the list in Table 9.1 which seems to include most of the math functions in glibc (I notice that the Bessel functions like y0 aren't included) and these all must be correctly rounded. I agree that a language may choose not to define all of the operations in 9.1, but that doesn't seem to give it the right to implement the function inaccurately.
For functions not in Table 9.1 (e.g., erf, gamma, etc.) the IEEE 754 obviously doesn't impose any requirements but the C spec that you linked to doesn't seem to allow for any inaccuracies either. For example, under section 7.12.8.1, it says "The erf functions compute the error function of x". If there is no requirement to be correctly rounded, is a conforming implementation just allowed to always return 0 (which would be many ULPs off for most inputs)?
I also don't see a special case for "sqrt" that you indicate. Why do you think that one is special?
By the way, I'm getting errors trying to access your PDFs that you list below.
A question about the cr_xxx functions. Is the intent that the programmer would call these explicitly or perhaps the compiler would map "sin" to "cr_sin" if a flag like "-fcorrect-math" were used? The latter would be friendlier I expect.
I really appreciate your help in improving my understanding here.
Michael
-----Original Message-----
From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Sent: Saturday, February 20, 2021 11:38 PM
To: Michael Morrell <mmorrell@tachyum.com>
Cc: libc-alpha@sourceware.org
Subject: Re: Fix the inaccuracy of j0f/y0f/j1f/y1f
Dear Michael,
> I'm curious about something. Where is it specified what the required accuracy is for a given math function? Before I discovered the libm-test-ulps file, I would have expected a C library to be required to computed the exact result for any defined function (possibly allowing some inaccuracy if -ffast-math were specified). I look through the POSIX spec, but couldn't find any place that mentioned an accuracy requirement.
to complete Joseph's answer:
* the IEEE 754 standard (even in its latest revision from 2019) does not
require "correct rounding" for the mathematical functions (except for
sqrt). More precisely it defines "Additional mathematical operations"
which shall return correctly rounded results, but these operations are
not mandatory. This means that mathematical libraries can do whatever
they want, and in practice they do [1].
* for GNU libc, the current requirement is an error of at most 9 ulps
(units in last place) for all rounding modes. (In some cases this is
16 ulps instead.) Any error larger than that is considered as a bug.
For j0f/y0f/j1f/y1f, the maximal error is millions of ulps for those
functions in GNU libc 2.33 [1]. My patch reduces that maximal error
to 9 ulps.
* there is some discussion in the C standard to reserve names like cr_sin
for correctly rounded functions. I have a project to implement such
functions, and contribute them to any interested mathematical library.
If you are interested in joining that project, please tell me.
Best regards,
Paul
[1] https://members.loria.fr/PZimmermann/papers/accuracy.pdf
[2] http://www.open-std.org/jtc1/sc22/wg14/www/docs/n2596.pdf (7.31.8) [3] https://members.loria.fr/PZimmermann/papers/coremath.pdf
^ permalink raw reply [flat|nested] 10+ messages in thread
* RE: Fix the inaccuracy of j0f/y0f/j1f/y1f
2021-02-23 19:07 ` Michael Morrell
@ 2021-02-23 19:34 ` Joseph Myers
0 siblings, 0 replies; 10+ messages in thread
From: Joseph Myers @ 2021-02-23 19:34 UTC (permalink / raw)
To: Michael Morrell; +Cc: Paul Zimmermann, libc-alpha
On Tue, 23 Feb 2021, Michael Morrell wrote:
> Paul,
>
> You state that the IEEE 754 standard doesn't require correct rounding,
> but I don't see where it says that. I see Section 9.2 ("Additional
> mathematical operations") and the list in Table 9.1 which seems to
> include most of the math functions in glibc (I notice that the Bessel
> functions like y0 aren't included) and these all must be correctly
> rounded. I agree that a language may choose not to define all of the
> operations in 9.1, but that doesn't seem to give it the right to
> implement the function inaccurately.
If you look at the current C2x draft N2596, 5.2.4.2.2 paragraph 8 says
that the accuracy of the library functions is implementation-defined.
Annex F says how C operations relate to operations defined in IEEE 754;
the table in F.3 paragraph 1 includes an entry for sqrt, for example, so
that corresponds to the IEEE squareRoot operation, but doesn't include
most of the math.h functions; paragraph 20 gives the mapping for other
functions and says that correct rounding is not required.
--
Joseph S. Myers
joseph@codesourcery.com
^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: Fix the inaccuracy of j0f/y0f/j1f/y1f
2021-02-19 9:02 Fix the inaccuracy of j0f/y0f/j1f/y1f Paul Zimmermann
2021-02-19 16:44 ` Michael Morrell
@ 2021-02-24 18:59 ` Adhemerval Zanella
2021-02-25 10:27 ` Paul Zimmermann
1 sibling, 1 reply; 10+ messages in thread
From: Adhemerval Zanella @ 2021-02-24 18:59 UTC (permalink / raw)
To: Paul Zimmermann, libc-alpha
On 19/02/2021 06:02, Paul Zimmermann wrote:
> Hi,
>
> I just sent a new version of patches to fix #14469, #14471, #14470, #14472.
> Following Wilco's advice, this new version uses the "reduce_large" function
> for argument reduction modulo pi/2.
>
> If anyone is willing to perform exhaustive tests on some exotic platforms,
> some code is available on
> <https://members.loria.fr/PZimmermann/software/check_exhaustive-555.c>.
> Warning: it takes about one hour on a 4-core computer for each function
> and each rounding mode.
>
> Paul
>
Paul,
I have checked your set on different architectures (all of them after a
make regen-ulps):
aarch64-linux-gnu:
FAIL: math/test-float-j1
FAIL: math/test-float32-j1
powerpc-linux-gnu/powerpc-linux-gnu-power4:
FAIL: math/test-float-j1
FAIL: math/test-float32-j1
FAIL: math/test-ldouble-j1
FAIL: math/test-ldouble-y0
FAIL: math/test-ldouble-y1
powerpc64-linux-gnu:
FAIL: math/test-float-j1
FAIL: math/test-float32-j1
FAIL: math/test-ldouble-j1
FAIL: math/test-ldouble-y0
FAIL: math/test-ldouble-y1
s390x-linux-gnu:
FAIL: math/test-float-j1
FAIL: math/test-float32-j1
All of them are maximum ulps higher than what regen-ulps expects (below).
Interesting enough, for sparc64, sparc32, and armhf power4) the regenerated
maximum ulps are within 9 ulps.
Also, for powerpc-linux-gnu-power4, by removing the SET_RESTORE_ROUNDF
I see a slight improvement:
powerpc-linux-gnu-power4$ cat math/test-float-j1.out
testing float (without inline functions)
Failure: Test: j1 (0x6.874828p+4)
Result:
is: -5.85878515e-05 -0x1.eb7874p-15
should be: -5.85878042e-05 -0x1.eb785ap-15
difference: 4.72937245e-11 0x1.a00000p-35
ulp : 13.0000
max.ulp : 9.0000
Failure: Test: j1_towardzero (0x6.874828p+4)
Result:
is: -5.85877605e-05 -0x1.eb7842p-15
should be: -5.85878042e-05 -0x1.eb785ap-15
difference: 4.36557456e-11 0x1.800000p-35
ulp : 12.0000
max.ulp : 8.0000
Test suite completed:
156 test cases plus 152 tests for exception flags and
152 tests for errno executed.
2 errors occurred.
---
aarch64-linux-gnu$ cat math/test-float-j1.out testing float (without inline functions) Failure: Test: j1 (0x6.874828p+4) Result: is: -5.85878515e-05 -0x1.eb7874p-15 should be: -5.85878042e-05 -0x1.eb785ap-15 difference: 4.72937245e-11 0x1.a00000p-35 ulp : 13.0000 max.ulp : 9.0000 Failure: Test: j1_downward (0x6.874828p+4) Result: is: -5.85878516e-05 -0x1.eb7874p-15 should be: -5.85878079e-05 -0x1.eb785cp-15 difference: 4.36557456e-11 0x1.800000p-35 ulp : 12.0000 max.ulp : 8.0000 Failure: Test: j1_towardzero (0x6.874828p+4) Result: is: -5.85878515e-05 -0x1.eb7874p-15 should be: -5.85878042e-05 -0x1.eb785ap-15 difference: 4.72937244e-11 0x1.a00000p-35 ulp : 13.0000 max.ulp : 8.0000 Failure: Test: j1_upward (0x6.874828p+4) Result: is: -5.85878515e-05 -0x1.eb7874p-15 should be: -5.85878042e-05 -0x1.eb785ap-15 difference: 4.72937245e-11 0x1.a00000p-35 ulp : 13.0000 max.ulp : 9.0000 Test suite completed: 156 test cases plus 152 tests for exception flags and 152 tests for errno executed. 4 errors occurred.
powerpc-linux-gnu-power4$ cat math/test-float-j1.out
testing float (without inline functions)
Failure: Test: j1 (0x6.874828p+4)
Result:
is: -5.85878515e-05 -0x1.eb7874p-15
should be: -5.85878042e-05 -0x1.eb785ap-15
difference: 4.72937245e-11 0x1.a00000p-35
ulp : 13.0000
max.ulp : 9.0000
Failure: Test: j1_downward (0x6.874828p+4)
Result:
is: -5.85878516e-05 -0x1.eb7874p-15
should be: -5.85878079e-05 -0x1.eb785cp-15
difference: 4.36557456e-11 0x1.800000p-35
ulp : 12.0000
max.ulp : 8.0000
Failure: Test: j1_towardzero (0x6.874828p+4)
Result:
is: -5.85878515e-05 -0x1.eb7874p-15
should be: -5.85878042e-05 -0x1.eb785ap-15
difference: 4.72937244e-11 0x1.a00000p-35
ulp : 13.0000
max.ulp : 8.0000
Failure: Test: j1_upward (0x6.874828p+4)
Result:
is: -5.85878515e-05 -0x1.eb7874p-15
should be: -5.85878042e-05 -0x1.eb785ap-15
difference: 4.72937245e-11 0x1.a00000p-35
ulp : 13.0000
max.ulp : 9.0000
Test suite completed:
156 test cases plus 152 tests for exception flags and
152 tests for errno executed.
4 errors occurred.
s390x-linux-gnu> cat math/test-float-j1.out
testing float (without inline functions)
Failure: Test: j1 (0x6.874828p+4)
Result:
is: -5.85878515e-05 -0x1.eb7874p-15
should be: -5.85878042e-05 -0x1.eb785ap-15
difference: 4.72937245e-11 0x1.a00000p-35
ulp : 13.0000
max.ulp : 9.0000
Failure: Test: j1_downward (0x6.874828p+4)
Result:
is: -5.85878516e-05 -0x1.eb7874p-15
should be: -5.85878079e-05 -0x1.eb785cp-15
difference: 4.36557456e-11 0x1.800000p-35
ulp : 12.0000
max.ulp : 8.0000
Failure: Test: j1_towardzero (0x6.874828p+4)
Result:
is: -5.85878515e-05 -0x1.eb7874p-15
should be: -5.85878042e-05 -0x1.eb785ap-15
difference: 4.72937244e-11 0x1.a00000p-35
ulp : 13.0000
max.ulp : 8.0000
Failure: Test: j1_upward (0x6.874828p+4)
Result:
is: -5.85878515e-05 -0x1.eb7874p-15
should be: -5.85878042e-05 -0x1.eb785ap-15
difference: 4.72937245e-11 0x1.a00000p-35
ulp : 13.0000
max.ulp : 9.0000
Test suite completed:
156 test cases plus 152 tests for exception flags and
152 tests for errno executed.
4 errors occurred.
^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: Fix the inaccuracy of j0f/y0f/j1f/y1f
2021-02-24 18:59 ` Adhemerval Zanella
@ 2021-02-25 10:27 ` Paul Zimmermann
2021-02-25 11:24 ` Adhemerval Zanella
0 siblings, 1 reply; 10+ messages in thread
From: Paul Zimmermann @ 2021-02-25 10:27 UTC (permalink / raw)
To: Adhemerval Zanella; +Cc: libc-alpha
Dear Adhemerval,
> I have checked your set on different architectures [...]
thank you, I can reproduce the failures on the GCC compile farm.
So far it is just due to some ranges around zeroes that are too
narrow. I guess it is due to the fact that computations on x86_64
are done using double extended precision. Btw is there a way to
disable that when building glibc?
Paul
^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: Fix the inaccuracy of j0f/y0f/j1f/y1f
2021-02-25 10:27 ` Paul Zimmermann
@ 2021-02-25 11:24 ` Adhemerval Zanella
2021-02-26 6:40 ` Paul Zimmermann
0 siblings, 1 reply; 10+ messages in thread
From: Adhemerval Zanella @ 2021-02-25 11:24 UTC (permalink / raw)
To: Paul Zimmermann; +Cc: libc-alpha
On 25/02/2021 07:27, Paul Zimmermann wrote:
> Dear Adhemerval,
>
>> I have checked your set on different architectures [...]
>
> thank you, I can reproduce the failures on the GCC compile farm.
> So far it is just due to some ranges around zeroes that are too
> narrow. I guess it is due to the fact that computations on x86_64
> are done using double extended precision. Btw is there a way to
> disable that when building glibc?
Afaik -mfpmath=sse is the default for x86_64, so I am not sure the
extension precision on x86_64 s the culprit here.
I noted that cosf shows a better precision on x86_64 compared to
other architectures (based on glibc own ulp files), so it might
be the origin of the difference.
^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: Fix the inaccuracy of j0f/y0f/j1f/y1f
2021-02-25 11:24 ` Adhemerval Zanella
@ 2021-02-26 6:40 ` Paul Zimmermann
0 siblings, 0 replies; 10+ messages in thread
From: Paul Zimmermann @ 2021-02-26 6:40 UTC (permalink / raw)
To: Adhemerval Zanella; +Cc: libc-alpha
Dear Adhemerval,
> Afaik -mfpmath=sse is the default for x86_64, so I am not sure the
> extension precision on x86_64 s the culprit here.
ok
> I noted that cosf shows a better precision on x86_64 compared to
> other architectures (based on glibc own ulp files), so it might
> be the origin of the difference.
indeed the Bessel functions depend in sincosf.
Paul
^ permalink raw reply [flat|nested] 10+ messages in thread
end of thread, other threads:[~2021-02-26 6:40 UTC | newest]
Thread overview: 10+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2021-02-19 9:02 Fix the inaccuracy of j0f/y0f/j1f/y1f Paul Zimmermann
2021-02-19 16:44 ` Michael Morrell
2021-02-19 21:15 ` Joseph Myers
2021-02-21 7:37 ` Paul Zimmermann
2021-02-23 19:07 ` Michael Morrell
2021-02-23 19:34 ` Joseph Myers
2021-02-24 18:59 ` Adhemerval Zanella
2021-02-25 10:27 ` Paul Zimmermann
2021-02-25 11:24 ` Adhemerval Zanella
2021-02-26 6:40 ` Paul Zimmermann
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).