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

Hi!

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.

Bootstrapped/regtested on x86_64-linux and i686-linux, ok for trunk?

2023-05-05  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.


	Jakub

On 5/5/23 22:53, Jakub Jelinek wrote:
> Hi!
> 
> 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.
> 
> Bootstrapped/regtested on x86_64-linux and i686-linux, ok for trunk?

It's getting to the point where I think my reviews are getting less 
useful for you.  The mathematical bits are beyond my area of expertise 
and I'm just limiting myself to commenting on general style in range-op, 
etc.

I have no problem with you going ahead with this patch, but it may be 
beneficial to get someone else's opinion on the math bits.  Up to you. 
I don't want to impede your progress here.

Thanks for your work in this area.

Aldy

> 
> 2023-05-05  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.
> 
> --- gcc/real.h.jj	2023-04-19 09:33:59.434350121 +0200
> +++ gcc/real.h	2023-05-05 16:36:35.606611170 +0200
> @@ -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);
> --- gcc/real.cc.jj	2023-04-20 09:36:09.066376175 +0200
> +++ gcc/real.cc	2023-05-05 16:39:25.244201299 +0200
> @@ -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)					\
> --- gcc/gimple-range-op.cc.jj	2023-05-05 16:02:48.174419009 +0200
> +++ gcc/gimple-range-op.cc	2023-05-05 19:44:27.292304968 +0200
> @@ -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:
> --- gcc/testsuite/gcc.dg/tree-ssa/range-sincos-2.c.jj	2023-05-05 19:54:11.705064161 +0200
> +++ gcc/testsuite/gcc.dg/tree-ssa/range-sincos-2.c	2023-05-05 20:13:11.859961631 +0200
> @@ -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 } } } } }
> 
> 	Jakub
>