[PATCH v3] sysdeps/ieee754: Fix remainder sign of zero for FE_DOWNWARD (BZ #32711)

[Adhemerval Zanella Netto]

On 25/02/25 15:05, Sergei Zimmerman wrote:
> 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 <bde@FreeBSD.org> in FreeBSD's
> a2ddfa5ea726c56dbf825763ad371c261b89b7c7.

LGTM, thanks.

Reviewed-by: Adhemerval Zanella  <adhemerval.zanella@linaro.org>

> ---
> Changes since v2:
> - Use XFAIL_ROUNDING_IBM128_LIBGCC for the added test cases.
> 
> Changes since v1:
> - Apply identical fix for incorrect sign to ieee754/ldbl-128/e_remainderl.c
> - Change the check for -0 in ieee754/flt-32/e_remainderf.c
> - Only enable new test cases for binary{32,64,128} formats
> 
>  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 1a1d3e9937..bca21af06b 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;
>  }