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.

Answer 3

The cos(x) == cos(x) procedure generated in release mode:

00DB101A  call        _CIcos (0DB1870h) 
00DB101F  fld         st(0) 
00DB1021  fucompp 

The value is computed once, and then cloned, then compared with itself - result will be ok

The same in debug mode:

00A51405  sub         esp,8 
00A51408  fld         qword ptr [x] 
00A5140B  fstp        qword ptr [esp] 
00A5140E  call        @ILT+270(_cos) (0A51113h) 
00A51413  fld         qword ptr [x] 
00A51416  fstp        qword ptr [esp] 
00A51419  fstp        qword ptr [ebp-0D8h] 
00A5141F  call        @ILT+270(_cos) (0A51113h) 
00A51424  add         esp,8 
00A51427  fld         qword ptr [ebp-0D8h] 
00A5142D  fucompp          

Now, strange things happen.
1. X is loaded on fstack (X, 0)
2. X is stored on normal stack (truncation)
3. Cosine is computed, result on float stack
4. X is loaded again
5. X is stored on normal stack (truncation, as for now, we are "symmetrical")
6. The result of 1st cosine that was on the stack is stored in memory, now, another truncation occurs for the 1st value
7. Cosine is computed, 2nd result if on the float-stack, but this value was truncated only once
8. 1st value is loaded onto the fstack, but this value was truncated twice (once before computing cosine, once after)
9. Those 2 values are compared - we're getting rounding errors.

Answer 4

You should never not compare doubles for equality in most cases. You may not get what you expect.

Floating point registers can have a different size than memory values (in current intel machines, FPU registers are 80 bit vs 64 bit doubles). If the compiler is generating code that calculates the first cosine, then stores the value into memory, calculates the second cosine and compares the value in memory from that in the register then the values could differ (due to rounding issues from 80 to 64 bits).

Floating point values are a bit tricky. Google for floating point comparissons.

Answer 5

The compiler might have generated code that ends up comparing a 64-bit double value with an 80-bit internal floating point register. Testing floating point values for equality is prone to these sorts of errors -- you're almost always better off doing a "fuzzy" comparison like (fabs(val1 - val2) < EPSILON) rather than (val1 == val2).

Answer 6

Incrementing and testing a float value as a loop control variable is generally a really Bad Idea. Create a separate int LCV just for looping on, if you have to.

In this case it's easier:

for ( int i = 90100; i<90120; i+=1 )    {
    if ( cos(i) == cos(i) )
        printf ("i==%d  OK!\n", i);
    else
        printf ("i==%d  FAIL!\n", i);
}

Answer 7

How to around problem? Modify if block:

if ( (float)cos(x) == (float)cos(x) )