Cosf (M_PI_2) without returning zero

Question

It started abruptly this morning.

The original lines were

float angle = (x+90)*(M_PI/180.0); float xx = cosf(angle); float yy = sinf(angle); 

After setting a breakpoint and turning the cursor. I get the correct answer for yy as 1. but xx is not zero.

I tried with cosf(M_PI_2); still no luck .. it worked fine until yesterday. I did not change the compiler settings, etc.

I am using the latest version of Xcode today

Answer 1

Unlike others, this is not an x87 coprocessor problem. Xcode uses SSE by default to calculate the default floating point (with the exception of long double arithmetic).

"Problem": when you write cosf(M_PI_2) , you are actually telling the cosf(M_PI_2) compiler (gcc or llvm-gcc or clang) to do the following:

• Locate the extension M_PI_2 in <math.h> . According to the POSIX standard, it is a double-precision literal that translates to a correctly rounded π / 2 value.
• Complete converted double precision value to single precision.
• Call cosf math library cosf with a single precision value.

Please note that during this process you are not working with the actual value of π / 2. Instead, you work with this value, rounded to a number representing a floating point number. While cos (π / 2) is zero, you are not telling the compiler to do this calculation. Instead, you tell the compiler to make cos (π / 2 + tiny), where is the tiny difference between the rounded value (float)M_PI_2 and the (non-representable) exact value of π / 2. If cos calculated without any error, the result of cos (π / 2 + tiny) is approximately -tiny. If it returns zero , that will be an error.

edit: step-by-step expansion of calculations on Intel layout using the current Xcode compiler:

M_PI_2 defined as

 1.57079632679489661923132169163975144 

but it’s not really a double precision number. When the compiler converts it to a double precision value, it becomes exactly

 1.5707963267948965579989817342720925807952880859375 

This is the closest double precision number to π / 2, but it differs from the actual mathematical value of π / 2 by about 6.12 * 10 ^ (- 17).

Step (2) rounds this number to a single point, which changes the value to exactly

 1.57079637050628662109375 

This is approximately π / 2 + 4.37 * 10 ^ (- 8). When we calculate the cosf this number, we get:

 -0.00000004371138828673792886547744274139404296875 

which is the almost exact value of the cosine estimated at this point:

 -0.00000004371139000186241438857289400265215231661... 

In fact, this is a properly rounded result; there is no value that calculation could return that would be more accurate. The only mistake here is that the calculation that you asked the compiler to do is different from the calculations that you thought you asked to do.

Answer 2

The first thing to notice is the use of float s. They are inherently inaccurate, and for most calculations you get only a close approximation of the mathematically correct answer. Assuming x in your code is 0 , angle will have a close approximation to π / 2. xx will therefore have an approximation to cos (π / 2). However, this is unlikely to be zero due to approximation and rounding problems.

If you were able to change your code to us double , not float , you are likely to get more accuracy and the answer is closer to zero. However, if it is important for this code to have a value of zero at the moment, you will have to rethink how you perform the calculations.

If this does not answer your specific problem, give us more detailed information and we will have a different opinion.

Answer 3

I suspect the answer is as close as damnit to 0, which is not worth worrying about.

If I run the same, I get the answer "-4.3711388e-008", which can also be written as "-0.000000043711388". Which is pretty damned close to 0. Definitely close enough not to worry about being in the eighth decimal place.

Edit: In addition to what LiraLuna says, I wrote the following x87 assembler snippet under visual studio

  float fRes; _asm { fld1 fld1 fadd st, st(1) fldpi fdiv st, st(1) fcos fstp [fRes] } char str[16]; sprintf( str, "%f", fRes ); 

This basically uses the x87 fcos instruction to create the cosine of pi / 2. the value held on str is "0.000000"

This, however, is not really what fcos returns. He STRONGLY returned 6.1230318e-017. This means that the error occurs in the 17th decimal place and, frankly, is much less significant than the standard debug cosf above.

Since SSE3 does not have a specific cosine command, I suspect (although I cannot confirm without seeing the assembler assembly) that it either uses its own taylor series extension or uses the fcos instruction. In any case, in my opinion, you are unlikely to get better accuracy than the error on the 17th decimal place.

Answer 4

The only thing I can think of is malicious macro substitution, i.e. M_PI_2 no longer 1.57079632679489661923.

Try calling cosf( 1.57079632679489661923 ) to check this.

Answer 5

The real thing you have to be careful about is the sign cosine. Make sure it is the same as you expected. For instance. if you work with angles between 0 and pi / 2. make sure that what you use as PI_2 is less than the actual pi / 2 value !

And the difference between 0.000001 and 0.0 less than you think.

Answer 6

Cause

What you're experiencing is the infamous x87 math co-processor float truncate 'bug' - or rather, a function. IEEE float have an amazing range of numbers, but at a price. They sacrifice high-level precession.

They are not inaccurate, as you think, although this is a half-myth created by the Intel x87 chip, which internally uses the internal representation of 80 bits for floats - they have a much higher precession, although a little slower.

When you perform a comparison with a float, x87 caches the float as an 80-bit float, then when it is full, it stores the 32-bit representation in RAM, reducing accuracy to a large extent.

Decision

x87 old, really old. This is a replacement for SSE . SSE computes 32-bit floats and 64-bit floats initially, which leads to minimal precession lost in math. Note that precession problems with floats still exist, but printf("%f\n", cosf(M_PI_2)); must be zero. Heck - even a smooth comparison with SSE is accurate again! (unlike x87).

Since the last Xcode is actually GCC 4.2.1, use the -msse3 -mfpmath=sse compiler switch and see how you get a perfectly round 0.00000 (Note: if you get -0.00000 , don’t worry, it’s fine and it’s equal to 0.00000 in the IEEE specification (more on this wikipedia ).

All Mac Macs are guaranteed SSE3 support (Mac OSx86 are excluded; if you want to support them, use -msse2 ).