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

  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.



	Jakub
  

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
  

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
  

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