Visual C++ math.h bug

Question

Hello! 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

Answer 2

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 3

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 4

How to around problem? Modify if block:

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

Answer 5

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 6

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 7

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);
}