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

[Aldy Hernandez <aldyh@redhat.com>]

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