[gcc r13-6549] range-op-float: Fix up reverse binary operations [PR109008]

[gcc r13-6549] range-op-float: Fix up reverse binary operations [PR109008]

[Jakub Jelinek]

https://gcc.gnu.org/g:bad177e848787258070415dbe002b2c6fba1c511

commit r13-6549-gbad177e848787258070415dbe002b2c6fba1c511
Author: Jakub Jelinek <jakub@redhat.com>
Date:   Thu Mar 9 09:51:37 2023 +0100

    range-op-float: Fix up reverse binary operations [PR109008]
    
    The following testcase is reduced from miscompilation of scipy package.
    If we have say lhs = [1., 1.] - [1., 1.] and want to compute the range
    of lhs from it, we correctly determine it is [0., 0.] (if computations
    are exact, we generally don't try to round them further in
    frange_arithmetic).  In the testcase it is about a reverse operation,
    [1., 1.] = op1 + [1., 1.] and we want to compute range of op1 from that.
    Right now we just perform the inverse operation (there are some corner
    cases about NaN and infinities handling) and so arrive to range
    [0., 0.] as well, and because it is a singleton, optimize return eps;
    to return 0.  That is incorrect though, for the reverse ops we need to
    take into account also rounding, the right exact range is
    [-0x1.0p-54, 0x1.0p-53] in this case when rounding to nearest, i.e.
    all numbers which added to 1. with round to nearest still produce 1.
    
    The problem isn't solely on singleton ranges, and isn't solely on
    results around zero.  We basically need to consider also values
    where the result is up to 0.5ulp away from the lhs range boundaries
    in each direction.
    
    The following patch fixes it by extending the lhs range for the
    reverse operations by 1ulp in each direction.  The PR contains
    a pseudo-random test generator I've used to generate 300000 tests
    of + and - and then used the same test with * and / instead of + and -
    together with a hack to print the discovered ranges by the patch in
    a form that another test could then verify the range is conservatively
    correct and how far it is from a minimal range.
    
    I believe the results are good enough for now, though plan to look
    incrementally into trying to do something better on the -XXX_MAX or
    XXX_MAX boundaries (where I think frange_nextafter will use -inf or +inf)
    and also try to increase the range just by 0.5ulp rather than 1ulp
    if !flag_rounding_math.  But dunno if either of those will be doable
    and will pass the testing, so I think it is worth committing this fix
    first.
    
    2023-03-09  Jakub Jelinek  <jakub@redhat.com>
                Richard Biener  <rguenther@suse.de>
    
            PR tree-optimization/109008
            * range-op-float.cc (float_widen_lhs_range): New function.
            (foperator_plus::op1_range, foperator_minus::op1_range,
            foperator_minus::op2_range, foperator_mult::op1_range,
            foperator_div::op1_range, foperator_div::op2_range): Use it.
    
            * gcc.c-torture/execute/ieee/pr109008.c: New test.

Diff:
---
 gcc/range-op-float.cc                              | 81 +++++++++++++++-------
 .../gcc.c-torture/execute/ieee/pr109008.c          | 18 +++++
 2 files changed, 75 insertions(+), 24 deletions(-)

diff --git a/gcc/range-op-float.cc b/gcc/range-op-float.cc
index ff42b95de4f..63b7ad1a008 100644
--- a/gcc/range-op-float.cc
+++ b/gcc/range-op-float.cc
@@ -2199,6 +2199,33 @@ zero_to_inf_range (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub, int signbit_known)
     }
 }
 
+/* Extend the LHS range by 1ulp in each direction.  For op1_range
+   or op2_range of binary operations just computing the inverse
+   operation on ranges isn't sufficient.  Consider e.g.
+   [1., 1.] = op1 + [1., 1.].  op1's range is not [0., 0.], but
+   [-0x1.0p-54, 0x1.0p-53] (when not -frounding-math), any value for
+   which adding 1. to it results in 1. after rounding to nearest.
+   So, for op1_range/op2_range extend the lhs range by 1ulp in each
+   direction.  See PR109008 for more details.  */
+
+static frange
+float_widen_lhs_range (tree type, const frange &lhs)
+{
+  frange ret = lhs;
+  if (lhs.known_isnan ())
+    return ret;
+  REAL_VALUE_TYPE lb = lhs.lower_bound ();
+  REAL_VALUE_TYPE ub = lhs.upper_bound ();
+  if (real_isfinite (&lb))
+    frange_nextafter (TYPE_MODE (type), lb, dconstninf);
+  if (real_isfinite (&ub))
+    frange_nextafter (TYPE_MODE (type), ub, dconstinf);
+  ret.set (type, lb, ub);
+  ret.clear_nan ();
+  ret.union_ (lhs);
+  return ret;
+}
+
 class foperator_plus : public range_operator_float
 {
   using range_operator_float::op1_range;
@@ -2214,8 +2241,9 @@ public:
     range_op_handler minus (MINUS_EXPR, type);
     if (!minus)
       return false;
-    return float_binary_op_range_finish (minus.fold_range (r, type, lhs, op2),
-					 r, type, lhs);
+    frange wlhs = float_widen_lhs_range (type, lhs);
+    return float_binary_op_range_finish (minus.fold_range (r, type, wlhs, op2),
+					 r, type, wlhs);
   }
   virtual bool op2_range (frange &r, tree type,
 			  const frange &lhs,
@@ -2260,9 +2288,10 @@ public:
   {
     if (lhs.undefined_p ())
       return false;
-    return float_binary_op_range_finish (fop_plus.fold_range (r, type, lhs,
+    frange wlhs = float_widen_lhs_range (type, lhs);
+    return float_binary_op_range_finish (fop_plus.fold_range (r, type, wlhs,
 							      op2),
-					 r, type, lhs);
+					 r, type, wlhs);
   }
   virtual bool op2_range (frange &r, tree type,
 			  const frange &lhs,
@@ -2271,8 +2300,9 @@ public:
   {
     if (lhs.undefined_p ())
       return false;
-    return float_binary_op_range_finish (fold_range (r, type, op1, lhs),
-					 r, type, lhs);
+    frange wlhs = float_widen_lhs_range (type, lhs);
+    return float_binary_op_range_finish (fold_range (r, type, op1, wlhs),
+					 r, type, wlhs);
   }
 private:
   void rv_fold (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub, bool &maybe_nan,
@@ -2338,13 +2368,14 @@ public:
     range_op_handler rdiv (RDIV_EXPR, type);
     if (!rdiv)
       return false;
-    bool ret = rdiv.fold_range (r, type, lhs, op2);
+    frange wlhs = float_widen_lhs_range (type, lhs);
+    bool ret = rdiv.fold_range (r, type, wlhs, op2);
     if (ret == false)
       return false;
-    if (lhs.known_isnan () || op2.known_isnan () || op2.undefined_p ())
-      return float_binary_op_range_finish (ret, r, type, lhs);
-    const REAL_VALUE_TYPE &lhs_lb = lhs.lower_bound ();
-    const REAL_VALUE_TYPE &lhs_ub = lhs.upper_bound ();
+    if (wlhs.known_isnan () || op2.known_isnan () || op2.undefined_p ())
+      return float_binary_op_range_finish (ret, r, type, wlhs);
+    const REAL_VALUE_TYPE &lhs_lb = wlhs.lower_bound ();
+    const REAL_VALUE_TYPE &lhs_ub = wlhs.upper_bound ();
     const REAL_VALUE_TYPE &op2_lb = op2.lower_bound ();
     const REAL_VALUE_TYPE &op2_ub = op2.upper_bound ();
     if ((contains_zero_p (lhs_lb, lhs_ub) && contains_zero_p (op2_lb, op2_ub))
@@ -2363,7 +2394,7 @@ public:
     // or if lhs must be zero and op2 doesn't include zero, it would be
     // UNDEFINED, while rdiv.fold_range computes a zero or singleton INF
     // range.  Those are supersets of UNDEFINED, so let's keep that way.
-    return float_binary_op_range_finish (ret, r, type, lhs);
+    return float_binary_op_range_finish (ret, r, type, wlhs);
   }
   virtual bool op2_range (frange &r, tree type,
 			  const frange &lhs,
@@ -2490,13 +2521,14 @@ public:
   {
     if (lhs.undefined_p ())
       return false;
-    bool ret = fop_mult.fold_range (r, type, lhs, op2);
+    frange wlhs = float_widen_lhs_range (type, lhs);
+    bool ret = fop_mult.fold_range (r, type, wlhs, op2);
     if (!ret)
       return ret;
-    if (lhs.known_isnan () || op2.known_isnan () || op2.undefined_p ())
-      return float_binary_op_range_finish (ret, r, type, lhs);
-    const REAL_VALUE_TYPE &lhs_lb = lhs.lower_bound ();
-    const REAL_VALUE_TYPE &lhs_ub = lhs.upper_bound ();
+    if (wlhs.known_isnan () || op2.known_isnan () || op2.undefined_p ())
+      return float_binary_op_range_finish (ret, r, type, wlhs);
+    const REAL_VALUE_TYPE &lhs_lb = wlhs.lower_bound ();
+    const REAL_VALUE_TYPE &lhs_ub = wlhs.upper_bound ();
     const REAL_VALUE_TYPE &op2_lb = op2.lower_bound ();
     const REAL_VALUE_TYPE &op2_ub = op2.upper_bound ();
     if ((contains_zero_p (lhs_lb, lhs_ub)
@@ -2512,7 +2544,7 @@ public:
 	zero_to_inf_range (lb, ub, signbit_known);
 	r.set (type, lb, ub);
       }
-    return float_binary_op_range_finish (ret, r, type, lhs);
+    return float_binary_op_range_finish (ret, r, type, wlhs);
   }
   virtual bool op2_range (frange &r, tree type,
 			  const frange &lhs,
@@ -2521,13 +2553,14 @@ public:
   {
     if (lhs.undefined_p ())
       return false;
-    bool ret = fold_range (r, type, op1, lhs);
+    frange wlhs = float_widen_lhs_range (type, lhs);
+    bool ret = fold_range (r, type, op1, wlhs);
     if (!ret)
       return ret;
-    if (lhs.known_isnan () || op1.known_isnan () || op1.undefined_p ())
-      return float_binary_op_range_finish (ret, r, type, lhs, true);
-    const REAL_VALUE_TYPE &lhs_lb = lhs.lower_bound ();
-    const REAL_VALUE_TYPE &lhs_ub = lhs.upper_bound ();
+    if (wlhs.known_isnan () || op1.known_isnan () || op1.undefined_p ())
+      return float_binary_op_range_finish (ret, r, type, wlhs, true);
+    const REAL_VALUE_TYPE &lhs_lb = wlhs.lower_bound ();
+    const REAL_VALUE_TYPE &lhs_ub = wlhs.upper_bound ();
     const REAL_VALUE_TYPE &op1_lb = op1.lower_bound ();
     const REAL_VALUE_TYPE &op1_ub = op1.upper_bound ();
     if ((contains_zero_p (lhs_lb, lhs_ub) && contains_zero_p (op1_lb, op1_ub))
@@ -2542,7 +2575,7 @@ public:
 	zero_to_inf_range (lb, ub, signbit_known);
 	r.set (type, lb, ub);
       }
-    return float_binary_op_range_finish (ret, r, type, lhs, true);
+    return float_binary_op_range_finish (ret, r, type, wlhs, true);
   }
 private:
   void rv_fold (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub, bool &maybe_nan,
diff --git a/gcc/testsuite/gcc.c-torture/execute/ieee/pr109008.c b/gcc/testsuite/gcc.c-torture/execute/ieee/pr109008.c
new file mode 100644
index 00000000000..610aed4cac0
--- /dev/null
+++ b/gcc/testsuite/gcc.c-torture/execute/ieee/pr109008.c
@@ -0,0 +1,18 @@
+/* PR tree-optimization/109008 */
+
+__attribute__((noipa)) double
+foo (double eps)
+{
+  double d = 1. + eps;
+  if (d == 1.)
+    return eps;
+  return 0.0;
+}
+
+int
+main ()
+{
+  if (foo (__DBL_EPSILON__ / 8.0) == 0.0)
+    __builtin_abort ();
+  return 0;
+}