cosf(M_PI_2) not returning zero

Question

This started suddenly today morning.

Original lines were this

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

After putting a breakpoint and hovering 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 was working fine till yesterday.. I did not change any compiler setting etc..

I am using Xcode latest version as of todays date

Answer 1

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

Try calling cosf( 1.57079632679489661923 ) to test this.

Answer 2

The first thing to notice is that you're using floats. These are inherently inaccurate, and for most calculations give you only a close approximation of the mathematically-correct answer. Assuming that x in your code has value 0, angle will have a close approximation to π/2. xx will therefore have an approximation to cos(π/2). However, this is unlikely to be exactly zero due to approximation and rounding issues.

If you were able to change your code to us doubles rather than floats you're likely to get more accuracy, and an answer nearer zero. However, if it is important for your code to produce a value of exactly zero at this point, you're going to have to rethink how you're doing the calculations.

If this doesn't answer your particular problem, give us some more details and we'll have another think.

Answer 3

I suspect the answer is as near as damnit to 0 as not to be worth worrying about.

If i run the same thing through I get the answer "-4.3711388e-008" which can also be written as "-0.000000043711388". Which is pretty damned close to 0. Definitely near enough to not worry about it being out at the 8th decimal place.

Edit: Further to what LiraLuna is saying I wrote the following piece of x87 assembler 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 );

Basically this uses the x87's fcos instruction to do a cosine of pi/2. the value held in str is "0.000000"

This, however, is not actually what fcos returned. It ACTUALLY returned 6.1230318e-017. This implies that the error occurs at the 17th decimal place and, lets be honest, thats far less significant than the standard debug cosf above.

As SSE3 has no specific cosine instruction I suspect (though i cannot confirm without seeing the assembler generated) that it is either using its own taylor series expansion or it is using the fcos instruction anyway. Either way you are still unlikely to get better precision than the error occurring at the 17th decimal place, in my opinion.

Answer 4

The real thing you should be careful about is the sign of cosine. Make sure it is the same as you expected. E.g. if you operate with angles between 0 and pi/2. make sure that what you use as PI_2 is less that actual value of pi/2!

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

Answer 5

The reason

What you are experiencing is the infamous x87 math co-processor float truncate 'bug' - or rather - a feature. IEEE floats have an amazing range of numbers, but at a cost. They sacrifice precession for high range.

They are not inaccurate as you think, though - this is a semi-myth generate by Intel's x87 chip design, that internally uses 80bit internal representation for floats - they have far superior precession though a bit slower.

When you perform a float comparison, x87 caches the float as an 80bit float, then when it's stack is full, it saves the 32bit representation in RAM, decreasing accuracy by a large degree.

The solution

x87 is old, really old. It's replacement is SSE. SSE computes 32bit floats and 64bit floats natively, leading to minimal precession lost on math. Please note that precession issues with floats still exist, but printf("%f\n", cosf(M_PI_2)); should be zero. Heck - even float comparison with SSE is accurate again! (unlike x87).

Since latest Xcode is actually GCC 4.2.1, use the compiler switch -msse3 -mfpmath=sse and see how you get a perfectly round 0.00000 (Note: if you get -0.00000, do not worry, it's perfectly fine and still equals 0.00000 under the IEEE spec (read more at this wikipedia article)).

All Intel macs are guaranteed to have SSE3 support (OSx86 Macs excluded, if you want to support those, use -msse2).

Answer 6

Contrary to what others have said, this is not an x87 co-processor issue. XCode uses SSE for floating-point computation on Intel by default (except for long double arithmetic).

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

1. Look up the expansion of M_PI_2 in <math.h>. Per the POSIX standard, it is a double precision literal that converts to the correctly rounded value of π/2.
2. Round the converted double precision value to single precision.
3. Call the math library function cosf on the single precision value.

Note that, throughout this process, you are not operating on the actual value of π/2. You are instead operating on that value rounded to a representable floating-point number. While cos(π/2) is exactly zero, you are not telling the compiler to do that computation. You are instead telling the compiler to do cos(π/2 + tiny), where tiny is the difference between the rounded value (float)M_PI_2 and the (unrepresentable) exact value of π/2. If cos is computed with no error at all, the result of cos(π/2 + tiny) is approximately -tiny. If it returned zero, that would be an error.

edit: a step-by-step expansion of the computation on an Intel mac with the current XCode compiler:

M_PI_2 is defined to be

1.57079632679489661923132169163975144

but that's not actually a representable 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 single-precision, which changes the value to exactly

1.57079637050628662109375

Which is approximately π/2 + 4.37*10^(-8). When we compute cosf of this number then, we get:

-0.00000004371138828673792886547744274139404296875

which is very nearly the exact value of cosine evaluated at that point:

-0.00000004371139000186241438857289400265215231661...

In fact, it is the correctly rounded result; there is no value that the computation could have returned that would be more accurate. The only error here is that the computation that you asked the compiler to perform is different from the computation that you thought you were asking it to do.