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