Buggy floating point handling in Visual Studio?

Question

There's an update in the end!

I have a little story.

I wanted to calculate the machine epsilon (the largest epsilon > 0 satisfying condition 1.0 + epsilon = 1.0) in a C program compiled by MS Visual Studio 2008 (running on Windows 7 in a 64 bit PC). Since I know that double and float have different precision I wanted to see the answer for both. For that reason I constructed 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 quite surprised to see that it gave the same answer for both types double and float: 2.22045e-16. That was strange since double occupies twice more memory than float and should be more precise. After that I looked into Wikipedia and took a sample code from there:

    #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 figured out the following fact: my program (the first one) starts to give the 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, that is a mystery you'd say. Oh, the real mystery is the following. Compare:

  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 answer which is incorrect for float and correct for double (2.22045e-16).

In fact that is totally wrong, the result should have been opposite. That is because by default constants like 1.0 are treated as double by the compiler (according to the standard) and if it is present in an arithmetic expression then all other operands are promoted to double. Conversely, when I write 1.0f all operands are float and no promotion should take place. And yet I get a completely different result.

After all these tests I tried to compile run the programs on Linux using gcc. No surprise, it printed exactly what I expected (correct answers). So I now guess that this is Visual Studio bug. To make you all laugh (if there are any people who have read my post till that moment what 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 of 1.19209e-07.

Please somebody tell me if I am right that the root of the problem is a buggy VS 2008 compiler (can you confirm the bug on your machines?). I would be also grateful if you tested the case in a newer version: MS VS 2010. Thanks.

UPDATE. With MS Visual Studio 2013 the first program I mentioned works without unexpected results -- it gives appropriate answers for float and double. I checked this with all floating point models (precise, strict and fast) and nothing changed. So it indeed seems that VS 2008 was buggy in this case.

Answer 1

By default, Visual Studio's floating-point setting is set to "precise". This means that it will try to make the result as precise as possible. One of the side-effects of this is that intermediates are promoted to double precision.

While I haven't looked at every bit of code you posted, I suspect the problem is here:

(float) (c + eps / 2.0f)

The c + eps / 2.0f is done using double-precision. Each of the 3 operands are promoted to double-precision and the entire expression is evaluated as such. It is only rounded down to a float when you cast it.

If you set the floating-point mode to "strict", it should work as you expect it to.