Created attachment 243634 [details] patch to fix half-cycle trigonometric functions Paul Zimmermann, a MPFR developer, contacted me about large errors in the half-cycle trigonometric functions. I've have investigated these issues and developed the attached patch. The float, double, and ld80 (long double) changes have been tested. Caveat emptor: The ld128 changes have not been compiled. The ld128 changes have not been tested. I do not have access to a system that uses ld128 floating point. Here is an itemized list of changes: * lib/msun/src/math_private.h: . Add fast floor macros to compute the integer part of |x| for 0 <= |x| 01xp(N-1), where N is the precision of the type of x. These macros are used in the half-cycle trigonometric functions (e.g., sinpi(x)). . The FFLOOR80 macros is used with the Intel 80-bit extended double functions. This macors corrects an off-by-one error, which led to enormous error for |x| > 0x1p32. * lib/msun/src/s_cospif.c: * lib/msun/src/s_cospi.c: * lib/msun/ld80/s_cospil.c: . Update Copyright years. . Use FFLOOR*() macro to get integer part of |x|. . Correct handle the range 0x1p(N-1) <= |x| < 0x1pN. Here, one needs to determine if the integral value of |x| is even or odd to choose +1 or -1. If |x| >= 0x1pN, always return +1. * lib/msun/src/s_sinpif.c: * lib/msun/src/s_sinpi.c: * lib/msun/ld80/s_sinpil.c: . Update Copyright years. . Use FFLOOR*() macro to get integer part of |x|. * lib/msun/src/s_tanpif.c: * lib/msun/src/s_tanpi.c: * lib/msun/ld80/s_tanpil.c: . Update Copyright years. . For +-0.5, return +-inf. Previously, tanpi[fl]() returned an NaN. . Use FFLOOR*() to get integer part of |x|. Need to determine if the integer part is even or odd. This is used to set +-0 for |x| integral and +-inf for (n+1/2). . For 0x1p(N-1) <= |x| < 0x1pN need to determine if x is an even or odd integer to select +0 or -0. For |x| >= 0x1pN, it is always an even integer, select 0. . Note, tanpi[fl](x) is an odd function, so one needs to consider tanpi[fl](-|x|) = - tanpi[fl](|x|). * lib/msun/ld128/s_cospil.c: * lib/msun/ld128/s_sinpil.c: * lib/msun/ld128/s_tanpil.c: . Update Copyright years. . These routines use an FFLOOR128 macros, which likely should be replaced by a bit twiddling algorithm. . The same considerations above are applied to 0x1p112 <= |x| < 0x1p113, and |x| >= 0x1p113 cases. . Note, even and odd determination used fmodl(x,2.), which is likely slow.
All, Paul supplied me with his test framework. Exhaustive testing of the float implementations in round-to-nearest yields Checking cospi MPFR library: 4.2.0-p9 MPFR header: 4.2.0-p9 (based on 4.2.0) Checking function cospif with MPFR_RNDN libm wrong by up to 5.01e-01 ulp(s) [1] for x=0x1.112f4ep-4 cospi gives 0x1.f4cd4cp-1 mpfr_cospi gives 0x1.f4cd4ap-1 Total: errors=110478 (0.00%) errors2=0 maxerr=5.01e-01 ulp(s) Checking sinpi MPFR library: 4.2.0-p9 MPFR header: 4.2.0-p9 (based on 4.2.0) Checking function sinpif with MPFR_RNDN libm wrong by up to 7.50e-01 ulp(s) [1] for x=0x1.72681p-129 sinpi gives 0x1.22eaa8p-127 mpfr_sinpi gives 0x1.22eaa4p-127 Total: errors=4614478 (0.11%) errors2=0 maxerr=7.50e-01 ulp(s) Checking tanpi MPFR library: 4.2.0-p9 MPFR header: 4.2.0-p9 (based on 4.2.0) Checking function tanpif with MPFR_RNDN libm wrong by up to 8.00e-01 ulp(s) [1] for x=0x1.4442bep-5 tanpi gives 0x1.fffd42p-4 mpfr_tanpi gives 0x1.fffd44p-4 Total: errors=63422266 (1.48%) errors2=0 maxerr=8.00e-01 ulp(s)
The following change is required for ld128 platforms to build: diff --git a/lib/msun/ld128/s_tanpil.c b/lib/msun/ld128/s_tanpil.c index b9ec3ef50f80..52003d0bf6c5 100644 --- a/lib/msun/ld128/s_tanpil.c +++ b/lib/msun/ld128/s_tanpil.c @@ -28,6 +28,7 @@ * See ../src/s_tanpi.c for implementation details. */ +#include "fpmath.h" #include "math.h" #include "math_private.h" @@ -71,6 +72,8 @@ tanpil(long double x) { long double ax, hi, lo, odd, t, xf; uint32_t ix; + uint64_t lx, llx; + int16_t hx; ax = fabsl(x); @@ -112,6 +115,8 @@ tanpil(long double x) if (isinf(x) || isnan(x)) return (vzero / vzero); + EXTRACT_LDBL128_WORDS(hx, lx, llx, x); + /* * For 0x1p112 <= |x| < 0x1p113 need to determine if x is an even * or odd integer to set t = +0 or -0. If you are fine with it, I will push both this PR patch, and PR 272765.
I don't understand why your patch is necessary. The ld128/s_tanpil.c does not do bit twiddling like the other routines. Ah, now, I see. There is a bug in the patch. :( /* - * |x| >= 0x1p112 is always an integer, so return +-0. + * For 0x1p112 <= |x| < 0x1p113 need to determine if x is an even + * or odd integer to set t = +0 or -0. + * For |x| >= 0x1p113, it is always an even integer, so t = 0. */ - return (copysignl(0, x)); + t = fmodl(xf, 2.L) == 1 ? 0 : (copysign(0, (lx & 1) ? -1 : 1)); + return ((hx & 0x80000000) ? -t : t); The copysign() is wrong. In fact, the whole evaluation of t is wrong. This can be t = fmodl(ax, 2.L) == 0 ? 0 : copysignl(0.L, -1,L); That is, t = 0 for even |x| and t = -0 for odd |x|. Apologies for the mess up. For future testing, I still have an account on freefall.freebsd.org. Is there an aarch64 reference system that I can ssh into for testing. PS: Thanks for looking at this bug report.
> Is there an aarch64 reference system that I can ssh into for testing. Yes, ref13-aarch64.freebsd.org (there's also a 12 jail, ref12-aarch64)
(In reply to Steve Kargl from comment #3) Then, what should be the return value? Because even with adjusted calculation for t, th return still looks at hx (the sign bit, I suppose)? Would you provide the updated patch?
I have successfully logged into ref13-aarch64. Konstantin, please hold off on the commit until I have time to do testing of the ld128 routines. It will likely take a few days as I need to construct a testing infrastructure (i.e., build mpfr and gmp and move my code).
Ed, It's been awhile since I ran tests on a cluster system. The last one was flame.freebsd.org, which was a sparc64 system. I have a patch I would like to test. Testing can consume CPU cycles. Are there any restrictions on how a developer should approach testing (other than common courtesy)?
Created attachment 243812 [details] New diff with corrected ld128 support The attached patch corrects issues with the ld128 function in the original patch. It should be applied to a clean source tree. @kib, apologies for the original sloppy patch.
A commit in branch main references this bug: URL: https://cgit.FreeBSD.org/src/commit/?id=2d3b0a687b910c84606e4bc36176945ad5c60406 commit 2d3b0a687b910c84606e4bc36176945ad5c60406 Author: Steve Kargl <kargl@FreeBSD.org> AuthorDate: 2023-07-31 22:34:48 +0000 Commit: Konstantin Belousov <kib@FreeBSD.org> CommitDate: 2023-08-03 04:27:58 +0000 Fixes for bugs in sinpi/cospi/tanpi patch to fix half-cycle trigonometric functions Paul Zimmermann, a MPFR developer, contacted me about large errors in the half-cycle trigonometric functions. I've have investigated these issues and developed the attached patch. The float, double, and ld80 (long double) changes have been tested. Caveat emptor: The ld128 changes have not been compiled. The ld128 changes have not been tested. I do not have access to a system that uses ld128 floating point. Here is an itemized list of changes: * lib/msun/src/math_private.h: . Add fast floor macros to compute the integer part of |x| for 0 <= |x| 01xp(N-1), where N is the precision of the type of x. These macros are used in the half-cycle trigonometric functions (e.g., sinpi(x)). . The FFLOOR80 macros is used with the Intel 80-bit extended double functions. This macors corrects an off-by-one error, which led to enormous error for |x| > 0x1p32. * lib/msun/src/s_cospif.c: * lib/msun/src/s_cospi.c: * lib/msun/ld80/s_cospil.c: . Update Copyright years. . Use FFLOOR*() macro to get integer part of |x|. . Correct handle the range 0x1p(N-1) <= |x| < 0x1pN. Here, one needs to determine if the integral value of |x| is even or odd to choose +1 or -1. If |x| >= 0x1pN, always return +1. * lib/msun/src/s_sinpif.c: * lib/msun/src/s_sinpi.c: * lib/msun/ld80/s_sinpil.c: . Update Copyright years. . Use FFLOOR*() macro to get integer part of |x|. * lib/msun/src/s_tanpif.c: * lib/msun/src/s_tanpi.c: * lib/msun/ld80/s_tanpil.c: . Update Copyright years. . For +-0.5, return +-inf. Previously, tanpi[fl]() returned an NaN. . Use FFLOOR*() to get integer part of |x|. Need to determine if the integer part is even or odd. This is used to set +-0 for |x| integral and +-inf for (n+1/2). . For 0x1p(N-1) <= |x| < 0x1pN need to determine if x is an even or odd integer to select +0 or -0. For |x| >= 0x1pN, it is always an even integer, select 0. . Note, tanpi[fl](x) is an odd function, so one needs to consider tanpi[fl](-|x|) = - tanpi[fl](|x|). * lib/msun/ld128/s_cospil.c: * lib/msun/ld128/s_sinpil.c: * lib/msun/ld128/s_tanpil.c: . Update Copyright years. . These routines use an FFLOOR128 macros, which likely should be replaced by a bit twiddling algorithm. . The same considerations above are applied to 0x1p112 <= |x| < 0x1p113, and |x| >= 0x1p113 cases. . Note, even and odd determination used fmodl(x,2.), which is likely slow. PR: 272742 MFC after: 1 week lib/msun/ld128/s_cospil.c | 50 ++++++++++++++++++++---------------------- lib/msun/ld128/s_sinpil.c | 38 ++++++++++++-------------------- lib/msun/ld128/s_tanpil.c | 39 ++++++++++++++++----------------- lib/msun/ld80/s_cospil.c | 22 ++++++------------- lib/msun/ld80/s_sinpil.c | 16 +++----------- lib/msun/ld80/s_tanpil.c | 34 ++++++++++++----------------- lib/msun/src/math_private.h | 53 +++++++++++++++++++++++++++++++++++++++++++++ lib/msun/src/s_cospi.c | 20 +++++------------ lib/msun/src/s_cospif.c | 16 ++++++-------- lib/msun/src/s_sinpi.c | 13 ++--------- lib/msun/src/s_sinpif.c | 10 +++------ lib/msun/src/s_tanpi.c | 41 +++++++++++++++++------------------ lib/msun/src/s_tanpif.c | 24 ++++++++++---------- 13 files changed, 183 insertions(+), 193 deletions(-)
A commit in branch stable/13 references this bug: URL: https://cgit.FreeBSD.org/src/commit/?id=6fe5d4d8c3c7cf3a2a81f38300766a733029c763 commit 6fe5d4d8c3c7cf3a2a81f38300766a733029c763 Author: Steve Kargl <kargl@FreeBSD.org> AuthorDate: 2023-07-31 22:34:48 +0000 Commit: Konstantin Belousov <kib@FreeBSD.org> CommitDate: 2023-08-10 02:57:29 +0000 Fixes for bugs in sinpi/cospi/tanpi PR: 272742 (cherry picked from commit 2d3b0a687b910c84606e4bc36176945ad5c60406) lib/msun/ld128/s_cospil.c | 50 ++++++++++++++++++++---------------------- lib/msun/ld128/s_sinpil.c | 38 ++++++++++++-------------------- lib/msun/ld128/s_tanpil.c | 39 ++++++++++++++++----------------- lib/msun/ld80/s_cospil.c | 22 ++++++------------- lib/msun/ld80/s_sinpil.c | 16 +++----------- lib/msun/ld80/s_tanpil.c | 34 ++++++++++++----------------- lib/msun/src/math_private.h | 53 +++++++++++++++++++++++++++++++++++++++++++++ lib/msun/src/s_cospi.c | 20 +++++------------ lib/msun/src/s_cospif.c | 16 ++++++-------- lib/msun/src/s_sinpi.c | 13 ++--------- lib/msun/src/s_sinpif.c | 10 +++------ lib/msun/src/s_tanpi.c | 41 +++++++++++++++++------------------ lib/msun/src/s_tanpif.c | 24 ++++++++++---------- 13 files changed, 183 insertions(+), 193 deletions(-)