range-op-float: Fix up reverse binary operations [PR109008]

  Hi!

The following testcase is reduced from miscompilation of scipy package.
If we have say lhs = [1., 1.] - [1., 1.] and want to compute the range
of lhs from it, we correctly determine it is [0., 0.] (if computations
are exact, we generally don't try to round them further in
frange_arithmetic).  In the testcase it is about a reverse operation,
[1., 1.] = op1 + [1., 1.] and we want to compute range of op1 from that.
Right now we just perform the inverse operation (there are some corner
cases about NaN and infinities handling) and so arrive to range
[0., 0.] as well, and because it is a singleton, optimize return eps;
to return 0.  That is incorrect though, for the reverse ops we need to
take into account also rounding, the right exact range is
[-0x1.0p-54, 0x1.0p-53] in this case when rounding to nearest, i.e.
all numbers which added to 1. with round to nearest still produce 1.

The problem isn't solely on singleton ranges, and isn't solely on
results around zero.  We basically need to consider also values
where the result is up to 0.5ulp away from the lhs range boundaries
in each direction.

The following patch fixes it by extending the lhs range for the
reverse operations by 1ulp in each direction.  The PR contains
a pseudo-random test generator I've used to generate 300000 tests
of + and - and then used the same test with * and / instead of + and -
together with a hack to print the discovered ranges by the patch in
a form that another test could then verify the range is conservatively
correct and how far it is from a minimal range.

I believe the results are good enough for now, though plan to look
incrementally into trying to do something better on the -XXX_MAX or
XXX_MAX boundaries (where I think frange_nextafter will use -inf or +inf)
and also try to increase the range just by 0.5ulp rather than 1ulp
if !flag_rounding_math.  But dunno if either of those will be doable
and will pass the testing, so I think it is worth committing this fix
first.

Bootstrapped/regtested on x86_64-linux and i686-linux, ok for trunk?

2023-03-09  Jakub Jelinek  <jakub@redhat.com>
	    Richard Biener  <rguenther@suse.de>

	PR tree-optimization/109008
	* range-op-float.cc (float_widen_lhs_range): New function.
	(foperator_plus::op1_range, foperator_minus::op1_range,
	foperator_minus::op2_range, foperator_mult::op1_range,
	foperator_div::op1_range, foperator_div::op2_range): Use it.

	* gcc.c-torture/execute/ieee/pr109008.c: New test.


	Jakub
  

On Thu, 9 Mar 2023, Jakub Jelinek wrote:

> Hi!
> 
> The following testcase is reduced from miscompilation of scipy package.
> If we have say lhs = [1., 1.] - [1., 1.] and want to compute the range
> of lhs from it, we correctly determine it is [0., 0.] (if computations
> are exact, we generally don't try to round them further in
> frange_arithmetic).  In the testcase it is about a reverse operation,
> [1., 1.] = op1 + [1., 1.] and we want to compute range of op1 from that.
> Right now we just perform the inverse operation (there are some corner
> cases about NaN and infinities handling) and so arrive to range
> [0., 0.] as well, and because it is a singleton, optimize return eps;
> to return 0.  That is incorrect though, for the reverse ops we need to
> take into account also rounding, the right exact range is
> [-0x1.0p-54, 0x1.0p-53] in this case when rounding to nearest, i.e.
> all numbers which added to 1. with round to nearest still produce 1.
> 
> The problem isn't solely on singleton ranges, and isn't solely on
> results around zero.  We basically need to consider also values
> where the result is up to 0.5ulp away from the lhs range boundaries
> in each direction.
> 
> The following patch fixes it by extending the lhs range for the
> reverse operations by 1ulp in each direction.  The PR contains
> a pseudo-random test generator I've used to generate 300000 tests
> of + and - and then used the same test with * and / instead of + and -
> together with a hack to print the discovered ranges by the patch in
> a form that another test could then verify the range is conservatively
> correct and how far it is from a minimal range.
> 
> I believe the results are good enough for now, though plan to look
> incrementally into trying to do something better on the -XXX_MAX or
> XXX_MAX boundaries (where I think frange_nextafter will use -inf or +inf)
> and also try to increase the range just by 0.5ulp rather than 1ulp
> if !flag_rounding_math.  But dunno if either of those will be doable
> and will pass the testing, so I think it is worth committing this fix
> first.

Sounds good.

> Bootstrapped/regtested on x86_64-linux and i686-linux, ok for trunk?

OK.

Thanks,
Richard.

> 2023-03-09  Jakub Jelinek  <jakub@redhat.com>
> 	    Richard Biener  <rguenther@suse.de>
> 
> 	PR tree-optimization/109008
> 	* range-op-float.cc (float_widen_lhs_range): New function.
> 	(foperator_plus::op1_range, foperator_minus::op1_range,
> 	foperator_minus::op2_range, foperator_mult::op1_range,
> 	foperator_div::op1_range, foperator_div::op2_range): Use it.
> 
> 	* gcc.c-torture/execute/ieee/pr109008.c: New test.
> 
> --- gcc/range-op-float.cc.jj	2023-03-08 12:33:44.641043477 +0100
> +++ gcc/range-op-float.cc	2023-03-08 13:13:09.015341002 +0100
> @@ -2199,6 +2199,33 @@ zero_to_inf_range (REAL_VALUE_TYPE &lb,
>      }
>  }
>  
> +/* Extend the LHS range by 1ulp in each direction.  For op1_range
> +   or op2_range of binary operations just computing the inverse
> +   operation on ranges isn't sufficient.  Consider e.g.
> +   [1., 1.] = op1 + [1., 1.].  op1's range is not [0., 0.], but
> +   [-0x1.0p-54, 0x1.0p-53] (when not -frounding-math), any value for
> +   which adding 1. to it results in 1. after rounding to nearest.
> +   So, for op1_range/op2_range extend the lhs range by 1ulp in each
> +   direction.  See PR109008 for more details.  */
> +
> +static frange
> +float_widen_lhs_range (tree type, const frange &lhs)
> +{
> +  frange ret = lhs;
> +  if (lhs.known_isnan ())
> +    return ret;
> +  REAL_VALUE_TYPE lb = lhs.lower_bound ();
> +  REAL_VALUE_TYPE ub = lhs.upper_bound ();
> +  if (real_isfinite (&lb))
> +    frange_nextafter (TYPE_MODE (type), lb, dconstninf);
> +  if (real_isfinite (&ub))
> +    frange_nextafter (TYPE_MODE (type), ub, dconstinf);
> +  ret.set (type, lb, ub);
> +  ret.clear_nan ();
> +  ret.union_ (lhs);
> +  return ret;
> +}
> +
>  class foperator_plus : public range_operator_float
>  {
>    using range_operator_float::op1_range;
> @@ -2214,8 +2241,9 @@ public:
>      range_op_handler minus (MINUS_EXPR, type);
>      if (!minus)
>        return false;
> -    return float_binary_op_range_finish (minus.fold_range (r, type, lhs, op2),
> -					 r, type, lhs);
> +    frange wlhs = float_widen_lhs_range (type, lhs);
> +    return float_binary_op_range_finish (minus.fold_range (r, type, wlhs, op2),
> +					 r, type, wlhs);
>    }
>    virtual bool op2_range (frange &r, tree type,
>  			  const frange &lhs,
> @@ -2260,9 +2288,10 @@ public:
>    {
>      if (lhs.undefined_p ())
>        return false;
> -    return float_binary_op_range_finish (fop_plus.fold_range (r, type, lhs,
> +    frange wlhs = float_widen_lhs_range (type, lhs);
> +    return float_binary_op_range_finish (fop_plus.fold_range (r, type, wlhs,
>  							      op2),
> -					 r, type, lhs);
> +					 r, type, wlhs);
>    }
>    virtual bool op2_range (frange &r, tree type,
>  			  const frange &lhs,
> @@ -2271,8 +2300,9 @@ public:
>    {
>      if (lhs.undefined_p ())
>        return false;
> -    return float_binary_op_range_finish (fold_range (r, type, op1, lhs),
> -					 r, type, lhs);
> +    frange wlhs = float_widen_lhs_range (type, lhs);
> +    return float_binary_op_range_finish (fold_range (r, type, op1, wlhs),
> +					 r, type, wlhs);
>    }
>  private:
>    void rv_fold (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub, bool &maybe_nan,
> @@ -2338,13 +2368,14 @@ public:
>      range_op_handler rdiv (RDIV_EXPR, type);
>      if (!rdiv)
>        return false;
> -    bool ret = rdiv.fold_range (r, type, lhs, op2);
> +    frange wlhs = float_widen_lhs_range (type, lhs);
> +    bool ret = rdiv.fold_range (r, type, wlhs, op2);
>      if (ret == false)
>        return false;
> -    if (lhs.known_isnan () || op2.known_isnan () || op2.undefined_p ())
> -      return float_binary_op_range_finish (ret, r, type, lhs);
> -    const REAL_VALUE_TYPE &lhs_lb = lhs.lower_bound ();
> -    const REAL_VALUE_TYPE &lhs_ub = lhs.upper_bound ();
> +    if (wlhs.known_isnan () || op2.known_isnan () || op2.undefined_p ())
> +      return float_binary_op_range_finish (ret, r, type, wlhs);
> +    const REAL_VALUE_TYPE &lhs_lb = wlhs.lower_bound ();
> +    const REAL_VALUE_TYPE &lhs_ub = wlhs.upper_bound ();
>      const REAL_VALUE_TYPE &op2_lb = op2.lower_bound ();
>      const REAL_VALUE_TYPE &op2_ub = op2.upper_bound ();
>      if ((contains_zero_p (lhs_lb, lhs_ub) && contains_zero_p (op2_lb, op2_ub))
> @@ -2363,7 +2394,7 @@ public:
>      // or if lhs must be zero and op2 doesn't include zero, it would be
>      // UNDEFINED, while rdiv.fold_range computes a zero or singleton INF
>      // range.  Those are supersets of UNDEFINED, so let's keep that way.
> -    return float_binary_op_range_finish (ret, r, type, lhs);
> +    return float_binary_op_range_finish (ret, r, type, wlhs);
>    }
>    virtual bool op2_range (frange &r, tree type,
>  			  const frange &lhs,
> @@ -2490,13 +2521,14 @@ public:
>    {
>      if (lhs.undefined_p ())
>        return false;
> -    bool ret = fop_mult.fold_range (r, type, lhs, op2);
> +    frange wlhs = float_widen_lhs_range (type, lhs);
> +    bool ret = fop_mult.fold_range (r, type, wlhs, op2);
>      if (!ret)
>        return ret;
> -    if (lhs.known_isnan () || op2.known_isnan () || op2.undefined_p ())
> -      return float_binary_op_range_finish (ret, r, type, lhs);
> -    const REAL_VALUE_TYPE &lhs_lb = lhs.lower_bound ();
> -    const REAL_VALUE_TYPE &lhs_ub = lhs.upper_bound ();
> +    if (wlhs.known_isnan () || op2.known_isnan () || op2.undefined_p ())
> +      return float_binary_op_range_finish (ret, r, type, wlhs);
> +    const REAL_VALUE_TYPE &lhs_lb = wlhs.lower_bound ();
> +    const REAL_VALUE_TYPE &lhs_ub = wlhs.upper_bound ();
>      const REAL_VALUE_TYPE &op2_lb = op2.lower_bound ();
>      const REAL_VALUE_TYPE &op2_ub = op2.upper_bound ();
>      if ((contains_zero_p (lhs_lb, lhs_ub)
> @@ -2512,7 +2544,7 @@ public:
>  	zero_to_inf_range (lb, ub, signbit_known);
>  	r.set (type, lb, ub);
>        }
> -    return float_binary_op_range_finish (ret, r, type, lhs);
> +    return float_binary_op_range_finish (ret, r, type, wlhs);
>    }
>    virtual bool op2_range (frange &r, tree type,
>  			  const frange &lhs,
> @@ -2521,13 +2553,14 @@ public:
>    {
>      if (lhs.undefined_p ())
>        return false;
> -    bool ret = fold_range (r, type, op1, lhs);
> +    frange wlhs = float_widen_lhs_range (type, lhs);
> +    bool ret = fold_range (r, type, op1, wlhs);
>      if (!ret)
>        return ret;
> -    if (lhs.known_isnan () || op1.known_isnan () || op1.undefined_p ())
> -      return float_binary_op_range_finish (ret, r, type, lhs, true);
> -    const REAL_VALUE_TYPE &lhs_lb = lhs.lower_bound ();
> -    const REAL_VALUE_TYPE &lhs_ub = lhs.upper_bound ();
> +    if (wlhs.known_isnan () || op1.known_isnan () || op1.undefined_p ())
> +      return float_binary_op_range_finish (ret, r, type, wlhs, true);
> +    const REAL_VALUE_TYPE &lhs_lb = wlhs.lower_bound ();
> +    const REAL_VALUE_TYPE &lhs_ub = wlhs.upper_bound ();
>      const REAL_VALUE_TYPE &op1_lb = op1.lower_bound ();
>      const REAL_VALUE_TYPE &op1_ub = op1.upper_bound ();
>      if ((contains_zero_p (lhs_lb, lhs_ub) && contains_zero_p (op1_lb, op1_ub))
> @@ -2542,7 +2575,7 @@ public:
>  	zero_to_inf_range (lb, ub, signbit_known);
>  	r.set (type, lb, ub);
>        }
> -    return float_binary_op_range_finish (ret, r, type, lhs, true);
> +    return float_binary_op_range_finish (ret, r, type, wlhs, true);
>    }
>  private:
>    void rv_fold (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub, bool &maybe_nan,
> --- gcc/testsuite/gcc.c-torture/execute/ieee/pr109008.c.jj	2023-03-08 21:30:19.158618157 +0100
> +++ gcc/testsuite/gcc.c-torture/execute/ieee/pr109008.c	2023-03-08 21:29:49.899039474 +0100
> @@ -0,0 +1,18 @@
> +/* PR tree-optimization/109008 */
> +
> +__attribute__((noipa)) double
> +foo (double eps)
> +{
> +  double d = 1. + eps;
> +  if (d == 1.)
> +    return eps;
> +  return 0.0;
> +}
> +
> +int
> +main ()
> +{
> +  if (foo (__DBL_EPSILON__ / 8.0) == 0.0)
> +    __builtin_abort ();
> +  return 0;
> +}
> 
> 	Jakub
> 
>