I found something very strange, although it might be just a misunderstanding.
As far as I know, the IEEE-754 standard defines round-to-nearest, tie-to-even
as follows:
- For rounding purposes, the operation must be performed as if it were done
with infinite precision
- Then, if the bit right of the LSB of the result (guard) is 0, do nothing.
- Otherwise, if *any* bit right of the guard bit is 1, add 1 to the result
- Otherwise, if the LSB of the result is 1, add 1 to the result
- Otherwise leave the result as it is
Assuming that the above is true, and that that is the default rounding
mode, then the following is most surprising. This bit of code:
#include <stdio.h>
#include <math.h>
main()
{
double a, b;
int i;
a = 1.0;
for ( i = -54 ; i > -106 ; i-- ) {
b = ldexp( 1.0, -53 ) + ldexp( 1.0, i );
printf( "%.13a + %.13a = %.13a\n", a, b, a+b );
}
}
generates this output on an Intel i5 core:
0x1.0000000000000p+0 + 0x1.8000000000000p-53 = 0x1.0000000000001p+0
0x1.0000000000000p+0 + 0x1.4000000000000p-53 = 0x1.0000000000001p+0
0x1.0000000000000p+0 + 0x1.2000000000000p-53 = 0x1.0000000000001p+0
[...]
0x1.0000000000000p+0 + 0x1.0100000000000p-53 = 0x1.0000000000001p+0
0x1.0000000000000p+0 + 0x1.0080000000000p-53 = 0x1.0000000000001p+0
0x1.0000000000000p+0 + 0x1.0040000000000p-53 = 0x1.0000000000001p+0
0x1.0000000000000p+0 + 0x1.0020000000000p-53 = 0x1.0000000000000p+0 <== ?????
0x1.0000000000000p+0 + 0x1.0010000000000p-53 = 0x1.0000000000000p+0
0x1.0000000000000p+0 + 0x1.0008000000000p-53 = 0x1.0000000000000p+0
0x1.0000000000000p+0 + 0x1.0004000000000p-53 = 0x1.0000000000000p+0
[...]
0x1.0000000000000p+0 + 0x1.0000000000004p-53 = 0x1.0000000000000p+0
0x1.0000000000000p+0 + 0x1.0000000000002p-53 = 0x1.0000000000000p+0
0x1.0000000000000p+0 + 0x1.0000000000001p-53 = 0x1.0000000000000p+0
which seems to indicate that the sticky bit contains only the first 10 bits
right of the guard, all the rest is thrown away silently (as if internally
the operations were done on 64 bits only).
I tried to pass -fexcess-precision=standard to gcc, but got the same result.
I wonder whether
- I know the IEEE rounding rules incorrectly
- I am missing a gcc flag
- gcc is doing something funny when setting up the FPU
- Intel's FPU is not standard compiant
Thanks,
Zoltan
