[PATCH] gimple-range-op: Improve handling of sqrt ranges
Aldy Hernandez
aldyh@redhat.com
Fri May 5 09:06:31 GMT 2023
On 5/5/23 10:00, Jakub Jelinek wrote:
> Hi!
>
> 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.
>
> Bootstrapped/regtested on x86_64-linux and i686-linux, ok for trunk?
>
> 2023-05-05 Jakub Jelinek <jakub@redhat.com>
>
> * 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/value-range.h.jj 2023-05-04 13:34:45.140336099 +0200
> +++ gcc/value-range.h 2023-05-04 16:28:18.286108178 +0200
> @@ -1294,5 +1294,8 @@ frange::nan_signbit_p (bool &signbit) co
>
> 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
> --- gcc/range-op-float.cc.jj 2023-05-04 13:34:45.139336114 +0200
> +++ gcc/range-op-float.cc 2023-05-04 16:28:18.285108192 +0200
> @@ -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,
> --- gcc/gimple-range-op.cc.jj 2023-05-04 13:34:45.139336114 +0200
> +++ gcc/gimple-range-op.cc 2023-05-04 19:58:44.842606865 +0200
> @@ -44,6 +44,7 @@ along with GCC; see the file COPYING3.
> #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)
> +{
Since you're returning a range of sorts [low, high], would it be cleaner
to return an frange, or is always calculating low/high too expensive? I
notice you avoid it when passing NULL.
Would you mind adding a typedef for the (*func) callback above? I
always find C callbacks a pain to read.
Thanks.
Aldy
> + 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,20 @@ 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);
> + r.intersect (r2);
> return true;
> }
> virtual bool op1_range (frange &r, tree type,
> @@ -455,27 +530,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));
> + else
> + 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 (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))
> {
> - // 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;
> + 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;
> }
> -
> - r.set_varying (type);
> + else
> + ub = dconstinf;
> + frange r2;
> + r2.set (type, lb, ub);
> + r.intersect (r2);
> return true;
> }
> } op_cfn_sqrt;
> --- gcc/testsuite/gcc.dg/tree-ssa/range-sqrt-2.c.jj 2023-05-04 20:05:26.780872035 +0200
> +++ gcc/testsuite/gcc.dg/tree-ssa/range-sqrt-2.c 2023-05-04 19:50:28.624689220 +0200
> @@ -0,0 +1,44 @@
> +// { dg-do compile }
> +// { dg-options "-O2 -fdump-tree-evrp -fno-thread-jumps" }
> +
> +#include <math.h>
> +
> +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 } } } } }
>
> Jakub
>
More information about the Gcc-patches
mailing list