diff mbox series

range-op, v2: Implement floating point multiplication fold_range [PR107569]

Message ID Y24NZRBs5H9In4Cr@tucnak
State New
Headers show
Series range-op, v2: Implement floating point multiplication fold_range [PR107569] | expand

Commit Message

Jakub Jelinek Nov. 11, 2022, 8:52 a.m. UTC 
On Thu, Nov 10, 2022 at 08:20:06PM +0100, Jakub Jelinek via Gcc-patches wrote:
> This made me think about it some more and I'll need to play around with it
> some more, perhaps the right thing is similarly to what I've attached for
> division to handle special cases upfront and call frange_arithmetic only
> for the safe cases.
> E.g. one case which the posted foperator_mult handles pessimistically is
> [0.0, 10.0] * [INF, INF].  This should be just [INF, INF] +-NAN IMHO,
> because the 0.0 * INF case will result in NAN, while
> nextafter (0.0, 1.0) * INF
> will be already INF and everything larger as well.
> I could in frange_mult be very conservative and for the 0 * INF cases
> set result_lb and result_ub to [0.0, INF] range (corresponding signs
> depending on the xor of sign of ops), but that would be quite pessimistic as
> well.  If one has:
> [0.0, 0.0] * [10.0, INF], the result should be just [0.0, 0.0] +-NAN,
> because again 0.0 * INF is NAN, but 0.0 * nextafter (INF, 0.0) is already 0.0.
> 
> Note, the is_square case doesn't suffer from all of this mess, the result
> is never NAN (unless operand is NAN).

Ok, here is the patch rewritten in the foperator_div style, with special
cases handled first and then the ordinary cases without problematic cases.
I guess if/once we have a plugin testing infrastructure, we could compare
the two versions of the patch, I think this one is more precise.
And, admittedly there are many similar spots with the foperator_div case
(but also with significant differences), so perhaps if foperator_{mult,div}
inherit from some derived class from range_operator_float and that class
would define various smaller helper static? methods, like this
discussed in the PR - contains_zero_p, singleton_nan_p, zero_p,
that
+           bool must_have_signbit_zero = false;
+           bool must_have_signbit_nonzero = false;
+           if (real_isneg (&lh_lb) == real_isneg (&lh_ub)
+               && real_isneg (&rh_lb) == real_isneg (&rh_ub))
+             {
+               if (real_isneg (&lh_lb) == real_isneg (&rh_ub))
+                 must_have_signbit_zero = true;
+               else
+                 must_have_signbit_nonzero = true;
+             }
returned as -1/0/1 int, and those set result (based on the above value) to
[+INF, +INF], [-INF, -INF] or [-INF, +INF]
or
[+0, +0], [-0, -0] or [-0, +0]
or
[+0, +INF], [-INF, -0] or [-INF, +INF]
and the
+    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]);
+      }
block, it could be smaller and more readable.

Thoughts?

This has been just compile tested so far.

2022-11-11  Jakub Jelinek  <jakub@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.
	(foperator_mult): New class.
	(floating_op_table::floating_op_table): Use foperator_mult for
	MULT_EXPR.



	Jakub
diff mbox series

Patch

--- gcc/range-op.h.jj	2022-11-11 08:15:20.952520590 +0100
+++ gcc/range-op.h	2022-11-11 08:48:27.649349048 +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-11 08:15:20.933520849 +0100
+++ gcc/range-op-float.cc	2022-11-11 09:39:14.950523368 +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,182 @@  class foperator_minus : public range_ope
   }
 } fop_minus;
 
+
+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
+  {
+    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));
+
+    maybe_nan = false;
+    // x * x never produces a new NAN and we only multiply the same
+    // values, so the 0 * INF problematic cases never appear there.
+    if (!is_square)
+      {
+	// [+-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, "", 0, TYPE_MODE (type));
+	    ub = lb;
+	    maybe_nan = true;
+	    return;
+	  }
+
+	// Otherwise, if one range includes zero and the other ends with +-INF,
+	// it is a maybe NAN.
+	if ((real_compare (LE_EXPR, &lh_lb, &dconst0)
+	     && real_compare (GE_EXPR, &lh_ub, &dconst0)
+	     && (real_isinf (&rh_lb) || real_isinf (&rh_ub)))
+	    || (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;
+
+	    bool must_have_signbit_zero = false;
+	    bool must_have_signbit_nonzero = false;
+	    if (real_isneg (&lh_lb) == real_isneg (&lh_ub)
+		&& real_isneg (&rh_lb) == real_isneg (&rh_ub))
+	      {
+		if (real_isneg (&lh_lb) == real_isneg (&rh_ub))
+		  must_have_signbit_zero = true;
+		else
+		  must_have_signbit_nonzero = true;
+	      }
+
+	    // If one of the ranges that includes INF is singleton
+	    // and the other range includes zero, the resulting
+	    // range is INF and NAN, because the 0 * INF boundary
+	    // case will be NAN, but already nextafter (0, 1) * INF
+	    // is INF.
+	    if ((real_isinf (&lh_lb)
+		 && real_isinf (&lh_ub, real_isneg (&lh_lb)))
+		|| (real_isinf (&rh_lb)
+		    && real_isinf (&rh_ub, real_isneg (&rh_lb))))
+	      {
+		// If all the boundary signs are the same, [+INF, +INF].
+		if (must_have_signbit_zero)
+		  ub = lb = dconstinf;
+		// If the two multiplicands have always different sign,
+		// [-INF, -INF].
+		else if (must_have_signbit_nonzero)
+		  ub = lb = dconstninf;
+		// Otherwise -> [-INF, +INF] (-INF or +INF).
+		else
+		  {
+		    lb = dconstninf;
+		    ub = dconstinf;
+		  }
+		return;
+	      }
+
+	    // If one of the multiplicands must be zero, the resulting
+	    // range is +-0 and NAN.
+	    if ((real_iszero (&lh_lb) && real_iszero (&lh_ub))
+		|| (real_iszero (&rh_lb) && real_iszero (&rh_ub)))
+	      {
+		ub = lb = dconst0;
+		// If all the boundary signs are the same, [+0.0, +0.0].
+		if (must_have_signbit_zero)
+		  ;
+		// If divisor and dividend must have different signs,
+		// [-0.0, -0.0].
+		else if (must_have_signbit_nonzero)
+		  ub = lb = real_value_negate (&dconst0);
+		// Otherwise -> [-0.0, +0.0].
+		else
+		  lb = real_value_negate (&dconst0);
+		return;
+	      }
+
+	    // Otherwise one of the multiplicands could be
+	    // [0.0, nextafter (0.0, 1.0)] and the [DBL_MAX, INF]
+	    // or similarly with different signs.  0.0 * DBL_MAX
+	    // is still 0.0, nextafter (0.0, 1.0) * INF is still INF,
+	    // so if the signs are always the same or always different,
+	    // result is [+0.0, +INF] or [-INF, -0.0], otherwise VARYING.
+	    if (must_have_signbit_zero)
+	      {
+		lb = dconst0;
+		ub = dconstinf;
+	      }
+	    else if (must_have_signbit_nonzero)
+	      {
+		lb = dconstninf;
+		ub = real_value_negate (&dconst0);
+	      }
+	    else
+	      {
+		lb = dconstninf;
+		ub = dconstinf;
+	      }
+	    return;
+	  }
+      }
+
+    REAL_VALUE_TYPE cp[8];
+    // Do a cross-product.
+    frange_arithmetic (MULT_EXPR, type, cp[0], lh_lb, rh_lb, dconstninf);
+    frange_arithmetic (MULT_EXPR, type, cp[4], lh_lb, rh_lb, dconstinf);
+    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_arithmetic (MULT_EXPR, type, cp[1], lh_lb, rh_ub, dconstninf);
+	frange_arithmetic (MULT_EXPR, type, cp[5], lh_lb, rh_ub, dconstinf);
+	frange_arithmetic (MULT_EXPR, type, cp[2], lh_ub, rh_lb, dconstninf);
+	frange_arithmetic (MULT_EXPR, type, cp[6], lh_ub, rh_lb, dconstinf);
+      }
+    frange_arithmetic (MULT_EXPR, type, cp[3], lh_ub, rh_ub, dconstninf);
+    frange_arithmetic (MULT_EXPR, type, cp[7], lh_ub, rh_ub, dconstinf);
+
+    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];
+
+  }
+} fop_mult;
+
 // Instantiate a range_op_table for floating point operations.
 static floating_op_table global_floating_table;
 
@@ -1942,6 +2120,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