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

[Aldy Hernandez <aldyh@redhat.com>]

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