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 <ja...@redhat.com>
> Richard Biener <rguent...@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 <rguent...@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)
