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

[Richard Biener <rguenther@suse.de>]

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

-- 
Richard Biener <rguenther@suse.de>
SUSE Software Solutions Germany GmbH, Frankenstrasse 146, 90461 Nuernberg,
Germany; GF: Ivo Totev, Andrew Myers, Andrew McDonald, Boudien Moerman;
HRB 36809 (AG Nuernberg)