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.
This commit is contained in:
Jakub Jelinek 2023-05-06 10:57:41 +02:00
parent 143e6695b2
commit b7fe38c14e
4 changed files with 178 additions and 17 deletions

View file

@ -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:

View file

@ -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) \

View file

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

View file

@ -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 } } } } }