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

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

[Jakub Jelinek]

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.

2022-11-10  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.
	(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.  */
+
+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

	Jakub

Re: [PATCH] range-op: Implement floating point multiplication fold_range [PR107569]

[Aldy Hernandez]

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  <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.
> 	(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

Re: [PATCH] range-op: Implement floating point multiplication fold_range [PR107569]

[Jakub Jelinek]

On Thu, Nov 10, 2022 at 03:50:47PM +0100, Aldy Hernandez wrote:
> > @@ -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.

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).

> 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.

I've tried following, but it suffers from various issues:
1) we don't handle __builtin_signbit (whatever) == 0 (or != 0) as guarantee
   that in the guarded code whatever has signbit 0 or 1
2) __builtin_isinf (x) > 0 is lowered to x > DBL_MAX, but unfortunately we don't
   infer from that [INF,INF] range, but [DBL_MAX, INF] range
3) what I wrote above, I think we don't handle [0, 2] * [INF, INF] right but
   due to 2) we can't see it

So, maybe for now a selftest will be better than a testcase, or
alternatively a plugin test which acts like a selftest.

/* { dg-do compile { target { ! { vax-*-* powerpc-*-*spe pdp11-*-* } } } } */
/* { dg-options "-O2 -fno-trapping-math -fno-signaling-nans -fsigned-zeros -fno-tree-fre -fno-tree-dominator-opts -fno-thread-jumps -fdump-tree-optimized" } */
/* { dg-add-options ieee } */

void
foo (double x, double y)
{
  const double inf = __builtin_inf ();
  const double minf = -inf;
  if (__builtin_isnan (x) || __builtin_isnan (y))
    return;
#define TEST(n, xl, xu, yl, yu, rl, ru, nan) \
  if ((__builtin_isinf (xl) > 0				\
       ? x > 0.0 && __builtin_isinf (x)			\
       : __builtin_isinf (xu) < 0			\
       ? x < 0.0 && __builtin_isinf (x)			\
       : x >= xl && x <= xu				\
	 && (xl != 0.0					\
	     || __builtin_signbit (xl)			\
	     || !__builtin_signbit (x))			\
	 && (xu != 0.0					\
	     || !__builtin_signbit (xu)			\
	     || __builtin_signbit (x)))			\
      && (__builtin_isinf (yl) > 0			\
	  ? y > 0.0 && __builtin_isinf (y)		\
	  : __builtin_isinf (yu) < 0			\
	  ? y < 0.0 && __builtin_isinf (y)		\
	  : y >= yl && y <= yu				\
	    && (yl != 0.0				\
		|| __builtin_signbit (yl)		\
		|| !__builtin_signbit (y))		\
	    && (yu != 0.0				\
		|| !__builtin_signbit (yu)		\
		|| __builtin_signbit (y))))		\
    {							\
      double r##n = x * y;				\
      if (nan == 2)					\
	{						\
	  if (!__builtin_isnan (r##n))			\
	    __builtin_abort ();				\
	}						\
      else if (nan == 1)				\
	{						\
	  if (!__builtin_isnan (r##n))			\
	    {						\
	      if (r##n < rl || r##n > ru)		\
		__builtin_abort ();			\
	    }						\
	}						\
      else						\
	{						\
	  if (__builtin_isnan (r##n))			\
	    __builtin_abort ();				\
	  if (r##n < rl || r##n > ru)			\
	    __builtin_abort ();				\
	}						\
    }
#define TEST2(n, xl, xu, rl, ru) \
  if (__builtin_isinf (xl) > 0				\
      ? x > 0.0 && __builtin_isinf (x)			\
      : __builtin_isinf (xu) < 0			\
      ? x < 0.0 && __builtin_isinf (x)			\
      : x >= xl && x <= xu				\
	&& (xl != 0.0					\
	    || __builtin_signbit (xl)			\
	    || !__builtin_signbit (x))			\
	&& (xu != 0.0					\
	    || !__builtin_signbit (xu)			\
	    || __builtin_signbit (x)))			\
    {							\
      double s##n = x * x;				\
      if (__builtin_isnan (s##n))			\
	__builtin_abort ();				\
      if (s##n < rl || s##n > ru)			\
	__builtin_abort ();				\
    }
  TEST (1,	2.0, 4.0,	6.0, 8.0,	12.0, 32.0, 0);
  TEST (2,	-2.0, 3.0,	-7.0, 4.0,	-21.0, 14.0, 0);
  TEST (3,	-9.0, 5.0,	8.0, 10.0,	-90.0, 50.0, 0);
  TEST (4,	-0.0, 0.0,	16.0, 32.0,	-0.0, 0.0, 0);
  TEST (5,	-0.0, 0.0,	16.0, 32.0,	-0.0, 0.0, 0);
  TEST (6,	0.0, 0.0,	16.0, 32.0,	0.0, 0.0, 0);
  TEST (7,	0.0, 0.0,	16.0, 32.0,	0.0, 0.0, 0);
  TEST (8,	-0.0, -0.0,	16.0, 32.0,	-0.0, -0.0, 0);
  TEST (9,	-0.0, -0.0,	16.0, 32.0,	-0.0, -0.0, 0);
  TEST (10,	minf, inf,	minf, inf,	minf, inf, 1);
  TEST (11,	-0.0, 0.0,	128.0, inf,	-0.0, 0.0, 1);
  TEST (12,	-0.0, 0.0,	inf, inf,	0.0, 0.0, 2);
  TEST (13,	minf, minf,	0.0, 0.0,	0.0, 0.0, 2);
  TEST (14,	0.0, 2.0,	inf, inf,	inf, inf, 1);
  TEST (15,	0.0, 2.0,	minf, minf,	minf, minf, 1);
  TEST (16,	inf, inf,	-0.0, 2.0,	minf, inf, 1);
  TEST (17,	minf, minf,	-0.0, 3.0,	minf, inf, 1);

  TEST2 (1,	2.0, 4.0,			4.0, 16.0);
  TEST2 (2,	-0.0, 0.0,			-0.0, 0.0);
  TEST2 (3,	0.0, inf,			0.0, inf);
  TEST2 (4,	inf, inf,			inf, inf);
}


	Jakub

Re: [PATCH] range-op: Implement floating point multiplication fold_range [PR107569]

[Jakub Jelinek]

On Thu, Nov 10, 2022 at 08:20:06PM +0100, Jakub Jelinek via Gcc-patches wrote:
> So, maybe for now a selftest will be better than a testcase, or
> alternatively a plugin test which acts like a selftest.

For a testsuite/g*.dg/plugin/ range-op-plugin.c test, would be
nice to write a short plugin that does:
1) register 2 new attributes, say gcc_range and gcc_expected_range,
   parse their arguments
2) registers some new pass (dunno where it would be best, say before evrp
   or so), which would:
   - for function parameters with gcc_range attributes set global? range
     on default def of that parameter
   - if function itself has a gcc_expected_range attribute, propagate
     ranges from arguments to the function result and compare against
     gcc_expected_range
   - if function itself has gcc_range attribute and one of the arguments
     gcc_expected_range itself, try to propagate range backwards from
     result to the argument
Then we could say write tests like:
__attribute__((gcc_expected_range (12.0, 32.0, 0))) double
foo (double x __attribute__((gcc_range (2.0, 4.0, 0))), double y __attribute__((gcc_range (6.0, 8.0, 0))))
{
  return x * y;
}
with for floating point types (parameter type or function result type)
the arguments of the attribute being (constant folded)
low bound, high bound, integer about NAN (say 0 meaning clear_nan,
bit 0 meaning +NAN, bit 1 -NAN, bit 2 meaning known NAN (then
the 2 bounds would be ignored)).
Eventually we could do something similar for integral types, pointer types
etc.
I think this would be far more useful compared to writing selftests for it,
and compared to the testcase I've posted would be easier to request or check
exact range rather than a rough range.  And we could easily
test not just very simple binary ops (in both directions), but builtin calls
etc.

Thoughts on this?

I can try to write a plugin that registers the attributes, parses them
and registers a pass, but would appreciate your or Andrew's help in filling
the actual pass.

	Jakub

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

[Jakub Jelinek]

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.

--- 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


	Jakub

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

[Jakub Jelinek]

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

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

[Aldy Hernandez]

On 11/11/22 11:01, Jakub Jelinek wrote:
> 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;
> +      }
> +  }

The above functions look like they could be useful outside of the 
mult/div implementation.  Perhaps put them in file scope, instead 
limiting it to foperator_mult_div_base?

 > +  static void zero_to_inf_range (REAL_VALUE_TYPE &lb, 
REAL_VALUE_TYPE &ub,
 > +				 int signbit_known)
 > +  {
 > +    if (signbit_known > 0)

The rest of frange uses bool for a sign.  Also, real_iszero, real_isinf, 
real_inf, etc all use bool sign.  Can you use a bool, or is there a 
reason for the int?

> +
> +  // 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);
>         }
>     }

BTW, looks a lot more readable.

Thanks.
Aldy

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

[Jakub Jelinek]

On Fri, Nov 11, 2022 at 11:12:01AM +0100, Aldy Hernandez wrote:
> > --- 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;
> > +      }
> > +  }
> 
> The above functions look like they could be useful outside of the mult/div
> implementation.  Perhaps put them in file scope, instead limiting it to
> foperator_mult_div_base?

Well, I didn't want to export them to everything and most of the file
works on franges, not on REAL_VALUE_TYPE pairs.  But sure, if there
are other uses, it can be moved elsewhere.

> > +  static void zero_to_inf_range (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE
> &ub,
> > +				 int signbit_known)
> > +  {
> > +    if (signbit_known > 0)
> 
> The rest of frange uses bool for a sign.  Also, real_iszero, real_isinf,
> real_inf, etc all use bool sign.  Can you use a bool, or is there a reason
> for the int?

I need a tristate.  signbit is known and clear (this happens when
all the 4 bounds have the same sign), signbit is known and set
(this happens when one operand has signbit clear and the other signbit
set, or vice versa), or the state of the resulting signbit is unknown
(at least one operand has some values in the range with clear and others
with set signbit).

	Jakub

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

[Aldy Hernandez]

On 11/11/22 11:01, Jakub Jelinek wrote:
> 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];
> +      }
> +  }
> +};

It seems find_range() is only applicable to mult/div.  If so, perhaps 
that one should go in foperator_mult_div_base.  Unless you think it'll 
be useful for implementing other operators.

> +
> +
> +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);
>         }
>     }

I have no objections to any of your mult/div patches.  You're obviously 
more qualified to determine what's right here.  Feel free to commit when 
you're ready.

Thanks for your work on this.
Aldy

Re: [PATCH] range-op: Implement floating point multiplication fold_range [PR107569]

[Aldy Hernandez]

On 11/10/22 20:20, Jakub Jelinek wrote:
> On Thu, Nov 10, 2022 at 03:50:47PM +0100, Aldy Hernandez wrote:
>>> @@ -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.
> 
> 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).
> 
>> 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.
> 
> I've tried following, but it suffers from various issues:
> 1) we don't handle __builtin_signbit (whatever) == 0 (or != 0) as guarantee
>     that in the guarded code whatever has signbit 0 or 1

We have a range-op entry for __builtin_signbit in cfn_signbit.  Is this 
a shortcoming of this code, or something else?

> 2) __builtin_isinf (x) > 0 is lowered to x > DBL_MAX, but unfortunately we don't
>     infer from that [INF,INF] range, but [DBL_MAX, INF] range
> 3) what I wrote above, I think we don't handle [0, 2] * [INF, INF] right but
>     due to 2) we can't see it

Doesn't this boil down to a representation issue?  I wonder if we should 
bite the bullet and tweak build_gt() and build_lt() to represent open 
ranges.  In theory it should be one more/less ULP, while adjusting for 
HONOR_INFINITIES.

If the signbit issue were resolved and we could represent > and < 
properly, would that allow us to write proper testcases without having 
to writing a plug-in (which I assume is a lot harder)?

Aldy

> 
> So, maybe for now a selftest will be better than a testcase, or
> alternatively a plugin test which acts like a selftest.
> 
> /* { dg-do compile { target { ! { vax-*-* powerpc-*-*spe pdp11-*-* } } } } */
> /* { dg-options "-O2 -fno-trapping-math -fno-signaling-nans -fsigned-zeros -fno-tree-fre -fno-tree-dominator-opts -fno-thread-jumps -fdump-tree-optimized" } */
> /* { dg-add-options ieee } */
> 
> void
> foo (double x, double y)
> {
>    const double inf = __builtin_inf ();
>    const double minf = -inf;
>    if (__builtin_isnan (x) || __builtin_isnan (y))
>      return;
> #define TEST(n, xl, xu, yl, yu, rl, ru, nan) \
>    if ((__builtin_isinf (xl) > 0				\
>         ? x > 0.0 && __builtin_isinf (x)			\
>         : __builtin_isinf (xu) < 0			\
>         ? x < 0.0 && __builtin_isinf (x)			\
>         : x >= xl && x <= xu				\
> 	 && (xl != 0.0					\
> 	     || __builtin_signbit (xl)			\
> 	     || !__builtin_signbit (x))			\
> 	 && (xu != 0.0					\
> 	     || !__builtin_signbit (xu)			\
> 	     || __builtin_signbit (x)))			\
>        && (__builtin_isinf (yl) > 0			\
> 	  ? y > 0.0 && __builtin_isinf (y)		\
> 	  : __builtin_isinf (yu) < 0			\
> 	  ? y < 0.0 && __builtin_isinf (y)		\
> 	  : y >= yl && y <= yu				\
> 	    && (yl != 0.0				\
> 		|| __builtin_signbit (yl)		\
> 		|| !__builtin_signbit (y))		\
> 	    && (yu != 0.0				\
> 		|| !__builtin_signbit (yu)		\
> 		|| __builtin_signbit (y))))		\
>      {							\
>        double r##n = x * y;				\
>        if (nan == 2)					\
> 	{						\
> 	  if (!__builtin_isnan (r##n))			\
> 	    __builtin_abort ();				\
> 	}						\
>        else if (nan == 1)				\
> 	{						\
> 	  if (!__builtin_isnan (r##n))			\
> 	    {						\
> 	      if (r##n < rl || r##n > ru)		\
> 		__builtin_abort ();			\
> 	    }						\
> 	}						\
>        else						\
> 	{						\
> 	  if (__builtin_isnan (r##n))			\
> 	    __builtin_abort ();				\
> 	  if (r##n < rl || r##n > ru)			\
> 	    __builtin_abort ();				\
> 	}						\
>      }
> #define TEST2(n, xl, xu, rl, ru) \
>    if (__builtin_isinf (xl) > 0				\
>        ? x > 0.0 && __builtin_isinf (x)			\
>        : __builtin_isinf (xu) < 0			\
>        ? x < 0.0 && __builtin_isinf (x)			\
>        : x >= xl && x <= xu				\
> 	&& (xl != 0.0					\
> 	    || __builtin_signbit (xl)			\
> 	    || !__builtin_signbit (x))			\
> 	&& (xu != 0.0					\
> 	    || !__builtin_signbit (xu)			\
> 	    || __builtin_signbit (x)))			\
>      {							\
>        double s##n = x * x;				\
>        if (__builtin_isnan (s##n))			\
> 	__builtin_abort ();				\
>        if (s##n < rl || s##n > ru)			\
> 	__builtin_abort ();				\
>      }
>    TEST (1,	2.0, 4.0,	6.0, 8.0,	12.0, 32.0, 0);
>    TEST (2,	-2.0, 3.0,	-7.0, 4.0,	-21.0, 14.0, 0);
>    TEST (3,	-9.0, 5.0,	8.0, 10.0,	-90.0, 50.0, 0);
>    TEST (4,	-0.0, 0.0,	16.0, 32.0,	-0.0, 0.0, 0);
>    TEST (5,	-0.0, 0.0,	16.0, 32.0,	-0.0, 0.0, 0);
>    TEST (6,	0.0, 0.0,	16.0, 32.0,	0.0, 0.0, 0);
>    TEST (7,	0.0, 0.0,	16.0, 32.0,	0.0, 0.0, 0);
>    TEST (8,	-0.0, -0.0,	16.0, 32.0,	-0.0, -0.0, 0);
>    TEST (9,	-0.0, -0.0,	16.0, 32.0,	-0.0, -0.0, 0);
>    TEST (10,	minf, inf,	minf, inf,	minf, inf, 1);
>    TEST (11,	-0.0, 0.0,	128.0, inf,	-0.0, 0.0, 1);
>    TEST (12,	-0.0, 0.0,	inf, inf,	0.0, 0.0, 2);
>    TEST (13,	minf, minf,	0.0, 0.0,	0.0, 0.0, 2);
>    TEST (14,	0.0, 2.0,	inf, inf,	inf, inf, 1);
>    TEST (15,	0.0, 2.0,	minf, minf,	minf, minf, 1);
>    TEST (16,	inf, inf,	-0.0, 2.0,	minf, inf, 1);
>    TEST (17,	minf, minf,	-0.0, 3.0,	minf, inf, 1);
> 
>    TEST2 (1,	2.0, 4.0,			4.0, 16.0);
>    TEST2 (2,	-0.0, 0.0,			-0.0, 0.0);
>    TEST2 (3,	0.0, inf,			0.0, inf);
>    TEST2 (4,	inf, inf,			inf, inf);
> }
> 
> 
> 	Jakub
> 

Re: [PATCH] range-op: Implement floating point multiplication fold_range [PR107569]

[Jakub Jelinek]

On Fri, Nov 11, 2022 at 11:01:38AM +0100, Aldy Hernandez wrote:
> > I've tried following, but it suffers from various issues:
> > 1) we don't handle __builtin_signbit (whatever) == 0 (or != 0) as guarantee
> >     that in the guarded code whatever has signbit 0 or 1
> 
> We have a range-op entry for __builtin_signbit in cfn_signbit.  Is this a
> shortcoming of this code, or something else?

Dunno, I admit I haven't investigated it much.  I just saw it when putting
a breakpoint on the mult fold_range.

> > 2) __builtin_isinf (x) > 0 is lowered to x > DBL_MAX, but unfortunately we don't
> >     infer from that [INF,INF] range, but [DBL_MAX, INF] range
> > 3) what I wrote above, I think we don't handle [0, 2] * [INF, INF] right but
> >     due to 2) we can't see it
> 
> Doesn't this boil down to a representation issue?  I wonder if we should
> bite the bullet and tweak build_gt() and build_lt() to represent open
> ranges.  In theory it should be one more/less ULP, while adjusting for
> HONOR_INFINITIES.

At least with the exception of MODE_COMPOSITE_P, I think we don't need
to introduce open ranges (and who cares about MODE_COMPOSITE_P if it is
conservatively correct).
For other floats, I think
x > cst
is always equivalent to
x >= nextafter (cst, inf)
and
x < cst
is always equivalent to
x <= nextafter (cst, -inf)
except for the signed zero cases which needs tiny bit more thought.
So, if we have
if (x > DBL_MAX)
then in code guarded by that we can just use [INF, INF] as range.

> If the signbit issue were resolved and we could represent > and < properly,
> would that allow us to write proper testcases without having to writing a
> plug-in (which I assume is a lot harder)?

I don't know, we'd need to see.
First work out on all the issues that result on the testcase the operand
ranges aren't exactly what we want (whether on the testcase side or on the
range-ops side or wherever) and once that looks ok, see if the ranges
on the rN/sN vars are correct and if so, watch what hasn't been folded away
and why.
I think the plugin would be 100-200 lines of code and then we could just
write multiple testcases against the plugin.

	Jakub

Re: [PATCH] range-op: Implement floating point multiplication fold_range [PR107569]

[Aldy Hernandez]

On 11/11/22 11:47, Jakub Jelinek wrote:
> On Fri, Nov 11, 2022 at 11:01:38AM +0100, Aldy Hernandez wrote:
>>> I've tried following, but it suffers from various issues:
>>> 1) we don't handle __builtin_signbit (whatever) == 0 (or != 0) as guarantee
>>>      that in the guarded code whatever has signbit 0 or 1
>>
>> We have a range-op entry for __builtin_signbit in cfn_signbit.  Is this a
>> shortcoming of this code, or something else?
> 
> Dunno, I admit I haven't investigated it much.  I just saw it when putting
> a breakpoint on the mult fold_range.

Could you send me a small testcase.  I can look into that.

> 
>>> 2) __builtin_isinf (x) > 0 is lowered to x > DBL_MAX, but unfortunately we don't
>>>      infer from that [INF,INF] range, but [DBL_MAX, INF] range
>>> 3) what I wrote above, I think we don't handle [0, 2] * [INF, INF] right but
>>>      due to 2) we can't see it
>>
>> Doesn't this boil down to a representation issue?  I wonder if we should
>> bite the bullet and tweak build_gt() and build_lt() to represent open
>> ranges.  In theory it should be one more/less ULP, while adjusting for
>> HONOR_INFINITIES.
> 
> At least with the exception of MODE_COMPOSITE_P, I think we don't need
> to introduce open ranges (and who cares about MODE_COMPOSITE_P if it is
> conservatively correct).
> For other floats, I think
> x > cst
> is always equivalent to
> x >= nextafter (cst, inf)
> and
> x < cst
> is always equivalent to
> x <= nextafter (cst, -inf)
> except for the signed zero cases which needs tiny bit more thought.
> So, if we have
> if (x > DBL_MAX)
> then in code guarded by that we can just use [INF, INF] as range.

Yeah, yeah.  That's exactly what I meant... using nextafter.  I'll look 
into that, as there seems there's more than one issue related to our 
lack of precision in representing < and >.

> 
>> If the signbit issue were resolved and we could represent > and < properly,
>> would that allow us to write proper testcases without having to writing a
>> plug-in (which I assume is a lot harder)?
> 
> I don't know, we'd need to see.
> First work out on all the issues that result on the testcase the operand
> ranges aren't exactly what we want (whether on the testcase side or on the
> range-ops side or wherever) and once that looks ok, see if the ranges
> on the rN/sN vars are correct and if so, watch what hasn't been folded away
> and why.
> I think the plugin would be 100-200 lines of code and then we could just
> write multiple testcases against the plugin.

If you think the plug-in will get better test coverage, by all means.  I 
was just trying to save you/us some work.

Andrew, do you have any thoughts on the plug-in?

Aldy