In the end, an update!
I have a little story.
I wanted to calculate the epsilon machine (the largest epsilon> 0, satisfying the condition 1.0 + epsilon = 1.0) in a C program compiled by MS Visual Studio 2008 (runs on Windows 7 on a 64-bit PC). Since I know that double and float have different precision, I wanted to see the answer for both. For this reason, I built the following program:
#include <stdio.h> typedef double float_type; int main() { float_type eps = 1.0; while ((float_type) 1.0 + eps / (float_type) 2.0 > (float_type) 1.0) eps = eps / (float_type) 2.0; printf("%g\n", eps); return 0; }
I was very surprised to see that it gave the same answer for both double and float types: 2.22045e-16. This was strange, since a double takes up twice as much memory as a floating one, and should be more accurate. After that, I looked at Wikipedia and took a sample code:
#include <stdio.h> int main(int argc, char **argv) { float machEps = 1.0f; do { machEps /= 2.0f; } while ((float)(1.0 + (machEps/2.0)) != 1.0); printf( "\nCalculated Machine epsilon: %G\n", machEps ); return 0; }
I was even more surprised when it worked correctly! After some attempts to understand the fundamental difference between the two programs, I realized the following fact: my program (first) starts giving a corrent answer for float (1.19209e-07) if I change the loop condition to
while ((float_type) (1.0 + eps / (float_type) 2.0) > (float_type) 1.0)
Well, this is the secret you will say. Oh, the real secret is this. For comparison:
while ((float) (1.0 + eps / 2.0f) > 1.0f)
which gave the correct answer (1.19209e-07) and
while ((float) (1.0f + eps / 2.0f) > 1.0f)
which gave an answer that is incorrect for float and correct for double (2.22045e-16).
In fact, this is absolutely wrong, the result should be the opposite. This is due to the fact that by default constants like 1.0 are processed by the compiler as double (according to the standard), and if it is present in an arithmetic expression, then all other operands are increased to two. Conversely, when I write 1.0f, all operands are floating, and promotions should not occur. And yet I get a completely different result.
After all these tests, I tried to compile programs on Linux using gcc. No wonder he printed exactly what I expected (correct answers). So now I assume this is a Visual Studio error. So that you all laugh (if there are people who have read my post up to this point, which is doubtful ^ _ ^), I will give you another comparison:
float c = 1.0; while ((float) (c + eps / 2.0f) > 1.0f)
This does not work properly in VS, but ...
const float c = 1.0;
gives the correct answer 1.19209e-07.
Please tell me if I am right that the root of the problem is the compiler with the VS 2008 error (can you acknowledge the error on your machines?). I would also appreciate if you tested the case in a newer version: MS VS 2010. Thank you.
UPDATE With MS Visual Studio 2013, the first program I mentioned works without unexpected results - it gives the corresponding answers for float and double. I checked this with all floating point models (accurate, rigorous and fast) and nothing has changed. Therefore, it really seems that VS 2008 was a mistake in this case.