Test case: // gcc main.c -std=c99 -lm -m32 #include <stdio.h> #include <fenv.h> #include <stdlib.h> // abs function for int argument extern float y0f(float); extern double y0(double); #define i32(a) (((int*)&(a))[0]) int main() { float a,r1,r2; i32(a)=0x4155c70d; //argument fesetround(FE_TONEAREST); r1 = y0f(a); //actual result r2 = y0(a); //actual result 2 printf("inputs: a = (0x%08x)\t%e\n", i32(a), a); printf("actual = (0x%08x)\t%.13e\n", i32(r1), r1); printf("expected = (0x%08x)\t%.13e\n", i32(r2), r2); printf("error = %d ulp\n", abs(i32(r1) - i32(r2))); return 0; } Result: (core 2 for example): -bash-4.2$ gcc main_y0f.c -lm -odo -m32 -bash-4.2$ ./do inputs: a = (0x4155c70d) 1.336110e+01 actual = (0xb4801c51) -2.3862460807322e-07 expected = (0xb47ff7c8) -2.3838867946324e-07 error = 9353 ulp
And x86_64 result (core2): -bash-4.2$ gcc main_y0f.c -lm -odo -bash-4.2$ ./do inputs: a = (0x4155c70d) 1.336110e+01 actual = (0xb4800855) -2.3847920260778e-07 expected = (0xb47ff7c8) -2.3838867946324e-07 error = 4237 ulp
>> represent the test argument using a hex float value Updated test: #include <stdio.h> #include <math.h> #include <fenv.h> extern float y0f(float); extern double y0(double); int main() { float a,r1; double r2; a=0x1.ab8e1ap3; //argument fesetround(FE_TONEAREST); // FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, and FE_TOWARDZERO r1 = y0f(a); //actual result r2 = y0(a); //actual result 2 printf("inputs: a = %.6a\n", a); printf("actual = %.6a\n", r1); printf("expected = %.6a\n", r2); printf("error = %f ulp\n", fabs(r1 - r2) * exp2(23.0-logbf(r1))); return 0; } Result x86: inputs: a = 0x1.ab8e1ap+3 actual = -0x1.0038a2p-22 expected = -0x1.ffef8fp-23 error = 8301.219724 ulp Result x86_64: inputs: a = 0x1.ab8e1ap+3 actual = -0x1.0010aap-22 expected = -0x1.ffef8fp-23 error = 3185.219724 ulp
another example: Testing function y0 for exponent 0 [seed=4519]. rounding mode MPFR_RNDZ: wrong underflow flag for x=e.4c1760@-1 library gives 0 mpfr gives 3.274468@-7 underflow: mpfr 0, library 0 The correct result given by mpfr (I hope) is far away from the subnormal range (1.17472840344135e-8 vs 1.17549435082229e-38) thus the error is huge, and moreover the glibc result 0 is inconsistent with the underflow flag.
here is the worst case found by exhaustive search on all binary32 values with glibc 2.31 on x86_64: #include <stdio.h> #include <math.h> #include <fenv.h> extern float y0f(float); extern double y0(double); int main() { float a,r1; double r2; a=0x25c54.84p0; //argument fesetround(FE_TONEAREST); // FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, and FE_TOWARDZERO r1 = y0f(a); //actual result r2 = y0(a); //actual result 2 printf("inputs: a = %.6a\n", a); printf("actual = %.6a\n", r1); printf("expected = %.6a\n", r2); printf("error = %f ulp\n", fabs(r1 - r2) * exp2(23.0-logbf(r2))); return 0; } inputs: a = 0x1.2e2a42p+17 actual = -0x1.c21d2ep-30 expected = 0x1.a48974p-40 error = 15117099194.163063 ulp
with glibc-2.32 the maximum error on x86_64 is now 4837658 ulps (for x=8.935769796e-01)
fixed by commit 9acda61