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

[Jakub Jelinek]

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
-------------- next part --------------
#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;
}