https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107569

--- Comment #18 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
Ok, just so that I don't just kibbitz/review frange stuff but also try to do
something, here is my so far untested attempt at normal multiplication
fold_range  (not the x * x stuff discussed elsewhere):

--- range-op-float.cc.jj1       2022-11-09 18:00:28.612247613 +0100
+++ range-op-float.cc   2022-11-09 19:06:11.075716000 +0100
@@ -1908,6 +1908,88 @@ 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.  */
+
+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))
+       real_value_negate (&result_lb);
+      result_ub = result_lb;
+    }
+  else if (real_isinf (&op1) && real_iszero (&op2))
+    {
+      result_lb = op2;
+      if (real_isneg (&op1))
+       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) const final override
+  {
+    REAL_VALUE_TYPE cp[8];
+    // Do a cross-product.
+    frange_mult (type, cp[0], cp[4], lh_lb, rh_lb);
+    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 < 3; ++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];
+
+    // [+-0, +-0] * [+INF,+INF] (or [-INF,-INF] or swapped is a known NaN.
+    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, NULL, 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 +2024,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

For
float
foo (float x, float y)
{
  if (x < -17.0f || x > 19.5f || y < -25.0f || y > 15.5f)
    return 0.0;
  return x * y;
}

float
bar (float x, float y)
{
  if (x > -17.0f && x < 19.5f && y > -25.0f && y < 15.5f)
    return x * y;
  return 0.0;
}
the product range is
[frange] float [-4.875e+2 (-0x0.f3cp+9), 4.25e+2 (0x0.d48p+9)] +-NAN
in the first case and
[frange] float [-4.875e+2 (-0x0.f3cp+9), 4.25e+2 (0x0.d48p+9)]
in the latter (where NAN can't appear).  And while +-0.0 is in the range of
even both operands, neither range includes any infinities.

Reply via email to