[PATCH] range-op-float: Fix up frange_arithmetic [PR107967]
Aldy Hernandez
aldyh@redhat.com
Wed Dec 7 15:26:14 GMT 2022
On 12/7/22 09:29, Jakub Jelinek wrote:
> Hi!
>
> The addition of PLUS/MINUS/MULT/RDIV_EXPR frange handlers causes
> miscompilation of some of the libm routines, resulting in lots of
> glibc test failures. A part of them is purely PR107608 fold-overflow-1.c
> etc. issues, say when the code does
> return -0.5 / 0.0;
> and expects division by zero to be emitted, but we propagate -Inf
> and avoid the operation.
> But there are also various tests where we end up with different computed
> value from the expected ones. All those cases are like:
> is: inf inf
> should be: 1.18973149535723176502e+4932 0xf.fffffffffffffff0p+16380
> is: inf inf
> should be: 1.18973149535723176508575932662800701e+4932 0x1.ffffffffffffffffffffffffffffp+16383
> is: inf inf
> should be: 1.7976931348623157e+308 0x1.fffffffffffffp+1023
> is: inf inf
> should be: 3.40282346e+38 0x1.fffffep+127
> and the corresponding source looks like:
> static const double huge = 1.0e+300;
> double whatever (...) {
> ...
> return huge * huge;
> ...
> }
> which for rounding to nearest or +inf should and does return +inf, but
> for rounding to -inf or 0 should instead return nextafter (inf, -inf);
> The rules IEEE754 has are that operations on +-Inf operands are exact
> and produce +-Inf (except for the invalid ones that produce NaN) regardless
> of rounding mode, while overflows:
> "a) roundTiesToEven and roundTiesToAway carry all overflows to ∞ with the
> sign of the intermediate result.
> b) roundTowardZero carries all overflows to the format’s largest finite
> number with the sign of the intermediate result.
> c) roundTowardNegative carries positive overflows to the format’s largest
> finite number, and carries negative overflows to −∞.
> d) roundTowardPositive carries negative overflows to the format’s most
> negative finite number, and carries positive overflows to +∞."
>
> The behavior around overflows to -Inf or nextafter (-inf, inf) was actually
> handled correctly, we'd construct [-INF, -MAX] ranges in those cases
> because !real_less (&value, &result) in that case - value is finite
> but larger in magnitude than what the format can represent (but GCC
> internal's format can), while result is -INF in that case.
> But for the overflows to +Inf or nextafter (inf, -inf) was handled
> incorrectly, it tested real_less (&result, &value) rather than
> !real_less (&result, &value), the former test is true when already the
> rounding value -> result rounded down and in that case we shouldn't
> round again, we should round down when it didn't.
>
> So, in theory this could be fixed just by adding one ! character,
> - if ((mode_composite || (real_isneg (&inf) ? real_less (&result, &value)
> + if ((mode_composite || (real_isneg (&inf) ? !real_less (&result, &value)
> : !real_less (&value, &result)))
> but the following patch goes further. The distance between
> nextafter (inf, -inf) and inf is large (infinite) and expressions like
> 1.0e+300 * 1.0e+300 always produce +inf in round to nearest mode by far,
> so I think having low bound of nextafter (inf, -inf) in that case is
> unnecessary. But if it isn't multiplication but say addition and we are
> inexact and very close to the boundary between rounding to nearest
> maximum representable vs. rounding to nearest +inf, still using [MAX, +INF]
> etc. ranges seems safer because we don't know exactly what we lost in the
> inexact computation.
>
> The following patch implements that.
>
> Bootstrapped/regtested on x86_64-linux and i686-linux, ok for trunk?
>
> 2022-12-07 Jakub Jelinek <jakub@redhat.com>
>
> PR tree-optimization/107967
> * range-op-float.cc (frange_arithmetic): Fix a thinko - if
> inf is negative, use nextafter if !real_less (&result, &value)
> rather than if real_less (&result, &value). If result is +-INF
> while value is finite and -fno-rounding-math, don't do rounding
> if !inexact or if result is significantly above max representable
> value or below min representable value.
>
> * gcc.dg/pr107967-1.c: New test.
> * gcc.dg/pr107967-2.c: New test.
> * gcc.dg/pr107967-3.c: New test.
>
> --- gcc/range-op-float.cc.jj 2022-12-06 10:25:16.594848892 +0100
> +++ gcc/range-op-float.cc 2022-12-06 20:53:47.751295689 +0100
> @@ -287,9 +287,64 @@ frange_arithmetic (enum tree_code code,
>
> // Be extra careful if there may be discrepancies between the
> // compile and runtime results.
> - if ((mode_composite || (real_isneg (&inf) ? real_less (&result, &value)
> - : !real_less (&value, &result)))
> - && (inexact || !real_identical (&result, &value)))
> + bool round = false;
> + if (mode_composite)
> + round = true;
> + else if (real_isneg (&inf))
> + {
> + round = !real_less (&result, &value);
> + if (real_isinf (&result, false)
> + && !real_isinf (&value)
> + && !flag_rounding_math)
> + {
> + // Use just [+INF, +INF] rather than [MAX, +INF]
> + // even if value is larger than MAX and rounds to
> + // nearest to +INF. Unless INEXACT is true, in
> + // that case we need some extra buffer.
> + if (!inexact)
> + round = false;
> + else
> + {
> + REAL_VALUE_TYPE tmp = result, tmp2;
> + frange_nextafter (mode, tmp, inf);
> + // TMP is at this point the maximum representable
> + // number.
> + real_arithmetic (&tmp2, MINUS_EXPR, &value, &tmp);
> + if (!real_isneg (&tmp2)
> + && (REAL_EXP (&tmp2) - REAL_EXP (&tmp)
> + >= 2 - REAL_MODE_FORMAT (mode)->p))
> + round = false;
> + }
> + }
> + }
This chunk...
> + else
> + {
> + round = !real_less (&value, &result);
> + if (real_isinf (&result, true)
> + && !real_isinf (&value)
> + && !flag_rounding_math)
> + {
> + // Use just [-INF, -INF] rather than [-INF, +MAX]
> + // even if value is smaller than -MAX and rounds to
> + // nearest to -INF. Unless INEXACT is true, in
> + // that case we need some extra buffer.
> + if (!inexact)
> + round = false;
> + else
> + {
> + REAL_VALUE_TYPE tmp = result, tmp2;
> + frange_nextafter (mode, tmp, inf);
> + // TMP is at this point the minimum representable
> + // number.
> + real_arithmetic (&tmp2, MINUS_EXPR, &value, &tmp);
> + if (real_isneg (&tmp2)
> + && (REAL_EXP (&tmp2) - REAL_EXP (&tmp)
> + >= 2 - REAL_MODE_FORMAT (mode)->p))
> + round = false;
> + }
> + }
> + }
...is quite similar to this one. Could you abstract this?
Also:
> + if (real_isinf (&result, true)
> + && !real_isinf (&value)
> + && !flag_rounding_math)
Either abstract this out into a predicate since it's also repeated, or
comment it.
Aldy
> + if (round && (inexact || !real_identical (&result, &value)))
> {
> if (mode_composite)
> {
> --- gcc/testsuite/gcc.dg/pr107967-1.c.jj 2022-12-06 20:02:22.844086729 +0100
> +++ gcc/testsuite/gcc.dg/pr107967-1.c 2022-12-06 20:03:59.444683025 +0100
> @@ -0,0 +1,35 @@
> +/* PR tree-optimization/107967 */
> +/* { dg-do compile { target float64 } } */
> +/* { dg-options "-O2 -frounding-math -fno-trapping-math -fdump-tree-optimized" } */
> +/* { dg-add-options float64 } */
> +/* { dg-final { scan-tree-dump-not "return\[ \t]\*-?Inf;" "optimized" } } */
> +
> +_Float64
> +foo (void)
> +{
> + const _Float64 huge = 1.0e+300f64;
> + return huge * huge;
> +}
> +
> +_Float64
> +bar (void)
> +{
> + const _Float64 huge = 1.0e+300f64;
> + return huge * -huge;
> +}
> +
> +_Float64
> +baz (void)
> +{
> + const _Float64 a = 0x1.fffffffffffffp+1023f64;
> + const _Float64 b = 0x1.fffffffffffffp+970f64;
> + return a + b;
> +}
> +
> +_Float64
> +qux (void)
> +{
> + const _Float64 a = 0x1.fffffffffffffp+1023f64;
> + const _Float64 b = 0x1.fffffffffffffp+969f64;
> + return a + b;
> +}
> --- gcc/testsuite/gcc.dg/pr107967-2.c.jj 2022-12-06 20:02:29.683987331 +0100
> +++ gcc/testsuite/gcc.dg/pr107967-2.c 2022-12-06 20:03:48.685839355 +0100
> @@ -0,0 +1,35 @@
> +/* PR tree-optimization/107967 */
> +/* { dg-do compile { target float64 } } */
> +/* { dg-options "-O2 -fno-rounding-math -fno-trapping-math -fdump-tree-optimized" } */
> +/* { dg-add-options float64 } */
> +/* { dg-final { scan-tree-dump-times "return\[ \t]\*-?Inf;" 3 "optimized" } } */
> +
> +_Float64
> +foo (void)
> +{
> + const _Float64 huge = 1.0e+300f64;
> + return huge * huge;
> +}
> +
> +_Float64
> +bar (void)
> +{
> + const _Float64 huge = 1.0e+300f64;
> + return huge * -huge;
> +}
> +
> +_Float64
> +baz (void)
> +{
> + const _Float64 a = 0x1.fffffffffffffp+1023f64;
> + const _Float64 b = 0x1.fffffffffffffp+970f64;
> + return a + b;
> +}
> +
> +_Float64
> +qux (void)
> +{
> + const _Float64 a = 0x1.fffffffffffffp+1023f64;
> + const _Float64 b = 0x1.fffffffffffffp+969f64;
> + return a + b;
> +}
> --- gcc/testsuite/gcc.dg/pr107967-3.c.jj 2022-12-06 20:29:35.243370388 +0100
> +++ gcc/testsuite/gcc.dg/pr107967-3.c 2022-12-06 20:53:16.553748313 +0100
> @@ -0,0 +1,53 @@
> +/* PR tree-optimization/107967 */
> +/* { dg-do compile { target float64 } } */
> +/* { dg-options "-O2 -fno-rounding-math -fno-trapping-math -fdump-tree-optimized" } */
> +/* { dg-add-options float64 } */
> +/* { dg-final { scan-tree-dump-times "return\[ \t]\*-?Inf;" 3 "optimized" } } */
> +
> +_Float64
> +foo (_Float64 x)
> +{
> + if (x >= 1.0e+300f64)
> + ;
> + else
> + __builtin_unreachable ();
> + return x * x;
> +}
> +
> +_Float64
> +bar (_Float64 x)
> +{
> + if (x >= 1.0e+300f64)
> + ;
> + else
> + __builtin_unreachable ();
> + return x * -x;
> +}
> +
> +_Float64
> +baz (_Float64 a, _Float64 b)
> +{
> + if (a >= 0x1.fffffffffffffp+1023f64)
> + ;
> + else
> + __builtin_unreachable ();
> + if (b >= 0x1.p+972f64)
> + ;
> + else
> + __builtin_unreachable ();
> + return a + b;
> +}
> +
> +_Float64
> +qux (_Float64 a, _Float64 b)
> +{
> + if (a >= 0x1.fffffffffffffp+1023f64)
> + ;
> + else
> + __builtin_unreachable ();
> + if (b >= 0x1.fffffffffffffp+969f64)
> + ;
> + else
> + __builtin_unreachable ();
> + return a + b;
> +}
>
> Jakub
>
More information about the Gcc-patches
mailing list