Double precision strange behaviour. Need an explanation

Question

Here's the code:

#include <stdio.h>
#include <math.h>

static double const x = 665857;
static double const y = 470832;

int main(){
    double z = x*x*x*x -(y*y*y*y*4+y*y*4);
    printf("%f \n",z);
    return 0;
}

Mysteriously (to me) this code prints "0.0" if compiled on 32 bits machines (or with the -m32 flag on 64 bits machines like in my case) with GCC 4.6. As far as I know about floating point operations, it is possible to overflow/underflow them or to lose precision with them, but... a 0? How?

Thanks in advance.

Answer 1

The problem is not that the numbers overflow. The problem is that doubles don't have enough precision to distinguish between the two operands of your subtraction.

The value of x*x*x*x is 196573006004558194713601.

The value of y*y*y*y*4+y*y*4 is 196573006004558194713600.

These numbers have 78 bits, and only the last bit is different. Double precision numbers only have 53 bits. Other numbers are rounded to only 53 bits.

In your case, the two operands are rounded to the same number, and so their difference is 0.

Even stranger things happen if you slightly rewrite your expression for z:

double z = x * x * x * x - ((y * y + 1) * y * y * 4);

With this change, you get 33554432! Why? Because the way intermediate results were rounded caused the last bit of the right operand to be different. The value of the last bit is 2^(78-53)=2^25.

Answer 2

Evaluating the expression with arbitrary precision integers:

Prelude> 665857^4 - 4*(470832^4 + 470832^2)
1

Since a double normally only has 53 bits of precision and the intermediate results have 78 bits, the precision isn't sufficient to calculate the result exactly, hence it is rounded, the last bits are forgotten at some point.