Infobox math broken in current master

[Linus Torvalds]

On Tue, Jan 21, 2014 at 9:36 AM, Dirk Hohndel <dirk at hohndel.org> wrote:
>
> But given where min and max come from, we know that they are either
> identical, or at least 0.1 apart.

Actually, you don't know that.

Try this program:

    volatile double mem;
    volatile double three = 3.0;

    int main(int argc, void **argv)
    {
            double a = 1.0 / three;
            double b;

            mem = a;
            b = mem;

            return a == b;
    }

On x86-64, it will work:

    [torvalds at i7 ~]$ gcc -O2 t.c
    [torvalds at i7 ~]$ ./a.out ; echo $?
    1

but on 32-bit x86 it fails:

    [torvalds at i7 ~]$ gcc -m32 -O2 t.c
    [torvalds at i7 ~]$ ./a.out ; echo $?
    0

even though *obviously* a and b have to be the same.

Why does it fail on 32-bit x86? It fails because the i387 floating
point registers are 80-bit wide, and math is generally done in the
full precision of the registers (you can tweak that to get closer to
ieee, but it's normally not done, and you'll never get exact ieee
semantics anyway). And going through memory (which can happen as part
of a function call etc - I just used a "volatile double" to make it
explicit) will round to the normal 64-bit double.

On x86-64, you won't see that, because the math is done in the xmm
registers that don't have that extended precision.

So even "knowing where min and max come from, we know that they are
either identical, or at least 0.1 apart" is actually not true. They
might come from the same exact value, but if one of them was saved to
memory and reloaded, and the other was not, you might get an
ineuqality anyway. Floating point equality checking is basically
always a bug.

And no, I do not think it is true that "a == b" and "isnan(1 / (a-b))"
always go together. I think 1.0/0.0 returns +Inf, not NaN, and if a
and b are really small, but not equal, then a-b migth be denormal, and
then 1/(a-b) overflows.

FP math is nasty.

                Linus