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 j1f(float); extern double j1(double); #define i32(a) (((int*)&(a))[0]) int main() { float a,r1,r2; i32(a)=0x4122c68a; //argument fesetround(FE_TONEAREST); r1 = j1f(a); //actual result r2 = j1(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; } Results: inputs: a = (0x4122c68a) 1.017347e+01 actual = (0xb55e2ad3) -8.2763762065952e-07 expected = (0xb55e40c6) -8.2795702383009e-07 error = 5619 ulp
In bug description result 5619 ulps was the result for x86_64 and for 32 bit mode result is: -bash-4.2$ gcc main_j1f.c -lm -odo -m32 -bash-4.2$ ./do inputs: a = (0x4122c68a) 1.017347e+01 actual = (0xb55e4d6d) -8.2814113966378e-07 expected = (0xb55e40c6) -8.2795702383009e-07 error = 3239 ulp
represent the test argument using a hex float value, test case update: #include <stdio.h> #include <math.h> #include <fenv.h> extern float j1f(float); extern double j1(double); int main() { float a,r1; double r2; a=0x1.458d14p3; //argument fesetround(FE_TONEAREST); r1 = j1f(a); //actual result r2 = j1(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.458d14p+3 actual = -0x1.bc9adap-21 expected = -0x1.bc818dp-21 error = 3238.561446 ulp Result x86_64: inputs: a = 0x1.458d14p+3 actual = -0x1.bc55a6p-21 expected = -0x1.bc818dp-21 error = 5619.438563 ulp
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 j1f(float); extern double j1(double); int main() { float a,r1; double r2; a=0x3ba13.e8p0; //argument fesetround(FE_TONEAREST); r1 = j1f(a); //actual result r2 = j1(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.dd09f4p+17 actual = 0x1.4d4036p-29 expected = -0x1.db6a7ep-35 error = 714456063.171884 ulp
with glibc-2.32, the maximum error on x86_64 is now 2246675 ulps (for x=1.892188179e+38)
fixed by commit 9acda61