Visual C++ math.h bug

Question

I was debugging my project and could not find a bug. Finally I located it. Look at the code. You think everything is OK, and result will be "OK! OK! OK!", don't you? Now compile it with VC (i've tried vs2005 and vs2008).

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


int main () {
    for ( double x = 90100.0; x<90120.0; x+=1 )
    {
        if ( cos(x) == cos(x) )
            printf ("x==%f  OK!\n", x);
        else
            printf ("x==%f  FAIL!\n", x);
    }

    getchar();
    return 0; 
}

The magic double constant is 90112.0. When x < 90112.0 everything is OK, when x > 90112.0 -- Nope! You can change cos to sin.

Any ideas? Dont forget that sin and cos are periodic.

Answer 1

Could be this: http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.18

I know it's hard to accept, but floating point arithmetic simply does not work like most people expect. Worse, some of the differences are dependent on the details of your particular computer's floating point hardware and/or the optimization settings you use on your particular compiler. You might not like that, but it's the way it is. The only way to "get it" is to set aside your assumptions about how things ought to behave and accept things as they actually do behave...

(with emphasis on the word "often"; the behavior depends on your hardware, compiler, etc.): floating point calculations and comparisons are often performed by special hardware that often contain special registers, and those registers often have more bits than a double. That means that intermediate floating point computations often have more bits than sizeof(double), and when a floating point value is written to RAM, it often gets truncated, often losing some bits of precision...

just remember this: floating point comparisons are tricky and subtle and fraught with danger. Be careful. The way floating point actually works is different from the way most programmers tend to think it ought to work. If you intend to use floating point, you need to learn how it actually works...

Answer 2

As others have noted, the VS math library is doing its computation on the x87 FPU, and generating 80-bit results even though the type is double.

Thus:

1. cos( ) is called, and returns with cos(x) on the top of the x87 stack as an 80bit float
2. cos(x) is popped off the x87 stack and stored to memory as a double; this causes it to be rounded to 64bit float, which changes its value
3. cos( ) is called, and returns with cos(x) on the top of the x87 stack as an 80bit float
4. the rounded value is loaded onto the x87 stack from memory
5. the rounded and unrounded values of cos(x) compare unequal.

Many math libraries and compilers protect you from this by either doing the computation in 64bit float in the SSE registers when available, or by forcing the values to be stored and rounded before the comparison, or by storing and reloading the final result in the actual computation of cos( ). The compiler/library combination you happen to be working with isn't so forgiving.