Implement range-op entry for sin/cos

On Tue, Apr 18, 2023 at 03:12:50PM +0200, Aldy Hernandez wrote:
> [I don't know why I keep poking at floats.  I must really like the pain.
>
> This is the range-op entry for sin/cos.  It is meant to serve as an
> example of what we can do for glibc math functions.  It is by no means
> exhaustive, just a stub to restrict the return range from sin/cos to
> [-1.0, 1.0] with appropriate smarts of NANs.
>
> As can be seen in the testcase, we see sin() as well as
> __builtin_sin() in the IL, and can resolve the resulting range
> accordingly.

Here is an updated version of the patch on top of the
Add targetm.libm_function_max_error
patch with all my comments incorporated into your patch (but still no
handling of sin/cos ranges shorter than 2*M_PI).

2023-04-28  Aldy Hernandez  <aldyh@redhat.com>
	    Jakub Jelinek  <jakub@redhat.com>

	* value-range.h (frange_nextafter): Declare.
	* gimple-range-op.cc (class cfn_sincos): New.
	(op_cfn_sin, op_cfn_cos): New variables.
	(gimple_range_op_handler::maybe_builtin_call): Handle
	CASE_CFN_{SIN,COS}{,_FN}.

	* gcc.dg/tree-ssa/range-sincos.c: New test.

3 files changed, 143 insertions, 0 deletions

diff --git a/gcc/gimple-range-op.cc b/gcc/gimple-range-op.cc
index 04e27d6..b66a4a0 100644
--- a/gcc/gimple-range-op.cc
+++ b/gcc/gimple-range-op.cc
@@ -401,6 +401,89 @@ public:
   }
 } op_cfn_copysign;
 
+class cfn_sincos : public range_operator_float
+{
+public:
+  using range_operator_float::fold_range;
+  using range_operator_float::op1_range;
+  cfn_sincos (combined_fn cfn) { m_cfn = cfn; }
+  virtual bool fold_range (frange &r, tree type,
+			   const frange &lh, const frange &,
+			   relation_trio) const final override
+  {
+    if (lh.undefined_p ())
+      return false;
+    if (lh.known_isnan () || lh.known_isinf ())
+      {
+	r.set_nan (type);
+	return true;
+      }
+    unsigned bulps = targetm.libm_function_max_error (m_cfn, TYPE_MODE (type),
+						      true);
+    if (bulps == ~0U)
+      r.set_varying (type);
+    else if (bulps == 0)
+      r.set (type, dconstm1, dconst1);
+    else
+      {
+	REAL_VALUE_TYPE boundmin, boundmax;
+	boundmax = dconst1;
+	while (bulps--)
+	  frange_nextafter (TYPE_MODE (type), boundmax, dconstinf);
+	real_arithmetic (&boundmin, NEGATE_EXPR, &boundmax, NULL);
+	r.set (type, boundmin, boundmax);
+      }
+    if (!lh.maybe_isnan () && !lh.maybe_isinf ())
+      r.clear_nan ();
+    return true;
+  }
+  virtual bool op1_range (frange &r, tree type,
+			  const frange &lhs, const frange &,
+			  relation_trio) const final override
+  {
+    if (lhs.undefined_p ())
+      return false;
+
+    // A known NAN means the input is [-INF,-INF][+INF,+INF] U +-NAN,
+    // which we can't currently represent.
+    if (lhs.known_isnan ())
+      {
+	r.set_varying (type);
+	return true;
+      }
+
+    // Results outside of [-1.0, +1.0] are impossible.
+    REAL_VALUE_TYPE lb = lhs.lower_bound ();
+    REAL_VALUE_TYPE ub = lhs.upper_bound ();
+    if (real_less (&ub, &dconstm1) || real_less (&dconst1, &lb))
+      {
+	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,
+	     [-INF,-INF][+INF,+INF] U +-NAN.  */
+	  r.set_varying (type);
+	return true;
+      }
+
+    if (!lhs.maybe_isnan ())
+      {
+	// If NAN is not valid result, the input cannot include either
+	// a NAN nor a +-INF.
+	lb = real_min_representable (type);
+	ub = real_max_representable (type);
+	r.set (type, lb, ub, nan_state (false, false));
+	return true;
+      }
+
+    r.set_varying (type);
+    return true;
+  }
+private:
+  combined_fn m_cfn;
+} op_cfn_sin (CFN_SIN), op_cfn_cos (CFN_COS);
+
 // Implement range operator for CFN_BUILT_IN_TOUPPER and CFN_BUILT_IN_TOLOWER.
 class cfn_toupper_tolower : public range_operator
 {
@@ -879,6 +962,20 @@ gimple_range_op_handler::maybe_builtin_call ()
       m_valid = true;
       break;
 
+    CASE_CFN_SIN:
+    CASE_CFN_SIN_FN:
+      m_op1 = gimple_call_arg (call, 0);
+      m_float = &op_cfn_sin;
+      m_valid = true;
+      break;
+
+    CASE_CFN_COS:
+    CASE_CFN_COS_FN:
+      m_op1 = gimple_call_arg (call, 0);
+      m_float = &op_cfn_cos;
+      m_valid = true;
+      break;
+
     case CFN_BUILT_IN_TOUPPER:
     case CFN_BUILT_IN_TOLOWER:
       // Only proceed If the argument is compatible with the LHS.
diff --git a/gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c b/gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c
new file mode 100644
index 0000000..337f9cd
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c
@@ -0,0 +1,43 @@
+// { 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 (__builtin_isnan (x))
+    __builtin_unreachable ();
+  x = sin (x);
+  if (x < -1.0 || x > 1.0)
+    link_error ();
+  use (x);
+}
+
+void
+bar (double x)
+{
+  if (!__builtin_isnan (sin (x)))
+    {
+      if (__builtin_isnan (x))
+	link_error ();
+      if (__builtin_isinf (x))
+	link_error ();
+    }
+}
+
+void
+stool (double x)
+{
+  double res1 = sin (x);
+  double res2 = __builtin_sin (x);
+  if (res1 < -1.0 || res2 < -1.0)
+    link_error ();
+  if (res1 > 1.0 || res2 > 1.0)
+    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 c37b6ef..0b61341 100644
--- a/gcc/value-range.h
+++ b/gcc/value-range.h
@@ -1388,4 +1388,7 @@ frange::nan_signbit_p (bool &signbit) const
   return true;
 }
 
+void frange_nextafter (enum machine_mode, REAL_VALUE_TYPE &,
+		       const REAL_VALUE_TYPE &);
+
 #endif // GCC_VALUE_RANGE_H