Visual C ++ math.h error

Question

I debugged my project and could not find the error. I finally found it. Look at the code. You think that everything is in order, and the result will be “OK, good! Good!”, Right? Now compile it with VC (I 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 fine, when x> 90112.0 - No! You can change cos to sin.

Any ideas? Remember that sin and cos are periodic.

Answer 1

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

I know this is hard to accept, but floating point arithmetic just doesn't work as many people expect. Worse, some differences depend on the details of your particular floating-point computer and / or the optimization settings that you use in your particular compiler. You may not like it, but the way it is. The only way to “get it” is to put aside your assumptions about how things should behave and accept things, how they 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 equipment, which often contains special registers, and these registers often have more bits than a double . This means that calculations of intermediate floating point calculations often have more bits than sizeof(double) , and when a floating point value is written to RAM, it is often truncated, often losing some precision bits ...

just keep that in mind: floating point comparisons are complex and subtle and fraught with danger. Be careful. The way the floating point works is different from the way most programmers tend to think this should work. If you intend to use floating points, you need to find out how it actually works ...

Answer 2

As others have noted, the VS math library performs its calculations on the x87 FPU and generates 80-bit results, even if the type is double.

In this way:

• cos () is called and returns with cos (x) at the top of the x87 stack as an 80-bit float
• cos (x) is popped from the x87 stack and stored in memory as double; this results in rounding to 64 bit float, which changes its value
• cos () is called and returns with cos (x) at the top of the x87 stack as an 80-bit float
• rounded value is loaded onto the x87 stack from memory
• rounded and non-circular values ​​of cos (x) are compared unevenly.

Many math libraries and compilers protect you from this by performing a 64-bit float calculation in SSE registers when available, or by force saving and rounding values ​​before comparison, or by saving and reloading the final result into the actual cos () calculation. The compiler / library combination you work with is not very forgiving.

Answer 3

You should never compare deuces for equality in most cases. You cannot get what you expect.

Floating point registers can have different sizes than memory values ​​(in modern intel machines, FPU registers are 80 bits and 64 bits). If the compiler generates code that calculates the first cosine, then stores the value in memory, computes the second cosine and compares the value in memory with the register, then the values ​​may differ (due to rounding errors from 80 to 64 bits),

Floating point values ​​are a bit complicated. Google for floating point comparison.

Answer 4

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

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

The value is calculated once, and then cloned, then compared with itself - the result will be OK

Same thing 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          

Strange things are happening now.
1. X is uploaded to fstack (X, 0)
2. X is stored on a regular stack (truncation)
3. Cosine is calculated, the result on the fleet stack
4. boot X again

5. X is stored on a regular stack (truncation, since now we are "symmetrical")
6. The result of the 1st cosine that was on the stack is stored in memory, now another truncation occurs for the 1st value

7. Cosine is calculated, the 2nd result, if on a float-stack, but this value was truncated only once
8. The first value is uploaded to fstack, but this value has been truncated twice (once before calculating the cosine, once after)
9. These 2 values ​​are compared - we get rounding errors.

Answer 5

The compiler may have generated code that finishes comparing a 64-bit double value with an 80-bit internal floating-point register. Testing floating point values ​​for equality is prone to similar errors - you are almost always better at performing a “fuzzy” comparison, for example (fabs (val1 - val2) <EPSILON), and not (val1 == val2).

Answer 6

Incrementing and testing a float value as a loop control variable is usually a really bad idea. Create a separate int LCV only for the loop if you need to.

In this case, it is simpler:

 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 get around the problem? Change if :

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