From 82aef047eda376f41eccf51ba1277e46160f6420 Mon Sep 17 00:00:00 2001 From: Jakub Jelinek Date: Sat, 6 May 2023 02:35:02 +0200 Subject: [PATCH] gimple-range-op: Improve handling of sqrt ranges The previous patch just added basic intrinsic ranges for sqrt ([-0.0, +Inf] +-NAN being the general result range of the function and [-0.0, +Inf] the general operand range if result isn't NAN etc.), the following patch intersects those ranges with particular range computed from argument or result's exact range with the expected error in ulps taken into account and adds a function (frange_arithmetic variant) which can be used by other functions as well as helper. 2023-05-06 Jakub Jelinek * value-range.h (frange_arithmetic): Declare. * range-op-float.cc (frange_arithmetic): No longer static. * gimple-range-op.cc (frange_mpfr_arg1): New function. (cfn_sqrt::fold_range): Intersect the generic boundaries range with range computed from sqrt of the particular bounds. (cfn_sqrt::op1_range): Intersect the generic boundaries range with range computed from squared particular bounds. * gcc.dg/tree-ssa/range-sqrt-2.c: New test. --- gcc/gimple-range-op.cc | 151 +++++++++++++++++-- gcc/range-op-float.cc | 2 +- gcc/testsuite/gcc.dg/tree-ssa/range-sqrt-2.c | 44 ++++++ gcc/value-range.h | 3 + 4 files changed, 183 insertions(+), 17 deletions(-) create mode 100644 gcc/testsuite/gcc.dg/tree-ssa/range-sqrt-2.c diff --git a/gcc/gimple-range-op.cc b/gcc/gimple-range-op.cc index 5d1f921ba40..fe11128c07f 100644 --- a/gcc/gimple-range-op.cc +++ b/gcc/gimple-range-op.cc @@ -44,6 +44,7 @@ along with GCC; see the file COPYING3. If not see #include "value-query.h" #include "gimple-range.h" #include "attr-fnspec.h" +#include "realmpfr.h" // Given stmt S, fill VEC, up to VEC_SIZE elements, with relevant ssa-names // on the statement. For efficiency, it is an error to not pass in enough @@ -403,6 +404,66 @@ public: } } op_cfn_copysign; +/* Compute FUNC (ARG) where FUNC is a mpfr function. If RES_LOW is non-NULL, + set it to low bound of possible range if the function is expected to have + ULPS precision and similarly if RES_HIGH is non-NULL, set it to high bound. + If the function returns false, the results weren't set. */ + +static bool +frange_mpfr_arg1 (REAL_VALUE_TYPE *res_low, REAL_VALUE_TYPE *res_high, + int (*func) (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t), + const REAL_VALUE_TYPE &arg, tree type, unsigned ulps) +{ + if (ulps == ~0U || !real_isfinite (&arg)) + return false; + machine_mode mode = TYPE_MODE (type); + const real_format *format = REAL_MODE_FORMAT (mode); + auto_mpfr m (format->p); + mpfr_from_real (m, &arg, MPFR_RNDN); + mpfr_clear_flags (); + bool inexact = func (m, m, MPFR_RNDN); + if (!mpfr_number_p (m) || mpfr_overflow_p () || mpfr_underflow_p ()) + return false; + + REAL_VALUE_TYPE value, result; + real_from_mpfr (&value, m, format, MPFR_RNDN); + if (!real_isfinite (&value)) + return false; + if ((value.cl == rvc_zero) != (mpfr_zero_p (m) != 0)) + inexact = true; + + real_convert (&result, format, &value); + if (!real_isfinite (&result)) + return false; + bool round_low = false; + bool round_high = false; + if (!ulps && flag_rounding_math) + ++ulps; + if (inexact || !real_identical (&result, &value)) + { + if (MODE_COMPOSITE_P (mode)) + round_low = round_high = true; + else + { + round_low = !real_less (&result, &value); + round_high = !real_less (&value, &result); + } + } + if (res_low) + { + *res_low = result; + for (unsigned int i = 0; i < ulps + round_low; ++i) + frange_nextafter (mode, *res_low, dconstninf); + } + if (res_high) + { + *res_high = result; + for (unsigned int i = 0; i < ulps + round_high; ++i) + frange_nextafter (mode, *res_high, dconstinf); + } + return true; +} + class cfn_sqrt : public range_operator_float { public: @@ -434,6 +495,21 @@ public: } if (!lh.maybe_isnan () && !real_less (&lh.lower_bound (), &dconst0)) r.clear_nan (); + + unsigned ulps + = targetm.libm_function_max_error (CFN_SQRT, TYPE_MODE (type), false); + if (ulps == ~0U) + return true; + REAL_VALUE_TYPE lb = lh.lower_bound (); + REAL_VALUE_TYPE ub = lh.upper_bound (); + if (!frange_mpfr_arg1 (&lb, NULL, mpfr_sqrt, lb, type, ulps)) + lb = dconstninf; + if (!frange_mpfr_arg1 (NULL, &ub, mpfr_sqrt, ub, type, ulps)) + ub = dconstinf; + frange r2; + r2.set (type, lb, ub); + r2.flush_denormals_to_zero (); + r.intersect (r2); return true; } virtual bool op1_range (frange &r, tree type, @@ -455,27 +531,70 @@ public: } // Results outside of [-0.0, +Inf] are impossible. - const REAL_VALUE_TYPE &ub = lhs.upper_bound (); - if (real_less (&ub, &dconstm0)) + unsigned bulps + = targetm.libm_function_max_error (CFN_SQRT, TYPE_MODE (type), true); + if (bulps != ~0U) { - if (!lhs.maybe_isnan ()) - r.set_undefined (); - else - // If lhs could be NAN and finite result is impossible, - // the range is like lhs.known_isnan () above. - goto known_nan; - return true; + const REAL_VALUE_TYPE &ub = lhs.upper_bound (); + REAL_VALUE_TYPE m0 = dconstm0; + while (bulps--) + frange_nextafter (TYPE_MODE (type), m0, dconstninf); + if (real_less (&ub, &m0)) + { + if (!lhs.maybe_isnan ()) + r.set_undefined (); + else + // If lhs could be NAN and finite result is impossible, + // the range is like lhs.known_isnan () above. + goto known_nan; + return true; + } } if (!lhs.maybe_isnan ()) - { - // If NAN is not valid result, the input cannot include either - // a NAN nor values smaller than -0. - r.set (type, dconstm0, dconstinf, nan_state (false, false)); - return true; - } + // If NAN is not valid result, the input cannot include either + // a NAN nor values smaller than -0. + r.set (type, dconstm0, dconstinf, nan_state (false, false)); + else + r.set_varying (type); - r.set_varying (type); + unsigned ulps + = targetm.libm_function_max_error (CFN_SQRT, TYPE_MODE (type), false); + if (ulps == ~0U) + return true; + REAL_VALUE_TYPE lb = lhs.lower_bound (); + REAL_VALUE_TYPE ub = lhs.upper_bound (); + if (!lhs.maybe_isnan () && real_less (&dconst0, &lb)) + { + for (unsigned i = 0; i < ulps; ++i) + frange_nextafter (TYPE_MODE (type), lb, dconstninf); + if (real_less (&dconst0, &lb)) + { + REAL_VALUE_TYPE op = lb; + frange_arithmetic (MULT_EXPR, type, lb, op, op, dconstninf); + } + else + lb = dconstninf; + } + else + lb = dconstninf; + if (real_isfinite (&ub) && real_less (&dconst0, &ub)) + { + for (unsigned i = 0; i < ulps; ++i) + frange_nextafter (TYPE_MODE (type), ub, dconstinf); + if (real_isfinite (&ub)) + { + REAL_VALUE_TYPE op = ub; + frange_arithmetic (MULT_EXPR, type, ub, op, op, dconstinf); + } + else + ub = dconstinf; + } + else + ub = dconstinf; + frange r2; + r2.set (type, lb, ub); + r.intersect (r2); return true; } } op_cfn_sqrt; diff --git a/gcc/range-op-float.cc b/gcc/range-op-float.cc index 9d184611ae9..a99a6b01ed8 100644 --- a/gcc/range-op-float.cc +++ b/gcc/range-op-float.cc @@ -305,7 +305,7 @@ frange_nextafter (enum machine_mode mode, // SF/DFmode (when storing into memory from the 387 stack). Maybe // this is ok as well though it is just occasionally more precise. ?? -static void +void frange_arithmetic (enum tree_code code, tree type, REAL_VALUE_TYPE &result, const REAL_VALUE_TYPE &op1, diff --git a/gcc/testsuite/gcc.dg/tree-ssa/range-sqrt-2.c b/gcc/testsuite/gcc.dg/tree-ssa/range-sqrt-2.c new file mode 100644 index 00000000000..a146d9d14a5 --- /dev/null +++ b/gcc/testsuite/gcc.dg/tree-ssa/range-sqrt-2.c @@ -0,0 +1,44 @@ +// { dg-do compile } +// { dg-options "-O2 -fdump-tree-evrp -fno-thread-jumps" } + +#include + +void use (double); +void link_error (); + +void +foo (double x) +{ + if (x < 1.0 || x > 9.0) + __builtin_unreachable (); + x = sqrt (x); + if (x < 0.875 || x > 3.125) + link_error (); + use (x); +} + +void +bar (double x) +{ + if (sqrt (x) >= 2.0 && sqrt (x) <= 4.0) + { + if (__builtin_isnan (x)) + link_error (); + if (x < 3.875 || x > 16.125) + link_error (); + } +} + +void +stool (double x) +{ + if (x >= 64.0) + { + double res1 = sqrt (x); + double res2 = __builtin_sqrt (x); + if (res1 < 7.875 || res2 < 7.875) + link_error (); + } +} + +// { dg-final { scan-tree-dump-not "link_error" "evrp" { target { { *-*-linux* } && { glibc } } } } } diff --git a/gcc/value-range.h b/gcc/value-range.h index ab8d7d3c5dc..22b0250b11b 100644 --- a/gcc/value-range.h +++ b/gcc/value-range.h @@ -1294,5 +1294,8 @@ frange::nan_signbit_p (bool &signbit) const void frange_nextafter (enum machine_mode, REAL_VALUE_TYPE &, const REAL_VALUE_TYPE &); +void frange_arithmetic (enum tree_code, tree, REAL_VALUE_TYPE &, + const REAL_VALUE_TYPE &, const REAL_VALUE_TYPE &, + const REAL_VALUE_TYPE &); #endif // GCC_VALUE_RANGE_H