On 11/10/22 14:44, Jakub Jelinek wrote:
Hi!

The following patch implements frange multiplication, including the
special case of x * x.  The callers don't tell us that it is x * x,
just that it is either z = x * x or if (x == y) z = x * y;
For irange that makes no difference, but for frange it can mean
x is -0.0 and y is 0.0 if they have the same range that includes both
signed and unsigned zeros, so we need to assume result could be -0.0.

The patch causes one regression:
+FAIL: gcc.dg/fold-overflow-1.c scan-assembler-times 2139095040 2
but that is already tracked in PR107608 and affects not just the newly
added multiplication, but addition and other floating point operations
(and doesn't seem like a ranger bug but dce or whatever else).

Bootstrapped/regtested on x86_64-linux and i686-linux, ok for trunk?

Once we have division and the reverse ops for all of these, perhaps
we can do some cleanups to share common code, but the way I have division
now partly written doesn't show up many commonalities.  Multiplication
is simple, division is a nightmare.

Thanks for tackling this. I'm happy you think multiplication is simple, cause all the cross product operators make my head spin.


2022-11-10  Jakub Jelinek  <ja...@redhat.com>

        PR tree-optimization/107569
        PR tree-optimization/107591
        * range-op.h (range_operator_float::rv_fold): Add relation_kind
        argument.
        * range-op-float.cc (range_operator_float::fold_range): Name
        last argument trio and pass trio.op1_op2 () as last argument to
        rv_fold.
        (range_operator_float::rv_fold): Add relation_kind argument.
        (foperator_plus::rv_fold, foperator_minus::rv_fold): Likewise.
        (frange_mult): New function.
        (foperator_mult): New class.
        (floating_op_table::floating_op_table): Use foperator_mult for
        MULT_EXPR.

--- gcc/range-op.h.jj   2022-11-10 00:55:09.430219763 +0100
+++ gcc/range-op.h      2022-11-10 11:30:33.594114939 +0100
@@ -123,7 +123,8 @@ public:
                        const REAL_VALUE_TYPE &lh_lb,
                        const REAL_VALUE_TYPE &lh_ub,
                        const REAL_VALUE_TYPE &rh_lb,
-                       const REAL_VALUE_TYPE &rh_ub) const;
+                       const REAL_VALUE_TYPE &rh_ub,
+                       relation_kind) const;
    // Unary operations have the range of the LHS as op2.
    virtual bool fold_range (irange &r, tree type,
                           const frange &lh,
--- gcc/range-op-float.cc.jj    2022-11-10 00:55:09.318221259 +0100
+++ gcc/range-op-float.cc       2022-11-10 11:31:29.040359082 +0100
@@ -51,7 +51,7 @@ along with GCC; see the file COPYING3.
  bool
  range_operator_float::fold_range (frange &r, tree type,
                                  const frange &op1, const frange &op2,
-                                 relation_trio) const
+                                 relation_trio trio) const
  {
    if (empty_range_varying (r, type, op1, op2))
      return true;
@@ -65,7 +65,7 @@ range_operator_float::fold_range (frange
    bool maybe_nan;
    rv_fold (lb, ub, maybe_nan, type,
           op1.lower_bound (), op1.upper_bound (),
-          op2.lower_bound (), op2.upper_bound ());
+          op2.lower_bound (), op2.upper_bound (), trio.op1_op2 ());
// Handle possible NANs by saturating to the appropriate INF if only
    // one end is a NAN.  If both ends are a NAN, just return a NAN.
@@ -103,8 +103,8 @@ range_operator_float::rv_fold (REAL_VALU
                               const REAL_VALUE_TYPE &lh_lb ATTRIBUTE_UNUSED,
                               const REAL_VALUE_TYPE &lh_ub ATTRIBUTE_UNUSED,
                               const REAL_VALUE_TYPE &rh_lb ATTRIBUTE_UNUSED,
-                              const REAL_VALUE_TYPE &rh_ub ATTRIBUTE_UNUSED)
-  const
+                              const REAL_VALUE_TYPE &rh_ub ATTRIBUTE_UNUSED,
+                              relation_kind) const
  {
    lb = dconstninf;
    ub = dconstinf;
@@ -1868,7 +1868,8 @@ class foperator_plus : public range_oper
                const REAL_VALUE_TYPE &lh_lb,
                const REAL_VALUE_TYPE &lh_ub,
                const REAL_VALUE_TYPE &rh_lb,
-               const REAL_VALUE_TYPE &rh_ub) const final override
+               const REAL_VALUE_TYPE &rh_ub,
+               relation_kind) const final override
    {
      frange_arithmetic (PLUS_EXPR, type, lb, lh_lb, rh_lb, dconstninf);
      frange_arithmetic (PLUS_EXPR, type, ub, lh_ub, rh_ub, dconstinf);
@@ -1892,7 +1893,8 @@ class foperator_minus : public range_ope
                const REAL_VALUE_TYPE &lh_lb,
                const REAL_VALUE_TYPE &lh_ub,
                const REAL_VALUE_TYPE &rh_lb,
-               const REAL_VALUE_TYPE &rh_ub) const final override
+               const REAL_VALUE_TYPE &rh_ub,
+               relation_kind) const final override
    {
      frange_arithmetic (MINUS_EXPR, type, lb, lh_lb, rh_ub, dconstninf);
      frange_arithmetic (MINUS_EXPR, type, ub, lh_ub, rh_lb, dconstinf);
@@ -1908,6 +1910,123 @@ class foperator_minus : public range_ope
    }
  } fop_minus;
+/* Wrapper around frange_arithmetics, that computes the result
+   if inexact rounded to both directions.  Also, if one of the
+   operands is +-0.0 and another +-INF, return +-0.0 rather than
+   NAN.  */

s/frange_arithmetics/frange_arithmetic/

Also, would you mind written a little blurb about why it's necessary not to compute INF*0.0 as NAN. I assume it's because you're using it for the cross product and you'll set maybe_nan separately, but it's nice to spell it out.

+
+static void
+frange_mult (tree type, REAL_VALUE_TYPE &result_lb, REAL_VALUE_TYPE &result_ub,
+            const REAL_VALUE_TYPE &op1, const REAL_VALUE_TYPE &op2)
+{
+  if (real_iszero (&op1) && real_isinf (&op2))
+    {
+      result_lb = op1;
+      if (real_isneg (&op2))
+       result_lb = real_value_negate (&result_lb);
+      result_ub = result_lb;
+    }
+  else if (real_isinf (&op1) && real_iszero (&op2))
+    {
+      result_lb = op2;
+      if (real_isneg (&op1))
+       result_lb = real_value_negate (&result_lb);
+      result_ub = result_lb;
+    }
+  else
+    {
+      frange_arithmetic (MULT_EXPR, type, result_lb, op1, op2, dconstninf);
+      frange_arithmetic (MULT_EXPR, type, result_ub, op1, op2, dconstinf);
+    }
+}
+
+class foperator_mult : public range_operator_float
+{
+  void rv_fold (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub, bool &maybe_nan,
+               tree type,
+               const REAL_VALUE_TYPE &lh_lb,
+               const REAL_VALUE_TYPE &lh_ub,
+               const REAL_VALUE_TYPE &rh_lb,
+               const REAL_VALUE_TYPE &rh_ub,
+               relation_kind kind) const final override
+  {
+    REAL_VALUE_TYPE cp[8];
+    bool is_square
+      = (kind == VREL_EQ
+        && real_equal (&lh_lb, &rh_lb)
+        && real_equal (&lh_ub, &rh_ub)
+        && real_isneg (&lh_lb) == real_isneg (&rh_lb)
+        && real_isneg (&lh_ub) == real_isneg (&rh_ub));
+    // Do a cross-product.
+    frange_mult (type, cp[0], cp[4], lh_lb, rh_lb);
+    if (is_square)
+      {
+       // For x * x we can just do max (lh_lb * lh_lb, lh_ub * lh_ub)
+       // as maximum and -0.0 as minimum if 0.0 is in the range,
+       // otherwise min (lh_lb * lh_lb, lh_ub * lh_ub).
+       // -0.0 rather than 0.0 because VREL_EQ doesn't prove that
+       // x and y are bitwise equal, just that they compare equal.
+       if (real_compare (LE_EXPR, &lh_lb, &dconst0)
+           && real_compare (GE_EXPR, &lh_ub, &dconst0))
+         cp[1] = real_value_negate (&dconst0);
+       else
+         cp[1] = cp[0];
+       cp[2] = cp[0];
+       cp[5] = cp[4];
+       cp[6] = cp[4];
+      }
+    else
+      {
+       frange_mult (type, cp[1], cp[5], lh_lb, rh_ub);
+       frange_mult (type, cp[2], cp[6], lh_ub, rh_lb);
+      }
+    frange_mult (type, cp[3], cp[7], lh_ub, rh_ub);
+    for (int i = 1; i < 4; ++i)
+      {
+       if (real_less (&cp[i], &cp[0])
+           || (real_iszero (&cp[0]) && real_isnegzero (&cp[i])))
+         std::swap (cp[i], cp[0]);
+       if (real_less (&cp[4], &cp[i + 4])
+           || (real_isnegzero (&cp[4]) && real_iszero (&cp[i + 4])))
+         std::swap (cp[i + 4], cp[4]);
+      }
+    lb = cp[0];
+    ub = cp[4];
+
+    // If both operands are the same, then we know it can be +-0.0, or +-INF,
+    // but not both at the same time, so it will never be invalid unless
+    // operand was already NAN.
+    if (is_square)
+      maybe_nan = false;
+    // [+-0, +-0] * [+INF,+INF] (or [-INF,-INF] or swapped is a known NAN.
+    else if ((real_iszero (&lh_lb)
+             && real_iszero (&lh_ub)
+             && real_isinf (&rh_lb)
+             && real_isinf (&rh_ub, real_isneg (&rh_lb)))
+            || (real_iszero (&rh_lb)
+                && real_iszero (&rh_ub)
+                && real_isinf (&lh_lb)
+                && real_isinf (&lh_ub, real_isneg (&lh_lb))))
+      {
+       real_nan (&lb, "", 0, TYPE_MODE (type));
+       ub = lb;
+       maybe_nan = true;
+      }
+    // Otherwise, if one range includes zero and the other ends with +-INF,
+    // it is a maybe NAN.
+    else if (real_compare (LE_EXPR, &lh_lb, &dconst0)
+            && real_compare (GE_EXPR, &lh_ub, &dconst0)
+            && (real_isinf (&rh_lb) || real_isinf (&rh_ub)))
+      maybe_nan = true;
+    else if (real_compare (LE_EXPR, &rh_lb, &dconst0)
+            && real_compare (GE_EXPR, &rh_ub, &dconst0)
+            && (real_isinf (&lh_lb) || real_isinf (&lh_ub)))
+      maybe_nan = true;
+    else
+      maybe_nan = false;
+  }
+} fop_mult;
+
  // Instantiate a range_op_table for floating point operations.
  static floating_op_table global_floating_table;
@@ -1942,6 +2061,7 @@ floating_op_table::floating_op_table ()
    set (NEGATE_EXPR, fop_negate);
    set (PLUS_EXPR, fop_plus);
    set (MINUS_EXPR, fop_minus);
+  set (MULT_EXPR, fop_mult);
  }
// Return a pointer to the range_operator_float instance, if there is

It'd be nice to have some testcases. For example, from what I can see, the original integer multiplication code came with some tests in gcc.dg/tree-ssa/vrp13.c (commit 9983270bec0a18). It'd be nice to have some sanity checks, especially because so many things can go wrong with floats.

I'll leave it to you to decide what tests to include.

Otherwise, LGTM.

Thanks.
Aldy

Reply via email to