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

[Jakub Jelinek via Gcc-patches]

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