On 4/21/23 18:40, Jakub Jelinek wrote:
On Thu, Apr 20, 2023 at 02:59:35PM +0200, Jakub Jelinek via Gcc-patches wrote:
Thanks for working on this.  Though expectedly here we are running
again into the discussions we had in November about math properties of the
functions vs. numeric properties in their implementations, how big maximum
error shall we expect for the functions (and whether we should hardcode
it for all implementations, or have some more fine-grained list of expected
ulp errors for each implementation), whether the implementations at least
guarantee the basic mathematical properties of the functions even if they
have some errors (say for sin/cos, if they really never return > 1.0 or <
-1.0) and the same questions repeated for -frounding-math, what kind of
extra errors to expect when using non-default rounding and whether say sin
could ever return nextafter (1.0, 2.0) or even larger value say when
using non-default rounding mode.
https://gcc.gnu.org/pipermail/gcc-patches/2022-November/606466.html
was my attempt to get at least some numbers on some targets, I'm afraid
for most implementations we aren't going to get any numerical proofs of
maximum errors and the like.  For sin/cos to check whether the implementation
really never returns > 1.0 or < -1.0 perhaps instead of using randomized
testing we could exhaustively check some range around 0, M_PI, 3*M_PI,
-M_PI, -3*M_PI, and say some much larger multiples of M_PI, say 50 ulps
in each direction about those points, and similarly for sin around M_PI/2
etc., in all arounding modes.

Appart from the ulps ranges which I plan to introduce a target hook next
week...

+    if (!lh.maybe_isnan ())

This condition isn't sufficient, even if lh can't be NAN, but just
may be +Inf or -Inf, the result needs to include maybe NAN.

I've incorporated my other comments into a patch form.
I also found missing undefined_p () checks in two spots, the
r.set_undefined (); stuff being misplaced (it was done in the
lhs.maybe_isnan () case, which is incorrect, if the lhs may be nan,
then even if the finite range is say [-30., -10.] or [1.5., 42.],
result shouldn't be invalid, because the result still could be NAN
and in that case operand of say +-Inf or NAN would be valid.
Actually thinking about it some more, perhaps we should do that check
before the maybe_isnan () check and if we find the impossible finite,
either use r.set_undefined (); of !maybe_isnan () or handle it like
known_isnan () otherwise.

Also, I think we should remember if it is SIN or COS, we'll need it both
for the ulps case and if we improve it for (ub-lb) < 2*PI ranges.
Now, talking about that, I'd very much like to avoid finding if some
multiple of PI/2 occurs inside of such ranges, the precision of real.cc
is clearly not sufficient for that, but perhaps we could use derivative
of sin (cos) or of cos (sin) to see if the function on the boundary is
increasing or decreasing and from that on both boundaries and their
approximate distance find out if the range needs to be extended to +1
or -1.

Could we do this PI stuff as a follow-up, or should it go in this patch?


So, just incremental WIP so far...

--- gcc/gimple-range-op.cc.jj   2023-04-21 17:09:48.250367999 +0200
+++ gcc/gimple-range-op.cc      2023-04-21 18:37:26.048325391 +0200
@@ -405,17 +405,20 @@ class cfn_sincos : public range_operator
  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;
        }
      r.set (type, dconstm1, dconst1);
-    if (!lh.maybe_isnan ())
+    if (!lh.maybe_isnan () && !lh.maybe_isinf ())
        r.clear_nan ();
      return true;
    }
@@ -423,15 +426,9 @@ public:
                          const frange &lhs, const frange &,
                          relation_trio) const final override
    {
-    if (!lhs.maybe_isnan ())
-      {
-       // If NAN is not valid result, the input cannot include either
-       // a NAN nor a +-INF.
-       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;
-      }
+    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 ())
@@ -439,20 +436,38 @@ public:
        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 (&lb, &dconstm1)
-       || real_less (&dconst1, &ub))
+    if (real_less (&lb, &dconstm1) || real_less (&dconst1, &ub))
        {
-       r.set_undefined ();
+       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));

I don't like the nan_state(false, false) idiom, and that's my fault. For that matter, the nan_state() default constructor is confusing because it's not clear whether that means a NAN or not a NAN. I'm fixing this in a patch I just posted.

        return true;
        }
r.set_varying (type);
      return true;
    }
-} op_cfn_sincos;
+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
@@ -934,10 +949,15 @@ gimple_range_op_handler::maybe_builtin_c
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_sincos;
+      m_float = &op_cfn_cos;
        m_valid = true;
        break;

Thanks.

I've merged your work with mine and am retesting. I have also swapped the lb/ub per Mikael's comment downthread.

Should we adjust the range by say 50 ulps in either direction, and do your target hook as a follow-up?

Aldy
From 065b8a1c1ec582eb3445ddb7f87ce9adc74394db Mon Sep 17 00:00:00 2001
From: Aldy Hernandez <al...@redhat.com>
Date: Wed, 29 Mar 2023 12:46:38 +0200
Subject: [PATCH] Implement range-op entry for sin/cos.

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.

gcc/ChangeLog:

	* gimple-range-op.cc (class cfn_sincos): New.
	(gimple_range_op_handler::maybe_builtin_call): Add case for sin/cos.

gcc/testsuite/ChangeLog:

	* gcc.dg/tree-ssa/range-sincos.c: New test.
---
 gcc/gimple-range-op.cc                       | 83 ++++++++++++++++++++
 gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c | 40 ++++++++++
 2 files changed, 123 insertions(+)
 create mode 100644 gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c

diff --git a/gcc/gimple-range-op.cc b/gcc/gimple-range-op.cc
index f7409e35a99..90305b8614a 100644
--- a/gcc/gimple-range-op.cc
+++ b/gcc/gimple-range-op.cc
@@ -400,6 +400,75 @@ 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;
+      }
+    r.set (type, dconstm1, dconst1);
+    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 a 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));
+	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
 {
@@ -878,6 +947,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 00000000000..50fbac1b68f
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c
@@ -0,0 +1,40 @@
+// { 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" } }
-- 
2.40.0

Reply via email to