[gcc r13-4549] range-op-float: Fix up frange_arithmetic [PR107967]

[gcc r13-4549] range-op-float: Fix up frange_arithmetic [PR107967]

[Jakub Jelinek]

https://gcc.gnu.org/g:2f187e9893916796992b60b278e743ec865f7418

commit r13-4549-g2f187e9893916796992b60b278e743ec865f7418
Author: Jakub Jelinek <jakub@redhat.com>
Date:   Thu Dec 8 10:34:26 2022 +0100

    range-op-float: Fix up frange_arithmetic [PR107967]
    
    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.
    
    2022-12-08  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.

Diff:
---
 gcc/range-op-float.cc             | 39 +++++++++++++++++++++++++---
 gcc/testsuite/gcc.dg/pr107967-1.c | 35 ++++++++++++++++++++++++++
 gcc/testsuite/gcc.dg/pr107967-2.c | 35 ++++++++++++++++++++++++++
 gcc/testsuite/gcc.dg/pr107967-3.c | 53 +++++++++++++++++++++++++++++++++++++++
 4 files changed, 159 insertions(+), 3 deletions(-)

diff --git a/gcc/range-op-float.cc b/gcc/range-op-float.cc
index c6c1137bafd..2c6026d48fa 100644
--- a/gcc/range-op-float.cc
+++ b/gcc/range-op-float.cc
@@ -287,9 +287,42 @@ frange_arithmetic (enum tree_code code, tree type,
 
   // 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
+    {
+      bool low = real_isneg (&inf);
+      round = (low ? !real_less (&result, &value)
+		   : !real_less (&value, &result));
+      if (real_isinf (&result, !low)
+	  && !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.  Similarly 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 maximum representable
+	      // number.
+	      real_arithmetic (&tmp2, MINUS_EXPR, &value, &tmp);
+	      if (real_isneg (&tmp2) != low
+		  && (REAL_EXP (&tmp2) - REAL_EXP (&tmp)
+		      >= 2 - REAL_MODE_FORMAT (mode)->p))
+		round = false;
+	    }
+	}
+    }
+  if (round && (inexact || !real_identical (&result, &value)))
     {
       if (mode_composite)
 	{
diff --git a/gcc/testsuite/gcc.dg/pr107967-1.c b/gcc/testsuite/gcc.dg/pr107967-1.c
new file mode 100644
index 00000000000..d099231024c
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/pr107967-1.c
@@ -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;
+}
diff --git a/gcc/testsuite/gcc.dg/pr107967-2.c b/gcc/testsuite/gcc.dg/pr107967-2.c
new file mode 100644
index 00000000000..554000cbac0
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/pr107967-2.c
@@ -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;
+}
diff --git a/gcc/testsuite/gcc.dg/pr107967-3.c b/gcc/testsuite/gcc.dg/pr107967-3.c
new file mode 100644
index 00000000000..9b36a1f7853
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/pr107967-3.c
@@ -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;
+}