Non-equivalent floating point numbers

Here’s an interesting paragraph from Handbook of Floating Point Arithmetic by Jean-Michel Muller et al.

Many common transformations of a code into some naively equivalent code become impossible if one takes into account special cases such as NaNs, signed zeros, and rounding modes. For instance, the expression x + 0 is not equivalent to the expression x if x is -0 and the rounding is to nearest.

Wait, what?!

First let’s see what this paragraph is not saying. What does the following code print?

    #include <stdio.h>

    int main()
    {
        double t = 1e-200;
        double x = -1e-200 * t; // x == -0
        double y = x + 0;

        printf( (x == y) ? "0 == 0\n" : "0 != 0\n");
    }

I compiled the code above with gcc using default options, which means rounding to nearest.

It prints 0 == 0 as you would have guessed if you weren’t spooked by the quotation above.

Let’s go back and read that paragraph more carefully. It speaks of equivalent code. In the code above, the variables x and y are equal in the sense that x == y returns true, but that does not mean that x and y are equivalent. The two variables are not equal bit-for-bit and so some code could behave differently when given one rather than the other.

Let’s define

    void hexdump(double x)
    {
        unsigned long* pu;
        pu = (unsigned long*)&x;
        printf("%lx\n", *pu);
    }

and add the lines

    hexdump(x);
    hexdump(y);

to main.

This prints the following.

    8000000000000000
    0

We see that x has its top bit set but y does not. More on why here.

Non-equivalent floating point numbers

More floating point posts

2 thoughts on “Non-equivalent floating point numbers”