[PATCH] range-op: Cleanup floating point multiplication and division fold_range [PR107569]

[Jakub Jelinek <jakub@redhat.com>]

On Fri, Nov 11, 2022 at 09:52:53AM +0100, Jakub Jelinek via Gcc-patches wrote:
> 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.

Here is an incremental patch on top of this and division patch,
which does that.

2022-11-11  Jakub Jelinek  <jakub@redhat.com>

	PR tree-optimization/107569
	* range-op-float.cc (foperator_mult_div_base): New class.
	(foperator_mult, foperator_div): Derive from that and use
	protected static methods from it to simplify the code.

--- gcc/range-op-float.cc.jj	2022-11-11 10:13:30.879410560 +0100
+++ gcc/range-op-float.cc	2022-11-11 10:55:57.602617289 +0100
@@ -1911,7 +1911,125 @@ class foperator_minus : public range_ope
 } fop_minus;
 
 
-class foperator_mult : public range_operator_float
+class foperator_mult_div_base : public range_operator_float
+{
+protected:
+  // True if [lb, ub] is [+-0, +-0].
+  static bool zero_p (const REAL_VALUE_TYPE &lb,
+		      const REAL_VALUE_TYPE &ub)
+  {
+    return real_iszero (&lb) && real_iszero (&ub);
+  }
+
+  // True if +0 or -0 is in [lb, ub] range.
+  static bool contains_zero_p (const REAL_VALUE_TYPE &lb,
+			       const REAL_VALUE_TYPE &ub)
+  {
+    return (real_compare (LE_EXPR, &lb, &dconst0)
+	    && real_compare (GE_EXPR, &ub, &dconst0));
+  }
+
+  // True if [lb, ub] is [-INF, -INF] or [+INF, +INF].
+  static bool singleton_inf_p (const REAL_VALUE_TYPE &lb,
+			       const REAL_VALUE_TYPE &ub)
+  {
+    return real_isinf (&lb) && real_isinf (&ub, real_isneg (&lb));
+  }
+
+  // Return -1 if binary op result must have sign bit set,
+  // 1 if binary op result must have sign bit clear,
+  // 0 otherwise.
+  // Sign bit of binary op result is exclusive or of the
+  // operand's sign bits.
+  static int signbit_known_p (const REAL_VALUE_TYPE &lh_lb,
+			      const REAL_VALUE_TYPE &lh_ub,
+			      const REAL_VALUE_TYPE &rh_lb,
+			      const REAL_VALUE_TYPE &rh_ub)
+  {
+    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))
+	  return 1;
+	else
+	  return -1;
+      }
+    return 0;
+  }
+
+  // Set [lb, ub] to [-0, -0], [-0, +0] or [+0, +0] depending on
+  // signbit_known.
+  static void zero_range (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub,
+			  int signbit_known)
+  {
+    ub = lb = dconst0;
+    if (signbit_known <= 0)
+      lb = real_value_negate (&dconst0);
+    if (signbit_known < 0)
+      ub = lb;
+  }
+
+  // Set [lb, ub] to [-INF, -INF], [-INF, +INF] or [+INF, +INF] depending on
+  // signbit_known.
+  static void inf_range (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub,
+			 int signbit_known)
+  {
+    if (signbit_known > 0)
+      ub = lb = dconstinf;
+    else if (signbit_known < 0)
+      ub = lb = dconstninf;
+    else
+      {
+	lb = dconstninf;
+	ub = dconstinf;
+      }
+  }
+
+  // Set [lb, ub] to [-INF, -0], [-INF, +INF] or [+0, +INF] depending on
+  // signbit_known.
+  static void zero_to_inf_range (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub,
+				 int signbit_known)
+  {
+    if (signbit_known > 0)
+      {
+	lb = dconst0;
+	ub = dconstinf;
+      }
+    else if (signbit_known < 0)
+      {
+	lb = dconstninf;
+	ub = real_value_negate (&dconst0);
+      }
+    else
+      {
+	lb = dconstninf;
+	ub = dconstinf;
+      }
+  }
+
+  // Given CP[0] to CP[3] floating point values rounded to -INF,
+  // set LB to the smallest of them (treating -0 as smaller to +0).
+  // Given CP[4] to CP[7] floating point values rounded to +INF,
+  // set UB to the largest of them (treating -0 as smaller to +0).
+  static void find_range (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub,
+			  const REAL_VALUE_TYPE (&cp)[8])
+  {
+    lb = cp[0];
+    ub = cp[4];
+    for (int i = 1; i < 4; ++i)
+      {
+	if (real_less (&cp[i], &lb)
+	    || (real_iszero (&lb) && real_isnegzero (&cp[i])))
+	  lb = cp[i];
+	if (real_less (&ub, &cp[i + 4])
+	    || (real_isnegzero (&ub) && real_iszero (&cp[i + 4])))
+	  ub = cp[i + 4];
+      }
+  }
+};
+
+
+class foperator_mult : public foperator_mult_div_base
 {
   void rv_fold (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub, bool &maybe_nan,
 		tree type,
@@ -1934,14 +2052,8 @@ class foperator_mult : public range_oper
     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))))
+	if ((zero_p (lh_lb, lh_ub) && singleton_inf_p (rh_lb, rh_ub))
+	    || (zero_p (rh_lb, rh_ub) && singleton_inf_p (lh_lb, lh_ub)))
 	  {
 	    real_nan (&lb, "", 0, TYPE_MODE (type));
 	    ub = lb;
@@ -1951,70 +2063,28 @@ class foperator_mult : public range_oper
 
 	// 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)
+	if ((contains_zero_p (lh_lb, lh_ub)
 	     && (real_isinf (&rh_lb) || real_isinf (&rh_ub)))
-	    || (real_compare (LE_EXPR, &rh_lb, &dconst0)
-		&& real_compare (GE_EXPR, &rh_ub, &dconst0)
+	    || (contains_zero_p (rh_lb, rh_ub)
 		&& (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;
-	      }
+	    int signbit_known = signbit_known_p (lh_lb, lh_ub, rh_lb, rh_ub);
 
 	    // 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 (singleton_inf_p (lh_lb, lh_ub)
+		|| singleton_inf_p (rh_lb, rh_ub))
+	      return inf_range (lb, ub, signbit_known);
 
 	    // 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;
-	      }
+	    if (zero_p (lh_lb, lh_ub) || zero_p (rh_lb, rh_ub))
+	      return zero_range (lb, ub, signbit_known);
 
 	    // Otherwise one of the multiplicands could be
 	    // [0.0, nextafter (0.0, 1.0)] and the [DBL_MAX, INF]
@@ -2022,27 +2092,13 @@ class foperator_mult : public range_oper
 	    // 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;
+	    return zero_to_inf_range (lb, ub, signbit_known);
 	  }
       }
 
     REAL_VALUE_TYPE cp[8];
-    // Do a cross-product.
+    // Do a cross-product.  At this point none of the multiplications
+    // should produce a NAN.
     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)
@@ -2052,9 +2108,13 @@ class foperator_mult : public range_oper
 	// 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);
+	if (contains_zero_p (lh_lb, lh_ub))
+	  {
+	    if (real_isneg (&lh_lb) == real_isneg (&lh_ub))
+	      cp[1] = dconst0;
+	    else
+	      cp[1] = real_value_negate (&dconst0);
+	  }
 	else
 	  cp[1] = cp[0];
 	cp[2] = cp[0];
@@ -2071,22 +2131,12 @@ class foperator_mult : public range_oper
     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];
-
+    find_range (lb, ub, cp);
   }
 } fop_mult;
 
-class foperator_div : public range_operator_float
+
+class foperator_div : public foperator_mult_div_base
 {
   void rv_fold (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub, bool &maybe_nan,
 		tree type,
@@ -2097,14 +2147,8 @@ class foperator_div : public range_opera
 		relation_kind) const final override
   {
     // +-0.0 / +-0.0 or +-INF / +-INF is a known NAN.
-    if ((real_iszero (&lh_lb)
-	 && real_iszero (&lh_ub)
-	 && real_iszero (&rh_lb)
-	 && real_iszero (&rh_ub))
-	|| (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 ((zero_p (lh_lb, lh_ub) && zero_p (rh_lb, rh_ub))
+	|| (singleton_inf_p (lh_lb, lh_ub) || singleton_inf_p (rh_lb, rh_ub)))
       {
 	real_nan (&lb, "", 0, TYPE_MODE (type));
 	ub = lb;
@@ -2112,84 +2156,31 @@ class foperator_div : public range_opera
 	return;
       }
 
-    bool both_maybe_zero = false;
-    bool both_maybe_inf = false;
-    bool must_have_signbit_zero = false;
-    bool must_have_signbit_nonzero = false;
-
     // If +-0.0 is in both ranges, it is a maybe NAN.
-    if (real_compare (LE_EXPR, &lh_lb, &dconst0)
-	&& real_compare (GE_EXPR, &lh_ub, &dconst0)
-	&& real_compare (LE_EXPR, &rh_lb, &dconst0)
-	&& real_compare (GE_EXPR, &rh_ub, &dconst0))
-      {
-	both_maybe_zero = true;
-	maybe_nan = true;
-      }
+    if (contains_zero_p (lh_lb, lh_ub) && contains_zero_p (rh_lb, rh_ub))
+      maybe_nan = true;
     // If +-INF is in both ranges, it is a maybe NAN.
     else if ((real_isinf (&lh_lb) || real_isinf (&lh_ub))
 	     && (real_isinf (&rh_lb) || real_isinf (&rh_ub)))
-      {
-	both_maybe_inf = true;
-	maybe_nan = true;
-      }
+      maybe_nan = true;
     else
       maybe_nan = 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;
-      }
+    int signbit_known = signbit_known_p (lh_lb, lh_ub, rh_lb, rh_ub);
 
     // If dividend must be zero, the range is just +-0
     // (including if the divisor is +-INF).
     // If divisor must be +-INF, the range is just +-0
     // (including if the dividend is zero).
-    if ((real_iszero (&lh_lb) && real_iszero (&lh_ub))
-	|| real_isinf (&rh_lb, false)
-	|| real_isinf (&rh_ub, true))
-      {
-	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;
-      }
+    if (zero_p (lh_lb, lh_ub) || singleton_inf_p (rh_lb, rh_ub))
+      return zero_range (lb, ub, signbit_known);
 
     // If divisor must be zero, the range is just +-INF
     // (including if the dividend is +-INF).
     // If dividend must be +-INF, the range is just +-INF
     // (including if the dividend is zero).
-    if ((real_iszero (&rh_lb) && real_iszero (&rh_ub))
-	|| real_isinf (&lh_lb, false)
-	|| real_isinf (&lh_ub, true))
-      {
-	// If all the boundary signs are the same, [+INF, +INF].
-	if (must_have_signbit_zero)
-	  ub = lb = dconstinf;
-	// If divisor and dividend must have different signs,
-	// [-INF, -INF].
-	else if (must_have_signbit_nonzero)
-	  ub = lb = dconstninf;
-	// Otherwise -> [-INF, +INF] (-INF or +INF).
-	else
-	  {
-	    lb = dconstninf;
-	    ub = dconstinf;
-	  }
-	return;
-      }
+    if (zero_p (rh_lb, rh_ub) || singleton_inf_p (lh_lb, lh_ub))
+      return inf_range (lb, ub, signbit_known);
 
     // Otherwise if both operands may be zero, divisor could be
     // nextafter(0.0, +-1.0) and dividend +-0.0
@@ -2204,30 +2195,12 @@ class foperator_div : public range_opera
     // signs of divisor and dividend are always the same we have
     // [+0.0, +INF], if they are always different we have
     // [-INF, -0.0].  If they vary, VARYING.
-    if (both_maybe_zero || both_maybe_inf)
-      {
-	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;
-      }
+    if (maybe_nan)
+      return zero_to_inf_range (lb, ub, signbit_known);
 
     REAL_VALUE_TYPE cp[8];
     // Do a cross-division.  At this point none of the divisions should
     // produce a NAN.
-    gcc_assert (!maybe_nan);
     frange_arithmetic (RDIV_EXPR, type, cp[0], lh_lb, rh_lb, dconstninf);
     frange_arithmetic (RDIV_EXPR, type, cp[1], lh_lb, rh_ub, dconstninf);
     frange_arithmetic (RDIV_EXPR, type, cp[2], lh_ub, rh_lb, dconstninf);
@@ -2237,27 +2210,16 @@ class foperator_div : public range_opera
     frange_arithmetic (RDIV_EXPR, type, cp[6], lh_ub, rh_lb, dconstinf);
     frange_arithmetic (RDIV_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];
+    find_range (lb, ub, cp);
 
     // If divisor may be zero (but is not known to be only zero),
     // and dividend can't be zero, the range can go up to -INF or +INF
     // depending on the signs.
-    if (real_compare (LE_EXPR, &rh_lb, &dconst0)
-	&& real_compare (GE_EXPR, &rh_ub, &dconst0))
+    if (contains_zero_p (rh_lb, rh_ub))
       {
-	if (!must_have_signbit_zero)
+	if (signbit_known <= 0)
 	  real_inf (&lb, true);
-	if (!must_have_signbit_nonzero)
+	if (signbit_known >= 0)
 	  real_inf (&ub, false);
       }
   }


	Jakub