* 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).