On a float rounding error

Question

I do not understand the output of the following program:

int main()
{
 float x = 14.567729f;
 float sqr = x * x;
 float diff1 = sqr - x * x;
 double diff2 = double(sqr) - double(x) * double(x);
 std::cout << diff1 << std::endl;
 std::cout << diff2 << std::endl;
 return 0;
}

Output:

6.63225e-006
6.63225e-006

I use VS2010, x86 compiler.

I expect to get a different output

0
6.63225e-006

Why diff1 is not equal to 0? To calculate sqr - x * x compiler increases float precision to double. Why?

• ereOn: Could it be that one compiler computes the result (or a part of it) at compile time using different "rules" than at runtime ?
• Alexey Malistov: @ereOn. It is not true. Assembly code does not confirm your suggestion.

Answer 1

The floating point registers are 80 bits (on most modern CPUs)

During an expression the result is an 80 bit value. It only gets truncated to 32 (float) or 64 (double) when it gets assigned to a location in memory. If you hold everything in registers (try compiling with -O3) you may see a different result.

Compiled with: -03:

> ./a.out
0
6.63225e-06

• David Hammen: Re "on most modern CPUs": Better stated as "on the family of CPUs that are most widely-used today". There are some machines that still treat floats and doubles differently all the way down to the machine level.
• Loduwijk: Isn't that assuming that the value is not being calculated in the register anymore by the time the code that represents "sqr-x*x" is being executed? I almost included this in my answer as well, but then decided not to since I thought the automatic promotion to double during arithmetic scenario seemed like a more likely culprit. I would have assumed that the compiler would continue to operate on the value in the register at that point. Although, now that I think about it, the question said "VS2010 compiler" and I have seen windows compilers be stupid about register usage before.
• edA-qa mort-ora-y: While your idea of truncation is correct you're wrong about register size. Common chips offer 128-bit registers (check SSIMD, and long double on 64-bit), as well as some having actual 64-bit registers. Compiler optimizers will mix and match the registers as they see fit.
• Loki Astari: @edA: Yes the standard registers are getting longer. Its been a while since I bothered to check but the last time the floating point instructions were still based of the x87 (not x86). Where the FPU floating point register stack was 80 bits wide (though with 64 bit processors now common I am sure you can find larger (but I have not checked)). The point of the answer being that that floating point register size is different from the floating point memory size and completely seprate from the standard registeres that you are talking about. en.wikipedia.org/wiki/X87
• Jan Hudec: @DavidHammen: I just happened to look in the relevant header and it seems x86-64 and x86 with SSE math both use precision of the operand. That makes "on most modern CPUs" false and "family of CPUs that are most widely used" false in near future (for Linux most code is compiled 64-bit where available these days, but for Windows people are still often compiling 32-bit).

Answer 2

float diff1 = sqr - x * x;
double diff2 = double(sqr) - double(x) * double(x);

Why diff1 is not equal to 0?

Because you have already cached sqr = x*x and forced its representation to be a float.

To calculate sqr - x * x compiler increases float precision to double. Why?

Because that is how C did things back before there was a C standard. I don't think modern compilers are bound to that convention, but many still do follow it. If this is the case, the right-hand sides of the calculations of diff1 and diff2 will be identical. The only difference is that after calculating the right-hand side of float diff1 = ..., the double result is converted back to a float.

• David Hammen: Found it in the C99 standard, but not in the C++: "The values of floating operands and of the results of floating expressions may be represented in greater precision and range than that required by the type; the types are not changed thereby." In other words, vendors are free to temporarily promote the floats to doubles in float diff1 = sqr - x*x;so long as the the types/values of the right-hand side variables aren't change and sol long as the end result is stored in diff1 as a float.

Answer 3

Apparently the standard allows floats to be automatically promoted to double in expressions like that. See here

Do a find on that page for "automatically promoted" and check out the first paragraph with that phrase in it.

If we go by that paragraph, as I understand it, your sqr=x*x is initially being treated as if it were a double as well, but once it is stored it is being rounded to a float. Then, in your diff1=sqr-x*x, x*x is again being treated like a double, and so is sqr although it's already rounded. Therefore, it yields the same result as casting them all to doubles: sqr is a double then but already rounded to float precision, and again x*x is double precision.

• KillianDS: also, do not forget optimizations, it could be that this is reduced to one big constant expression, where floats indeed can be promoted to doubles

Answer 4

On x86/x64 architectures it is common for compilers to promote all 32-bit floats to 64-bit doubles for computations; check the output assembly to see if the two variants produce the same instructions. The only difference between the types is the storage.