[PATCH] gimple-range-op: Improve handling of sqrt ranges

[PATCH] gimple-range-op: Improve handling of sqrt ranges

[Jakub Jelinek]

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)
+{
+  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

Re: [PATCH] gimple-range-op: Improve handling of sqrt ranges

[Aldy Hernandez]

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
> 

Re: [PATCH] gimple-range-op: Improve handling of sqrt ranges

[Jakub Jelinek]

On Fri, May 05, 2023 at 11:06:31AM +0200, Aldy Hernandez wrote:
> > +/* 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.

The point was that the caller can tell which bound it needs, low, high or
both and we don't waste time calculating ones we don't need (especially with
larger values of ulps).  E.g. for the sqrt case we only need one of them,
but when I thought about the sin/cos case, I'll probably need both and
calling the function twice would mean repeating the even more expensive mpfr
call.

> Would you mind adding a typedef for the (*func) callback above?  I always
> find C callbacks a pain to read.

I can, what I used comes from elsewhere (builtins.cc/fold-const-call.cc
which use it like that).

	Jakub

Re: [PATCH] gimple-range-op: Improve handling of sqrt ranges

[Aldy Hernandez]

On Fri, May 5, 2023 at 11:14 AM Jakub Jelinek <jakub@redhat.com> wrote:
>
> On Fri, May 05, 2023 at 11:06:31AM +0200, Aldy Hernandez wrote:
> > > +/* 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.
>
> The point was that the caller can tell which bound it needs, low, high or
> both and we don't waste time calculating ones we don't need (especially with
> larger values of ulps).  E.g. for the sqrt case we only need one of them,
> but when I thought about the sin/cos case, I'll probably need both and
> calling the function twice would mean repeating the even more expensive mpfr
> call.
>
> > Would you mind adding a typedef for the (*func) callback above?  I always
> > find C callbacks a pain to read.
>
> I can, what I used comes from elsewhere (builtins.cc/fold-const-call.cc
> which use it like that).

It would be my preference to have a typedef in
builtins.cc/fold-const-call.cc as well, as we could clean everything
up.  But I defer to you as a global maintainer, whether we want to do
that or not.  If not, then don't bother with just cleaning up
gimple-range-op.cc.

LGTM.

Aldy

Re: [PATCH] gimple-range-op: Improve handling of sqrt ranges

[Mikael Morin]

Hello,

Le 05/05/2023 à 10:00, Jakub Jelinek via Gcc-patches a écrit :
> 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?
> 
I think there is something wrong with the handling of NAN and lower 
bounds in op1_range.  If the lhs has positive lower bound and can be 
nan, op1 should have -inf lower bound regardless of lhs's lower bound value.

Maybe it would be less error prone to combine the ordered range part 
(r2) and the unordered one (r) with union instead of intersection?

The rest looks good.  There is the case of slightly negative lhs bounds 
coming (through ulps) from slightly positive op1 that could get better 
than infinite bounds for op1 range, but at least it is conservatively 
correct, and I'm not sure it's that important anyway.

Mikael

Re: [PATCH] gimple-range-op: Improve handling of sqrt ranges

[Jakub Jelinek]

On Fri, May 05, 2023 at 01:38:48PM +0200, Mikael Morin wrote:
> I think there is something wrong with the handling of NAN and lower bounds
> in op1_range.  If the lhs has positive lower bound and can be nan, op1
> should have -inf lower bound regardless of lhs's lower bound value.

Oops, you're right, will retest with
--- gcc/gimple-range-op.cc	2023-05-04 19:58:44.842606865 +0200
+++ gcc/gimple-range-op.cc	2023-05-05 13:53:58.742406160 +0200
@@ -508,6 +508,7 @@
       ub = dconstinf;
     frange r2;
     r2.set (type, lb, ub);
+    r2.flush_denormals_to_zero ();
     r.intersect (r2);
     return true;
   }
@@ -563,7 +564,7 @@
       return true;
     REAL_VALUE_TYPE lb = lhs.lower_bound ();
     REAL_VALUE_TYPE ub = lhs.upper_bound ();
-    if (real_less (&dconst0, &lb))
+    if (!lhs.maybe_isnan () && real_less (&dconst0, &lb))
       {
 	for (unsigned i = 0; i < ulps; ++i)
 	  frange_nextafter (TYPE_MODE (type), lb, dconstninf);
incremental change (the first hunk because fold_range is forward operation
and that could in some configurations flush denormals to zero).
Thanks for catching that.

> Maybe it would be less error prone to combine the ordered range part (r2)
> and the unordered one (r) with union instead of intersection?

The reason for the intersection rather than tweaking the bounds immediately
is for 2 reasons, ulps might be larger than bulps (the boundary ones) but
the boundary ulps should take priority on the boundaries, so if say bulps is
2ulps and ulps 10ulps and sqrt (lb) is 0.0 + 4ulps, then the range should be
[-0.0 - 2ulps, ...] rather than [-0.0 - 6ulps, ...], and to help readers
determine this is the range from the intrinsic boundaries, this is the range
from specific values, one can put breakpoints around the intersection and
print both ranges and the result etc.

> The rest looks good.  There is the case of slightly negative lhs bounds
> coming (through ulps) from slightly positive op1 that could get better than
> infinite bounds for op1 range, but at least it is conservatively correct,
> and I'm not sure it's that important anyway.

Currently bulps for sqrt in glibc is always 0, haven't discovered any case of
sqrt actually returning < -0.0 value, so the sqrt bulps other than ~0U
(don't know, currently all non-glibc targets) is always 0.  That is just in
case some other implementation returns such values.  For sin/cos
unfortunately even glibc sometimes returns > 1.0 or < -1.0 values in
non-default rounding modes.

	Jakub