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