diff mbox series

Implement range-op entry for sin/cos.

Message ID 20230418131250.310916-1-aldyh@redhat.com
State New
Headers show
Series Implement range-op entry for sin/cos. | expand

Commit Message

Aldy Hernandez April 18, 2023, 1:12 p.m. UTC 
[I don't know why I keep poking at floats.  I must really like the pain.
Jakub, are you OK with this patch for trunk?]

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                       | 63 ++++++++++++++++++++
 gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c | 40 +++++++++++++
 2 files changed, 103 insertions(+)
 create mode 100644 gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c

Comments

Jakub Jelinek April 20, 2023, 12:59 p.m. UTC | #1 
On Tue, Apr 18, 2023 at 03:12:50PM +0200, Aldy Hernandez wrote:
> [I don't know why I keep poking at floats.  I must really like the pain.
> Jakub, are you OK with this patch for trunk?]

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.

> 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                       | 63 ++++++++++++++++++++
>  gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c | 40 +++++++++++++
>  2 files changed, 103 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 4ca32a7b5d5..36390f2645e 100644
> --- a/gcc/gimple-range-op.cc
> +++ b/gcc/gimple-range-op.cc
> @@ -402,6 +402,60 @@ public:
>    }
>  } op_cfn_copysign;
>  
> +class cfn_sincos : public range_operator_float
> +{
> +public:
> +  using range_operator_float::fold_range;
> +  using range_operator_float::op1_range;
> +  virtual bool fold_range (frange &r, tree type,
> +			   const frange &lh, const frange &,
> +			   relation_trio) const final override
> +  {
> +    if (lh.known_isnan () || lh.known_isinf ())
> +      {
> +	r.set_nan (type);
> +	return true;
> +      }
> +    r.set (type, dconstm1, dconst1);

See above, are we sure we can use [-1., 1.] range safely, or should that be
[-1.-Nulps, 1.+Nulps] for some kind of expected worse error margin of the
implementation?  And ditto for -frounding-math, shall we increase that
interval in that case, or is [-1., 1.] going to be ok?

> +    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.

> +      r.clear_nan ();
> +    return true;

Incrementally, if we decide what to do with the maximum allowed errors in
ulps, if lh's range is smaller than 2*M_PI (upper_bound () - lower_bound
()), we could narrow it down further by computing the exact values
for the bounds and any local maximums or minimums in between if any
and creating a range out of that.

> +  }
> +  virtual bool op1_range (frange &r, tree type,
> +			  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;
> +      }
> +    // 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 (&lb, &dconstm1)
> +	|| real_less (&dconst1, &ub))

This condition could fit on one line.

> +      {
> +	r.set_undefined ();
> +	return true;
> +      }
> +
> +    r.set_varying (type);

I'm afraid we can't do really much better for the reverse case, even
for a singleton range between -1.0 and +1.0 we'd have infinite number of
results.  In theory we could perhaps narrow it down from the minimum/maximum
representable bounds, but even 1ulp there is much larger than 2*M_PI, so not
sure it is worth bothering with it.

	Jakub
Siddhesh Poyarekar April 20, 2023, 1:17 p.m. UTC | #2 
On 2023-04-20 08:59, Jakub Jelinek via Gcc-patches wrote:
>> +    r.set (type, dconstm1, dconst1);
> 
> See above, are we sure we can use [-1., 1.] range safely, or should that be
> [-1.-Nulps, 1.+Nulps] for some kind of expected worse error margin of the
> implementation?  And ditto for -frounding-math, shall we increase that
> interval in that case, or is [-1., 1.] going to be ok?

Do any math implementations generate results outside of [-1., 1.]?  If 
yes, then it's a bug in those implementations IMO, not in the range 
assumption.  It feels wrong to cater for what ought to be trivially 
fixable in libraries if they ever happen to generate such results.

Thanks,
Sid
Jakub Jelinek April 20, 2023, 2:02 p.m. UTC | #3 
On Thu, Apr 20, 2023 at 09:17:10AM -0400, Siddhesh Poyarekar wrote:
> On 2023-04-20 08:59, Jakub Jelinek via Gcc-patches wrote:
> > > +    r.set (type, dconstm1, dconst1);
> > 
> > See above, are we sure we can use [-1., 1.] range safely, or should that be
> > [-1.-Nulps, 1.+Nulps] for some kind of expected worse error margin of the
> > implementation?  And ditto for -frounding-math, shall we increase that
> > interval in that case, or is [-1., 1.] going to be ok?
> 
> Do any math implementations generate results outside of [-1., 1.]?  If yes,

Clearly they do.

> then it's a bug in those implementations IMO, not in the range assumption.
> It feels wrong to cater for what ought to be trivially fixable in libraries
> if they ever happen to generate such results.

So, I wrote following test.

On x86_64-linux with glibc 2.35, I see
for i in FLOAT DOUBLE LDOUBLE FLOAT128; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o /tmp/sincos{,.c} -lm; /tmp/sincos || echo $i $j; done; done
Aborted (core dumped)
FLOAT UPWARD
Aborted (core dumped)
FLOAT DOWNWARD
On sparc-sun-solaris2.11 I see
for i in FLOAT DOUBLE LDOUBLE; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o sincos{,.c} -lm; ./sincos || echo $i $j; done; done
Abort (core dumped)
DOUBLE UPWARD
Abort (core dumped)
DOUBLE DOWNWARD
Haven't tried anything else.  So that shows (but doesn't prove) that
maybe [-1., 1.] interval is fine for -fno-rounding-math on those, but not
for -frounding-math.

	Jakub
#define _GNU_SOURCE
#include <math.h>
#include <fenv.h>

#ifdef FLOAT
#define TYPE float
#define SIN sinf
#define COS cosf
#ifdef M_PI_2f
#define PI2 M_PI_2f
#else
#define PI2 1.570796326794896619231321691639751442f
#endif
#define NEXTAFTER nextafterf
#elif defined DOUBLE
#define TYPE double
#define SIN sin
#define COS cos
#ifdef M_PI_2
#define PI2 M_PI_2
#else
#define PI2 1.570796326794896619231321691639751442f
#endif
#define NEXTAFTER nextafter
#elif defined LDOUBLE
#define TYPE long double
#define SIN sinl
#define COS cosl
#ifdef M_PI_2l
#define PI2 M_PI_2l
#else
#define PI2 1.570796326794896619231321691639751442f
#endif
#define NEXTAFTER nextafterl
#elif defined FLOAT128
#define TYPE _Float128
#define SIN sinf128
#define COS cosf128
#ifdef 
#define PI2 M_PI_2f128
#else
#define PI2 1.570796326794896619231321691639751442f
#endif
#define NEXTAFTER nextafterf128
#endif

int
main ()
{
#ifdef ROUND
  fesetround (ROUND);
#endif
  for (int i = -1024; i <= 1024; i++)
    for (int j = -1; j <= 1; j += 2)
      {
        TYPE val = ((TYPE) i) * PI2;
        TYPE inf = j * __builtin_inf ();
        for (int k = 0; k < 1000; k++)
          {
            TYPE res = SIN (val);
            if (res < (TYPE) -1.0 || res > (TYPE) 1.0)
	      __builtin_abort ();
	    res = COS (val);
            if (res < (TYPE) -1.0 || res > (TYPE) 1.0)
	      __builtin_abort ();
	    val = NEXTAFTER (val, inf);
          }
      }
}
Jakub Jelinek April 20, 2023, 2:20 p.m. UTC | #4 
On Thu, Apr 20, 2023 at 04:02:02PM +0200, Jakub Jelinek via Gcc-patches wrote:
> So, I wrote following test.

Slightly adjusted to see more info:

x86_64-linux glibc 2.35:
for i in FLOAT DOUBLE LDOUBLE FLOAT128; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o /tmp/sincos{,.c} -lm; /tmp/sincos || echo $i $j; done; done
sin -0x1.2d97c800000000000000p+2 0x1.00000200000000000000p+0
sin -0x1.2d97c800000000000000p+2 0x1.00000200000000000000p+0
sin 0x1.f9cbe200000000000000p+7 0x1.00000200000000000000p+0
FLOAT UPWARD
cos -0x1.f9cbe200000000000000p+8 -0x1.00000200000000000000p+0
sin -0x1.f9cbe200000000000000p+7 -0x1.00000200000000000000p+0
sin 0x1.2d97c800000000000000p+2 -0x1.00000200000000000000p+0
sin 0x1.2d97c800000000000000p+2 -0x1.00000200000000000000p+0
cos 0x1.f9cbe200000000000000p+8 -0x1.00000200000000000000p+0
cos 0x1.f9cbe200000000000000p+8 -0x1.00000200000000000000p+0
FLOAT DOWNWARD
sparc-sun-solaris2.11 results are too large to post in detail, but are
sin -0x1.2d97c7f3321d20000000p+2 0x1.00000000000010000000p+0
...
sin 0x1.f6a7a29553c450000000p+2 0x1.00000000000010000000p+0
DOUBLE UPWARD
sin -0x1.f6a7a2955385e0000000p+2 -0x1.00000000000010000000p+0
...
sin 0x1.2d97c7f3325b90000000p+2 -0x1.00000000000010000000p+0
DOUBLE DOWNWARD
where all the DOUBLE UPWARD values have 0x1.00000000000010000000p+0
results and all DOUBLE DOWNWARD values -0x1.00000000000010000000p+0.
So, I think that is 1ulp in all cases in both directions for
-frounding-math.

	Jakub
#define _GNU_SOURCE
#include <math.h>
#include <fenv.h>
#include <stdio.h>
#include <stdlib.h>

#ifdef FLOAT
#define TYPE float
#define SIN sinf
#define COS cosf
#ifdef M_PI_2f
#define PI2 M_PI_2f
#else
#define PI2 1.570796326794896619231321691639751442f
#endif
#define PRINT(str) printf ("%s %.20a %.20a\n", str, val, res)
#define NEXTAFTER nextafterf
#elif defined DOUBLE
#define TYPE double
#define SIN sin
#define COS cos
#ifdef M_PI_2
#define PI2 M_PI_2
#else
#define PI2 1.570796326794896619231321691639751442f
#endif
#define NEXTAFTER nextafter
#define PRINT(str) printf ("%s %.20a %.20a\n", str, val, res)
#elif defined LDOUBLE
#define TYPE long double
#define SIN sinl
#define COS cosl
#ifdef M_PI_2l
#define PI2 M_PI_2l
#else
#define PI2 1.570796326794896619231321691639751442f
#endif
#define NEXTAFTER nextafterl
#define PRINT(str) printf ("%s %.20La %.20La\n", str, val, res)
#elif defined FLOAT128
#define TYPE _Float128
#define SIN sinf128
#define COS cosf128
#ifdef M_PI_2f128
#define PI2 M_PI_2f128
#else
#define PI2 1.570796326794896619231321691639751442f
#endif
#define NEXTAFTER nextafterf128
#define PRINT(str) __builtin_abort ()
#endif

int
main ()
{
  int ret = 0;
#ifdef ROUND
  fesetround (ROUND);
#endif
  for (int i = -1024; i <= 1024; i++)
    for (int j = -1; j <= 1; j += 2)
      {
        TYPE val = ((TYPE) i) * PI2;
        TYPE inf = j * __builtin_inf ();
        for (int k = 0; k < 1000; k++)
          {
            TYPE res = SIN (val);
            if (res < (TYPE) -1.0 || res > (TYPE) 1.0)
              { PRINT ("sin"); ret = 1; }
	    res = COS (val);
            if (res < (TYPE) -1.0 || res > (TYPE) 1.0)
	      { PRINT ("cos"); ret = 1; }
	    val = NEXTAFTER (val, inf);
          }
      }
  return ret;
}
Siddhesh Poyarekar April 20, 2023, 3:22 p.m. UTC | #5 
On 2023-04-20 10:02, Jakub Jelinek wrote:
> On x86_64-linux with glibc 2.35, I see
> for i in FLOAT DOUBLE LDOUBLE FLOAT128; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o /tmp/sincos{,.c} -lm; /tmp/sincos || echo $i $j; done; done
> Aborted (core dumped)
> FLOAT UPWARD
> Aborted (core dumped)
> FLOAT DOWNWARD
> On sparc-sun-solaris2.11 I see
> for i in FLOAT DOUBLE LDOUBLE; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o sincos{,.c} -lm; ./sincos || echo $i $j; done; done
> Abort (core dumped)
> DOUBLE UPWARD
> Abort (core dumped)
> DOUBLE DOWNWARD
> Haven't tried anything else.  So that shows (but doesn't prove) that
> maybe [-1., 1.] interval is fine for -fno-rounding-math on those, but not
> for -frounding-math.

Would there be a reason to not consider these as bugs?  I feel like 
these should be fixed in glibc, or any math implementation that ends up 
doing this.

I suppose one reason could be the overhead of an additional branch to 
check for result bounds, but is that serious enough to allow this 
imprecision?  The alternative of output range being defined as 
[-1.0-ulp, 1.0+ulp] avoids that conversation I guess.

Thanks,
Sid
Jakub Jelinek April 20, 2023, 3:52 p.m. UTC | #6 
On Thu, Apr 20, 2023 at 11:22:24AM -0400, Siddhesh Poyarekar wrote:
> On 2023-04-20 10:02, Jakub Jelinek wrote:
> > On x86_64-linux with glibc 2.35, I see
> > for i in FLOAT DOUBLE LDOUBLE FLOAT128; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o /tmp/sincos{,.c} -lm; /tmp/sincos || echo $i $j; done; done
> > Aborted (core dumped)
> > FLOAT UPWARD
> > Aborted (core dumped)
> > FLOAT DOWNWARD
> > On sparc-sun-solaris2.11 I see
> > for i in FLOAT DOUBLE LDOUBLE; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o sincos{,.c} -lm; ./sincos || echo $i $j; done; done
> > Abort (core dumped)
> > DOUBLE UPWARD
> > Abort (core dumped)
> > DOUBLE DOWNWARD
> > Haven't tried anything else.  So that shows (but doesn't prove) that
> > maybe [-1., 1.] interval is fine for -fno-rounding-math on those, but not
> > for -frounding-math.
> 
> Would there be a reason to not consider these as bugs?  I feel like these
> should be fixed in glibc, or any math implementation that ends up doing
> this.

Why?  Unless an implementation guarantees <= 0.5ulps errors, it can be one
or more ulps off, why is an error at or near 1.0 or -1.0 error any worse
than similar errors for other values?
Similarly for other functions which have other ranges, perhaps not with so
nice round numbers.  Say asin has [-pi/2, pi/2] range, those numbers aren't
exactly representable, but is it any worse to round those values to -inf or
+inf or worse give something 1-5 ulps further from that interval comparing
to other 1-5ulps errors?

	Jakub
Siddhesh Poyarekar April 20, 2023, 5:57 p.m. UTC | #7 
On 2023-04-20 11:52, Jakub Jelinek wrote:
> Why?  Unless an implementation guarantees <= 0.5ulps errors, it can be one
> or more ulps off, why is an error at or near 1.0 or -1.0 error any worse
> than similar errors for other values?

In a general sense, maybe not, but in the sense of breaching the bounds 
of admissible values, especially when it can be reasonably corrected, it 
seems worse IMO to let the error slide.

> Similarly for other functions which have other ranges, perhaps not with so
> nice round numbers.  Say asin has [-pi/2, pi/2] range, those numbers aren't
> exactly representable, but is it any worse to round those values to -inf or
> +inf or worse give something 1-5 ulps further from that interval comparing
> to other 1-5ulps errors?

I agree the argument in favour of allowing errors breaching the 
mathematical bounds is certainly stronger for bounds that are not 
exactly representable.  I just feel like the implementation ought to 
take the additional effort when the bounds are representable and make it 
easier for the compiler.

For bounds that aren't representable, one could get error bounds from 
libm-test-ulps data in glibc, although I reckon those won't be 
exhaustive.  From a quick peek at the sin/cos data, the arc target seems 
to be among the worst performers at about 7ulps, although if you include 
the complex routines we get close to 13 ulps.  The very worst 
imprecision among all math routines (that's gamma) is at 16 ulps for 
power in glibc tests, so maybe allowing about 25-30 ulps error in bounds 
might work across the board.

But yeah, it's probably going to be guesswork.

Thanks,
Sid
Siddhesh Poyarekar April 21, 2023, 1:14 a.m. UTC | #8 
On 2023-04-20 13:57, Siddhesh Poyarekar wrote:
> For bounds that aren't representable, one could get error bounds from 
> libm-test-ulps data in glibc, although I reckon those won't be 
> exhaustive.  From a quick peek at the sin/cos data, the arc target seems 
> to be among the worst performers at about 7ulps, although if you include 
> the complex routines we get close to 13 ulps.  The very worst 
> imprecision among all math routines (that's gamma) is at 16 ulps for 
> power in glibc tests, so maybe allowing about 25-30 ulps error in bounds 
> might work across the board.

I was thinking about this a bit more and it seems like limiting ranges 
to targets that can generate sane results (i.e. error bounds within, 
say, 5-6 ulps) and for the rest, avoid emitting the ranges altogether. 
Emitting a bad range for all architectures seems like a net worse 
solution again.

Thanks,
Sid
Jakub Jelinek April 21, 2023, 6:52 a.m. UTC | #9 
On Thu, Apr 20, 2023 at 09:14:10PM -0400, Siddhesh Poyarekar wrote:
> On 2023-04-20 13:57, Siddhesh Poyarekar wrote:
> > For bounds that aren't representable, one could get error bounds from
> > libm-test-ulps data in glibc, although I reckon those won't be
> > exhaustive.  From a quick peek at the sin/cos data, the arc target seems
> > to be among the worst performers at about 7ulps, although if you include
> > the complex routines we get close to 13 ulps.  The very worst
> > imprecision among all math routines (that's gamma) is at 16 ulps for
> > power in glibc tests, so maybe allowing about 25-30 ulps error in bounds
> > might work across the board.
> 
> I was thinking about this a bit more and it seems like limiting ranges to
> targets that can generate sane results (i.e. error bounds within, say, 5-6
> ulps) and for the rest, avoid emitting the ranges altogether. Emitting a bad
> range for all architectures seems like a net worse solution again.

Well, at least for basic arithmetics when libm functions aren't involved,
there is no point in disabling ranges altogether.
And, for libm functions, my plan was to introduce a target hook, which
would have combined_fn argument to tell which function is queried,
machine_mode to say which floating point format and perhaps a bool whether
it is ulps for these basic math boundaries or results somewhere in between,
and would return in unsigned int ulps, 0 for 0.5ulps precision.
So, we could say for CASE_CFN_SIN: CASE_CFN_COS: in the glibc handler
say that ulps is say 3 inside of the ranges and 0 on the boundaries if
!flag_rounding_math and 6 and 2 with flag_rounding_math or whatever.
And in the generic code don't assume anything if ulps is say 100 or more.
The hooks would need to be a union of precision of supported versions of
the library through the history, including say libmvec because function
calls could be vectorized.
And default could be that infinite precision.
Back in November I've posted a proglet that can generate some ulps from
random number testing, plus on glibc we could pick maximums from ulps files.
And if needed, say powerpc*-linux could override the generic glibc
version for some subset of functions and call default otherwise (say at
least for __ibm128).

	Jakub
Siddhesh Poyarekar April 21, 2023, 11:20 a.m. UTC | #10 
On 2023-04-21 02:52, Jakub Jelinek wrote:
> On Thu, Apr 20, 2023 at 09:14:10PM -0400, Siddhesh Poyarekar wrote:
>> On 2023-04-20 13:57, Siddhesh Poyarekar wrote:
>>> For bounds that aren't representable, one could get error bounds from
>>> libm-test-ulps data in glibc, although I reckon those won't be
>>> exhaustive.  From a quick peek at the sin/cos data, the arc target seems
>>> to be among the worst performers at about 7ulps, although if you include
>>> the complex routines we get close to 13 ulps.  The very worst
>>> imprecision among all math routines (that's gamma) is at 16 ulps for
>>> power in glibc tests, so maybe allowing about 25-30 ulps error in bounds
>>> might work across the board.
>>
>> I was thinking about this a bit more and it seems like limiting ranges to
>> targets that can generate sane results (i.e. error bounds within, say, 5-6
>> ulps) and for the rest, avoid emitting the ranges altogether. Emitting a bad
>> range for all architectures seems like a net worse solution again.
> 
> Well, at least for basic arithmetics when libm functions aren't involved,
> there is no point in disabling ranges altogether.

Oh yeah, I did mean only franges for math function call results.

> And, for libm functions, my plan was to introduce a target hook, which
> would have combined_fn argument to tell which function is queried,
> machine_mode to say which floating point format and perhaps a bool whether
> it is ulps for these basic math boundaries or results somewhere in between,
> and would return in unsigned int ulps, 0 for 0.5ulps precision.
> So, we could say for CASE_CFN_SIN: CASE_CFN_COS: in the glibc handler
> say that ulps is say 3 inside of the ranges and 0 on the boundaries if
> !flag_rounding_math and 6 and 2 with flag_rounding_math or whatever.
> And in the generic code don't assume anything if ulps is say 100 or more.
> The hooks would need to be a union of precision of supported versions of
> the library through the history, including say libmvec because function
> calls could be vectorized.
> And default could be that infinite precision.
> Back in November I've posted a proglet that can generate some ulps from
> random number testing, plus on glibc we could pick maximums from ulps files.
> And if needed, say powerpc*-linux could override the generic glibc
> version for some subset of functions and call default otherwise (say at
> least for __ibm128).

Ack, that sounds like a plan.

Thanks,
Sid
Jakub Jelinek April 21, 2023, 4:40 p.m. UTC | #11 
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
Mikael Morin April 21, 2023, 8:43 p.m. UTC | #12 
Hello,

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

Shouldn't lb and ub be swapped in this condition?
If I understand correctly, we are looking for ranges like [whatever,x] 
where x < -1.0 or [y, whatever] where 1.0 < y.

Mikael
Jakub Jelinek April 21, 2023, 8:45 p.m. UTC | #13 
On Fri, Apr 21, 2023 at 10:43:44PM +0200, Mikael Morin wrote:
> Hello,
> 
> > --- 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
> > @@ -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))
> >         {
> 
> Shouldn't lb and ub be swapped in this condition?

Yes, they should.

> If I understand correctly, we are looking for ranges like [whatever,x] where
> x < -1.0 or [y, whatever] where 1.0 < y.

	Jakub
Jakub Jelinek April 24, 2023, 4:03 p.m. UTC | #14 
On Thu, Apr 20, 2023 at 01:57:55PM -0400, Siddhesh Poyarekar wrote:
> > Similarly for other functions which have other ranges, perhaps not with so
> > nice round numbers.  Say asin has [-pi/2, pi/2] range, those numbers aren't
> > exactly representable, but is it any worse to round those values to -inf or
> > +inf or worse give something 1-5 ulps further from that interval comparing
> > to other 1-5ulps errors?

I've extended my hack^^^test to include sqrt and this time it seems
that the boundary in that case holds even for non-standard rounding modes
for glibc:
for i in FLOAT DOUBLE LDOUBLE FLOAT128; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o /tmp/sincossqrt{,.c} -lm; /tmp/sincossqrt || echo $i $j; done; done
sin -0x1.2d97c800000000000000p+2 0x1.00000200000000000000p+0
sin -0x1.2d97c800000000000000p+2 0x1.00000200000000000000p+0
sin 0x1.f9cbe200000000000000p+7 0x1.00000200000000000000p+0
FLOAT UPWARD
cos -0x1.f9cbe200000000000000p+8 -0x1.00000200000000000000p+0
sin -0x1.f9cbe200000000000000p+7 -0x1.00000200000000000000p+0
sin 0x1.2d97c800000000000000p+2 -0x1.00000200000000000000p+0
sin 0x1.2d97c800000000000000p+2 -0x1.00000200000000000000p+0
cos 0x1.f9cbe200000000000000p+8 -0x1.00000200000000000000p+0
cos 0x1.f9cbe200000000000000p+8 -0x1.00000200000000000000p+0
FLOAT DOWNWARD

	Jakub
#define _GNU_SOURCE
#include <math.h>
#include <fenv.h>
#include <stdio.h>
#include <stdlib.h>

#ifdef FLOAT
#define TYPE float
#define SIN sinf
#define COS cosf
#define SQRT sqrtf
#ifdef M_PI_2f
#define PI2 M_PI_2f
#else
#define PI2 1.570796326794896619231321691639751442f
#endif
#define PRINT(str) printf ("%s %.20a %.20a\n", str, val, res)
#define NEXTAFTER nextafterf
#elif defined DOUBLE
#define TYPE double
#define SIN sin
#define COS cos
#define SQRT sqrt
#ifdef M_PI_2
#define PI2 M_PI_2
#else
#define PI2 1.570796326794896619231321691639751442f
#endif
#define NEXTAFTER nextafter
#define PRINT(str) printf ("%s %.20a %.20a\n", str, val, res)
#elif defined LDOUBLE
#define TYPE long double
#define SIN sinl
#define COS cosl
#define SQRT sqrtl
#ifdef M_PI_2l
#define PI2 M_PI_2l
#else
#define PI2 1.570796326794896619231321691639751442f
#endif
#define NEXTAFTER nextafterl
#define PRINT(str) printf ("%s %.20La %.20La\n", str, val, res)
#elif defined FLOAT128
#define TYPE _Float128
#define SIN sinf128
#define COS cosf128
#define SQRT sqrtf128
#ifdef M_PI_2f128
#define PI2 M_PI_2f128
#else
#define PI2 1.570796326794896619231321691639751442f
#endif
#define NEXTAFTER nextafterf128
#define PRINT(str) __builtin_abort ()
#endif

int
main ()
{
  int ret = 0;
#ifdef ROUND
  fesetround (ROUND);
#endif
  for (int i = -1024; i <= 1024; i++)
    for (int j = -1; j <= 1; j += 2)
      {
        TYPE val = ((TYPE) i) * PI2;
        TYPE inf = j * __builtin_inf ();
        for (int k = 0; k < 1000; k++)
          {
            TYPE res = SIN (val);
            if (res < (TYPE) -1.0 || res > (TYPE) 1.0)
              { PRINT ("sin"); ret = 1; }
	    res = COS (val);
            if (res < (TYPE) -1.0 || res > (TYPE) 1.0)
	      { PRINT ("cos"); ret = 1; }
	    val = NEXTAFTER (val, inf);
          }
      }
  for (int j = -1; j <= 1; j += 2)
    {
      TYPE val = 0;
      TYPE inf = j * __builtin_inf ();
      if (j < 0)
	val = -val;
      for (int k = 0; k < 1000; k++)
        {
          TYPE res = SQRT (val);
          if (res < (TYPE) -0.0)
	    { PRINT ("sqrt"); ret = 1; }
        }
    }
  return ret;
}
Siddhesh Poyarekar April 24, 2023, 4:05 p.m. UTC | #15 
On 2023-04-24 12:03, Jakub Jelinek wrote:
> On Thu, Apr 20, 2023 at 01:57:55PM -0400, Siddhesh Poyarekar wrote:
>>> Similarly for other functions which have other ranges, perhaps not with so
>>> nice round numbers.  Say asin has [-pi/2, pi/2] range, those numbers aren't
>>> exactly representable, but is it any worse to round those values to -inf or
>>> +inf or worse give something 1-5 ulps further from that interval comparing
>>> to other 1-5ulps errors?
> 
> I've extended my hack^^^test to include sqrt and this time it seems
> that the boundary in that case holds even for non-standard rounding modes
> for glibc:

IIRC the standard _requires_ sqrt to be correctly rounded.

Sid
Jakub Jelinek April 24, 2023, 4:09 p.m. UTC | #16 
On Mon, Apr 24, 2023 at 12:05:43PM -0400, Siddhesh Poyarekar wrote:
> On 2023-04-24 12:03, Jakub Jelinek wrote:
> > On Thu, Apr 20, 2023 at 01:57:55PM -0400, Siddhesh Poyarekar wrote:
> > > > Similarly for other functions which have other ranges, perhaps not with so
> > > > nice round numbers.  Say asin has [-pi/2, pi/2] range, those numbers aren't
> > > > exactly representable, but is it any worse to round those values to -inf or
> > > > +inf or worse give something 1-5 ulps further from that interval comparing
> > > > to other 1-5ulps errors?
> > 
> > I've extended my hack^^^test to include sqrt and this time it seems
> > that the boundary in that case holds even for non-standard rounding modes
> > for glibc:
> 
> IIRC the standard _requires_ sqrt to be correctly rounded.

Well, we still need to make some effort to verify it is the case.

BTW, what my hacks didn't check yet and I should eventually is libmvec,
because if ranger makes some assumptions on libm calls and then we vectorize
it, it will be libmvec rather than libm.

	Jakub
Jeff Law April 24, 2023, 4:33 p.m. UTC | #17 
On 4/24/23 10:05, Siddhesh Poyarekar wrote:
> On 2023-04-24 12:03, Jakub Jelinek wrote:
>> On Thu, Apr 20, 2023 at 01:57:55PM -0400, Siddhesh Poyarekar wrote:
>>>> Similarly for other functions which have other ranges, perhaps not 
>>>> with so
>>>> nice round numbers.  Say asin has [-pi/2, pi/2] range, those numbers 
>>>> aren't
>>>> exactly representable, but is it any worse to round those values to 
>>>> -inf or
>>>> +inf or worse give something 1-5 ulps further from that interval 
>>>> comparing
>>>> to other 1-5ulps errors?
>>
>> I've extended my hack^^^test to include sqrt and this time it seems
>> that the boundary in that case holds even for non-standard rounding modes
>> for glibc:
> 
> IIRC the standard _requires_ sqrt to be correctly rounded.
Correct.  I spent an inordinate amount of time dealing with this in the 
past ;-)

jeff
Aldy Hernandez April 25, 2023, 8:59 a.m. UTC | #18 
On 4/21/23 08:52, Jakub Jelinek wrote:
> On Thu, Apr 20, 2023 at 09:14:10PM -0400, Siddhesh Poyarekar wrote:
>> On 2023-04-20 13:57, Siddhesh Poyarekar wrote:
>>> For bounds that aren't representable, one could get error bounds from
>>> libm-test-ulps data in glibc, although I reckon those won't be
>>> exhaustive.  From a quick peek at the sin/cos data, the arc target seems
>>> to be among the worst performers at about 7ulps, although if you include
>>> the complex routines we get close to 13 ulps.  The very worst
>>> imprecision among all math routines (that's gamma) is at 16 ulps for
>>> power in glibc tests, so maybe allowing about 25-30 ulps error in bounds
>>> might work across the board.
>>
>> I was thinking about this a bit more and it seems like limiting ranges to
>> targets that can generate sane results (i.e. error bounds within, say, 5-6
>> ulps) and for the rest, avoid emitting the ranges altogether. Emitting a bad
>> range for all architectures seems like a net worse solution again.
> 
> Well, at least for basic arithmetics when libm functions aren't involved,
> there is no point in disabling ranges altogether.
> And, for libm functions, my plan was to introduce a target hook, which
> would have combined_fn argument to tell which function is queried,
> machine_mode to say which floating point format and perhaps a bool whether
> it is ulps for these basic math boundaries or results somewhere in between,
> and would return in unsigned int ulps, 0 for 0.5ulps precision.

I wonder if we could export frange_nextafter() to take a final argument 
for the number of ulps, making it possible to extend in either direction 
by a number of ulps.  And perhaps add a generic function that calls the 
target hook and extends the range accordingly:

extern int TARGET_ULP_ERROR (combined_fn, machine_mode);
extern void frange_nextafter (machine_mode,
			      REAL_VALUE_TYPE &value,
			      const REAL_VALUE_TYPE &direction,
			      int ulps);
void
adjust_frange_for_op (frange &r, combined_fn fn)
{
   ...
   int ulps = TARGET_ULP_ERROR (fn, mode);
   frange_nextafter (mode, lb, frange_val_min (type), ulps);
   frange_nextafter (mode, ub, frange_val_max (type), ulps);
   ...
}

bool
cfn_sincos::fold_range (...)
{
...
  if (cond)
     {
       r.set (type, dconstm1, const1);
       adjust_frange_for_op (r, m_cfn);
     }
...
}

Or if adjusting by ULPS is a common enough idiom, we could promote it to 
an frange method:

	void frange::adjust_ulps (int ulps);

Would that work, or did you have something else in mind?

Aldy
Aldy Hernandez April 25, 2023, 9:08 a.m. UTC | #19 
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
Aldy Hernandez April 25, 2023, 9:10 a.m. UTC | #20 
On 4/21/23 22:43, Mikael Morin wrote:
> Hello,
> 
>> --- 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
>> @@ -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))
>>         {
> 
> Shouldn't lb and ub be swapped in this condition?
> If I understand correctly, we are looking for ranges like [whatever,x] 
> where x < -1.0 or [y, whatever] where 1.0 < y.
> 
> Mikael
> 

Thanks.  Merged into the latest revision.
Aldy
diff mbox series

Patch

diff --git a/gcc/gimple-range-op.cc b/gcc/gimple-range-op.cc
index 4ca32a7b5d5..36390f2645e 100644
--- a/gcc/gimple-range-op.cc
+++ b/gcc/gimple-range-op.cc
@@ -402,6 +402,60 @@  public:
   }
 } op_cfn_copysign;
 
+class cfn_sincos : public range_operator_float
+{
+public:
+  using range_operator_float::fold_range;
+  using range_operator_float::op1_range;
+  virtual bool fold_range (frange &r, tree type,
+			   const frange &lh, const frange &,
+			   relation_trio) const final override
+  {
+    if (lh.known_isnan () || lh.known_isinf ())
+      {
+	r.set_nan (type);
+	return true;
+      }
+    r.set (type, dconstm1, dconst1);
+    if (!lh.maybe_isnan ())
+      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.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;
+      }
+    // 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 (&lb, &dconstm1)
+	|| real_less (&dconst1, &ub))
+      {
+	r.set_undefined ();
+	return true;
+      }
+
+    r.set_varying (type);
+    return true;
+  }
+} op_cfn_sincos;
+
 // Implement range operator for CFN_BUILT_IN_TOUPPER and CFN_BUILT_IN_TOLOWER.
 class cfn_toupper_tolower : public range_operator
 {
@@ -880,6 +934,15 @@  gimple_range_op_handler::maybe_builtin_call ()
       m_valid = true;
       break;
 
+    CASE_CFN_SIN:
+    CASE_CFN_SIN_FN:
+    CASE_CFN_COS:
+    CASE_CFN_COS_FN:
+      m_op1 = gimple_call_arg (call, 0);
+      m_float = &op_cfn_sincos;
+      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" } }