(42) Double precision summation and DOT

[Albert Chan]

Double precision trick assumed round-to-nearest (so sign bit can hold a bit), halfway-to-even.

However, for decimal version, this also meant sign can only hold 1 bit.
Thus, decimal version with odd digits precision cannot be split accurately this way ...

Also, no internal extra precision (unless it were rounded with half-way-to-odd), to avoid double-rounding issues.

Code:

#include <stdio.h>

int main() // test.c
{
    double x = 1e16;
    double y = x + 2.9999;    // hit by double-rounding
    printf("%g %g %g", 2.9999+1e16-1e16, 2.9999+x-x, y-x);
    return 0;
}

Above code demonstrated problems of double-rounding.
Compile code with gcc v.7.1.0 (included with Strawberry Perl), Win7, no optimization

> gcc -otest test.c
> test
2 3 4

The reason for the discrepancy can be seen from the bits:

b = 1e16 + 2.9999 = 0b100011100001101111001001101111110000010000000000000010.1111111111111001 ...

round(b, 53-bits) = 1e16 + 2
round(b, 64-bits) = 1e16 + 3
round(round(b, 64-bits), 53-bits) = 1e16 + 4

Reference:
To SSE, or not to SSE, that is the question
GNU Linux and the Extended Precision on x86 Processors