Message ID | Y5BO16zJ9vReV+Af@tucnak |
---|---|
State | New |
Headers | show |
Series | range-op-float: Fix up frange_arithmetic [PR107967] | expand |
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 >
--- 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; + } + } + } + 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; + } + } + } + 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; +}