Re: [PATCH] Implement range-op entry for sin/cos.

[Jakub Jelinek <jakub@redhat.com>]

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.

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


	Jakub