gimple-range-op: Improve handling of sin/cos ranges

Similarly to the earlier sqrt patch, this patch attempts to improve
sin/cos ranges.  As the functions are periodic, for the reverse range
there is not much we can do (but I've discovered I forgot to take
into account the boundary ulps for the discovery of impossible result
ranges).  For fold_range, we can do something only if the range is
narrow enough (narrower than 2*pi).  The patch computes the value of
the functions (taking ulps into account) and also computes the derivative
to find out if the function is growing or declining on the boundaries and
from that it figures out if the result range should be
[min (fn (lb), fn (ub)), max (fn (lb), fn (ub))] or if it needs to be
extended to 1 (actually using +Inf) and/or -1 (actually using -Inf) because
there must be a local minimum and/or maximum in the range.

2023-05-06  Jakub Jelinek  <jakub@redhat.com>

	* real.h (dconst_pi): Define.
	(dconst_e_ptr): Formatting fix.
	(dconst_pi_ptr): Declare.
	* real.cc (dconst_pi_ptr): New function.
	* gimple-range-op.cc (cfn_sincos::fold_range): Intersect the generic
	boundaries range with range computed from sin/cos of the particular
	bounds if the argument range is shorter than 2*pi.
	(cfn_sincos::op1_range): Take bulps into account when determining
	which result ranges are always invalid or behave like known NAN.

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

4 files changed, 178 insertions, 17 deletions

diff --git a/gcc/gimple-range-op.cc b/gcc/gimple-range-op.cc
index fe11128..594bd30 100644
--- a/gcc/gimple-range-op.cc
+++ b/gcc/gimple-range-op.cc
@@ -633,6 +633,98 @@ public:
       }
     if (!lh.maybe_isnan () && !lh.maybe_isinf ())
       r.clear_nan ();
+
+    unsigned ulps
+      = targetm.libm_function_max_error (m_cfn, TYPE_MODE (type), false);
+    if (ulps == ~0U)
+      return true;
+    REAL_VALUE_TYPE lb = lh.lower_bound ();
+    REAL_VALUE_TYPE ub = lh.upper_bound ();
+    REAL_VALUE_TYPE diff;
+    real_arithmetic (&diff, MINUS_EXPR, &ub, &lb);
+    if (!real_isfinite (&diff))
+      return true;
+    REAL_VALUE_TYPE pi = dconst_pi ();
+    REAL_VALUE_TYPE pix2;
+    real_arithmetic (&pix2, PLUS_EXPR, &pi, &pi);
+    // We can only try to narrow the range further if ub-lb < 2*pi.
+    if (!real_less (&diff, &pix2))
+      return true;
+    REAL_VALUE_TYPE lb_lo, lb_hi, ub_lo, ub_hi;
+    REAL_VALUE_TYPE lb_deriv_lo, lb_deriv_hi, ub_deriv_lo, ub_deriv_hi;
+    if (!frange_mpfr_arg1 (&lb_lo, &lb_hi,
+			   m_cfn == CFN_SIN ? mpfr_sin : mpfr_cos, lb,
+			   type, ulps)
+	|| !frange_mpfr_arg1 (&ub_lo, &ub_hi,
+			      m_cfn == CFN_SIN ? mpfr_sin : mpfr_cos, ub,
+			      type, ulps)
+	|| !frange_mpfr_arg1 (&lb_deriv_lo, &lb_deriv_hi,
+			      m_cfn == CFN_SIN ? mpfr_cos : mpfr_sin, lb,
+			      type, 0)
+	|| !frange_mpfr_arg1 (&ub_deriv_lo, &ub_deriv_hi,
+			      m_cfn == CFN_SIN ? mpfr_cos : mpfr_sin, ub,
+			      type, 0))
+      return true;
+    if (m_cfn == CFN_COS)
+      {
+	// Derivative of cos is -sin, so negate.
+	lb_deriv_lo.sign ^= 1;
+	lb_deriv_hi.sign ^= 1;
+	ub_deriv_lo.sign ^= 1;
+	ub_deriv_hi.sign ^= 1;
+      }
+
+    if (real_less (&lb_lo, &ub_lo))
+      lb = lb_lo;
+    else
+      lb = ub_lo;
+    if (real_less (&lb_hi, &ub_hi))
+      ub = ub_hi;
+    else
+      ub = lb_hi;
+
+    // The range between the function result on the boundaries may need
+    // to be extended to +1 (+Inf) or -1 (-Inf) or both depending on the
+    // derivative or length of the argument range (diff).
+
+    // First handle special case, where the derivative has different signs,
+    // so the bound must be roughly -1 or +1.
+    if (real_isneg (&lb_deriv_lo) != real_isneg (&lb_deriv_hi))
+      {
+	if (real_isneg (&lb_lo))
+	  lb = dconstninf;
+	else
+	  ub = dconstinf;
+      }
+    if (real_isneg (&ub_deriv_lo) != real_isneg (&ub_deriv_hi))
+      {
+	if (real_isneg (&ub_lo))
+	  lb = dconstninf;
+	else
+	  ub = dconstinf;
+      }
+
+    // If derivative at lower_bound and upper_bound have the same sign,
+    // the function grows or declines on the whole range if diff < pi, so
+    // [lb, ub] is correct, and if diff >= pi the result range must include
+    // both the minimum and maximum.
+    if (real_isneg (&lb_deriv_lo) == real_isneg (&ub_deriv_lo))
+      {
+	if (!real_less (&diff, &pi))
+	  return true;
+      }
+    // If function declines at lower_bound and grows at upper_bound,
+    // the result range must include the minimum, so set lb to -Inf.
+    else if (real_isneg (&lb_deriv_lo))
+      lb = dconstninf;
+    // If function grows at lower_bound and declines at upper_bound,
+    // the result range must include the maximum, so set ub to +Inf.
+    else
+      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,
@@ -651,31 +743,42 @@ public:
       }
 
     // 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))
+    unsigned bulps
+      = targetm.libm_function_max_error (m_cfn, 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,
-	     [-INF,-INF][+INF,+INF] U +-NAN.  */
-	  r.set_varying (type);
-	return true;
+	const REAL_VALUE_TYPE &lb = lhs.lower_bound ();
+	const REAL_VALUE_TYPE &ub = lhs.upper_bound ();
+	REAL_VALUE_TYPE m1 = dconstm1;
+	REAL_VALUE_TYPE p1 = dconst1;
+	while (bulps--)
+	  {
+	    frange_nextafter (TYPE_MODE (type), m1, dconstninf);
+	    frange_nextafter (TYPE_MODE (type), p1, dconstinf);
+	  }
+	if (real_less (&ub, &m1) || real_less (&p1, &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);
+	REAL_VALUE_TYPE lb = real_min_representable (type);
+	REAL_VALUE_TYPE ub = real_max_representable (type);
 	r.set (type, lb, ub, nan_state (false, false));
-	return true;
       }
-
-    r.set_varying (type);
+    else
+      r.set_varying (type);
     return true;
   }
 private:
diff --git a/gcc/real.cc b/gcc/real.cc
index cf164e5..fc6b3a7 100644
--- a/gcc/real.cc
+++ b/gcc/real.cc
@@ -2475,6 +2475,26 @@ dconst_e_ptr (void)
   return &value;
 }
 
+/* Returns the special REAL_VALUE_TYPE corresponding to 'pi'.  */
+
+const REAL_VALUE_TYPE *
+dconst_pi_ptr (void)
+{
+  static REAL_VALUE_TYPE value;
+
+  /* Initialize mathematical constants for constant folding builtins.
+     These constants need to be given to at least 160 bits precision.  */
+  if (value.cl == rvc_zero)
+    {
+      auto_mpfr m (SIGNIFICAND_BITS);
+      mpfr_set_si (m, -1, MPFR_RNDN);
+      mpfr_acos (m, m, MPFR_RNDN);
+      real_from_mpfr (&value, m, NULL_TREE, MPFR_RNDN);
+
+    }
+  return &value;
+}
+
 /* Returns a cached REAL_VALUE_TYPE corresponding to 1/n, for various n.  */
 
 #define CACHED_FRACTION(NAME, N)					\
diff --git a/gcc/real.h b/gcc/real.h
index 9e02139..9ed6c37 100644
--- a/gcc/real.h
+++ b/gcc/real.h
@@ -480,9 +480,13 @@ extern REAL_VALUE_TYPE dconstninf;
 #define dconst_sixth() (*dconst_sixth_ptr ())
 #define dconst_ninth() (*dconst_ninth_ptr ())
 #define dconst_sqrt2() (*dconst_sqrt2_ptr ())
+#define dconst_pi() (*dconst_pi_ptr ())
 
 /* Function to return the real value special constant 'e'.  */
-extern const REAL_VALUE_TYPE * dconst_e_ptr (void);
+extern const REAL_VALUE_TYPE *dconst_e_ptr (void);
+
+/* Function to return the real value special constant 'pi'.  */
+extern const REAL_VALUE_TYPE *dconst_pi_ptr (void);
 
 /* Returns a cached REAL_VALUE_TYPE corresponding to 1/n, for various n.  */
 extern const REAL_VALUE_TYPE *dconst_third_ptr (void);
diff --git a/gcc/testsuite/gcc.dg/tree-ssa/range-sincos-2.c b/gcc/testsuite/gcc.dg/tree-ssa/range-sincos-2.c
new file mode 100644
index 0000000..fd578aa
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/tree-ssa/range-sincos-2.c
@@ -0,0 +1,34 @@
+// { dg-do compile }
+// { dg-options "-O2 -fdump-tree-evrp -fno-thread-jumps" }
+
+void use (double);
+void link_error ();
+
+void
+foo (double x)
+{
+  double y;
+  if (x >= 0.5 && x <= 1.3)
+    {
+      y = __builtin_sin (x);
+      if (y < 0.45 || y > 0.97)
+	link_error ();
+      use (y);
+    }
+  if (x >= 0.5 && x < 1.75)
+    {
+      y = __builtin_sin (x);
+      if (y < 0.45 || y > 1.05)
+	link_error ();
+      use (y);
+    }
+  if (x >= -1.57 && x <= 1.57)
+    {
+      y = __builtin_cos (x);
+      if (y < 0.0007 || y > 1.05)
+	link_error ();
+      use (y);
+    }
+}
+
+// { dg-final { scan-tree-dump-not "link_error" "evrp" { target { { *-*-linux* } && { glibc } } } } }