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

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,

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

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