From 9e51ae3cd0c7f65bdeba93b7f1d780cdb21fc269 Mon Sep 17 00:00:00 2001 From: Sergei Zimmerman Date: Tue, 25 Feb 2025 18:05:40 +0000 Subject: sysdeps/ieee754: Fix remainder sign of zero for FE_DOWNWARD (BZ #32711) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Single-precision remainderf() and quad-precision remainderl() implementation derived from Sun is affected by an issue when the result is +-0. IEEE754 requires that if remainder(x, y) = 0, its sign shall be that of x regardless of the rounding direction. The implementation seems to have assumed that x - x = +0 in all rounding modes, which is not the case. When rounding direction is roundTowardNegative the sign of an exact zero sum (or difference) is −0. Regression tests that triggered this erroneous behavior are added to math/libm-test-remainder.inc. Tested for cross riscv64 and powerpc. Original fix by: Bruce Evans in FreeBSD's a2ddfa5ea726c56dbf825763ad371c261b89b7c7. Reviewed-by: Adhemerval Zanella --- math/libm-test-remainder.inc | 4 ++++ sysdeps/ieee754/flt-32/e_remainderf.c | 3 +++ sysdeps/ieee754/ldbl-128/e_remainderl.c | 5 ++++- 3 files changed, 11 insertions(+), 1 deletion(-) diff --git a/math/libm-test-remainder.inc b/math/libm-test-remainder.inc index 0b846b2ddf..02e1e220db 100644 --- a/math/libm-test-remainder.inc +++ b/math/libm-test-remainder.inc @@ -173,6 +173,10 @@ static const struct test_ff_f_data remainder_test_data[] = TEST_ff_f (remainder, 2, -1, plus_zero, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED), TEST_ff_f (remainder, -2, 1, minus_zero, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED), TEST_ff_f (remainder, -2, -1, minus_zero, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED), + + TEST_ff_f (remainder, 0x1.08001cp-2, 0x1p-24, plus_zero, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED|XFAIL_ROUNDING_IBM128_LIBGCC), + TEST_ff_f (remainder, -0x1.003ffep-126, -0x1.8p-148, minus_zero, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED|XFAIL_ROUNDING_IBM128_LIBGCC), + TEST_ff_f (remainder, 0x1.2c3224p+17, 0x1p-5, plus_zero, NO_INEXACT_EXCEPTION|ERRNO_UNCHANGED|XFAIL_ROUNDING_IBM128_LIBGCC), }; static void diff --git a/sysdeps/ieee754/flt-32/e_remainderf.c b/sysdeps/ieee754/flt-32/e_remainderf.c index 81781363af..bb066b7162 100644 --- a/sysdeps/ieee754/flt-32/e_remainderf.c +++ b/sysdeps/ieee754/flt-32/e_remainderf.c @@ -56,6 +56,9 @@ __ieee754_remainderf(float x, float p) } } GET_FLOAT_WORD(hx,x); + /* Make sure x is not -0. This can occur only when x = p + and rounding direction is towards negative infinity. */ + if (hx==0x80000000) hx = 0; SET_FLOAT_WORD(x,hx^sx); return x; } diff --git a/sysdeps/ieee754/ldbl-128/e_remainderl.c b/sysdeps/ieee754/ldbl-128/e_remainderl.c index 07a15c2459..1e8605f258 100644 --- a/sysdeps/ieee754/ldbl-128/e_remainderl.c +++ b/sysdeps/ieee754/ldbl-128/e_remainderl.c @@ -64,7 +64,10 @@ __ieee754_remainderl(_Float128 x, _Float128 p) if(x>=p_half) x -= p; } } - GET_LDOUBLE_MSW64(hx,x); + GET_LDOUBLE_WORDS64(hx,lx,x); + /* Make sure x is not -0. This can occur only when x = p + and rounding direction is towards negative infinity. */ + if ((hx==0x8000000000000000ULL)&&(lx==0)) hx = 0; SET_LDOUBLE_MSW64(x,hx^sx); return x; } -- cgit v1.2.3