Different floating point result with optimization enabled - compiler bug?
Different compilers have different optimization settings. Some of those faster optimization settings do not maintain strict floating-point rules according to IEEE 754. Visual Studio has a specific setting, /fp:strict, /fp:precise, /fp:fast, where /fp:fast violates the standard on what can be done. You might find that this flag is what controls the optimization in such settings. You may also find a similar setting in GCC which changes the behaviour.
If this is the case then the only thing that's different between the compilers is that GCC would look for the fastest floating point behaviour by default on higher optimisations, whereas Visual Studio does not change the floating point behaviour with higher optimization levels. Thus it might not necessarily be an actual bug, but intended behaviour of an option you didn't know you were turning on.
Output should be: 4.5 4.6 That's what the output would be if you had infinite precision, or if you were working with a device that used a decimal-based rather than binary-based floating point representation. But, you aren't. Most computers use the binary IEEE floating point standard.
As Maxim Yegorushkin already noted in his answer, part of the problem is that internally your computer is using an 80 bit floating point representation. This is just part of the problem, though. The basis of the problem is that any number of the form n.nn5 does not have an exact binary floating representation. Those corner cases are always inexact numbers.
If you really want your rounding to be able to reliably round these corner cases, you need a rounding algorithm that addresses the fact that n.n5, n.nn5, or n.nnn5, etc. (but not n.5) is always inexact. Find the corner case that determines whether some input value rounds up or down and return the rounded-up or rounded-down value based on a comparison to this corner case. And you do need to take care that a optimizing compiler will not put that found corner case in an extended precision register.
See How does Excel successfully Rounds Floating numbers even though they are imprecise? for such an algorithm.
Or you can just live with the fact that the corner cases will sometimes round erroneously.
To those who can't reproduce the bug: do not uncomment the commented out debug stmts, they affect the result.
This implies that the problem is related to the debug statements. And it looks like there's a rounding error caused by loading the values into registers during the output statements, which is why others found that you can fix this with -ffloat-store
Further question:
I was wondering, should I always turn on
-ffloat-storeoption?
To be flippant, there must be a reason that some programmers don't turn on -ffloat-store, otherwise the option wouldn't exist (likewise, there must be a reason that some programmers do turn on -ffloat-store). I wouldn't recommend always turning it on or always turning it off. Turning it on prevents some optimizations, but turning it off allows for the kind of behavior you're getting.
But, generally, there is some mismatch between binary floating point numbers (like the computer uses) and decimal floating point numbers (that people are familiar with), and that mismatch can cause similar behavior to what your getting (to be clear, the behavior you're getting is not caused by this mismatch, but similar behavior can be). The thing is, since you already have some vagueness when dealing with floating point, I can't say that -ffloat-store makes it any better or any worse.
Instead, you may want to look into other solutions to the problem you're trying to solve (unfortunately, Koenig doesn't point to the actual paper, and I can't really find an obvious "canonical" place for it, so I'll have to send you to Google).
If you're not rounding for output purposes, I would probably look at std::modf() (in cmath) and std::numeric_limits<double>::epsilon() (in limits). Thinking over the original round() function, I believe it would be cleaner to replace the call to std::floor(d + .5) with a call to this function:
// this still has the same problems as the original rounding function
int round_up(double d)
{
// return value will be coerced to int, and truncated as expected
// you can then assign the int to a double, if desired
return d + 0.5;
}
I think that suggests the following improvement:
// this won't work for negative d ...
// this may still round some numbers up when they should be rounded down
int round_up(double d)
{
double floor;
d = std::modf(d, &floor);
return floor + (d + .5 + std::numeric_limits<double>::epsilon());
}
A simple note: std::numeric_limits<T>::epsilon() is defined as "the smallest number added to 1 that creates a number not equal to 1." You usually need to use a relative epsilon (i.e., scale epsilon somehow to account for the fact that you're working with numbers other than "1"). The sum of d, .5 and std::numeric_limits<double>::epsilon() should be near 1, so grouping that addition means that std::numeric_limits<double>::epsilon() will be about the right size for what we're doing. If anything, std::numeric_limits<double>::epsilon() will be too large (when the sum of all three is less than one) and may cause us to round some numbers up when we shouldn't.
Nowadays, you should consider std::nearbyint().
Intel x86 processors use 80-bit extended precision internally, whereas double is normally 64-bit wide. Different optimization levels affect how often floating point values from CPU get saved into memory and thus rounded from 80-bit precision to 64-bit precision.
Use the -ffloat-store gcc option to get the same floating point results with different optimization levels.
Alternatively, use the long double type, which is normally 80-bit wide on gcc to avoid rounding from 80-bit to 64-bit precision.
man gcc says it all:
-ffloat-store
Do not store floating point variables in registers, and inhibit
other options that might change whether a floating point value is
taken from a register or memory.
This option prevents undesirable excess precision on machines such
as the 68000 where the floating registers (of the 68881) keep more
precision than a "double" is supposed to have. Similarly for the
x86 architecture. For most programs, the excess precision does
only good, but a few programs rely on the precise definition of
IEEE floating point. Use -ffloat-store for such programs, after
modifying them to store all pertinent intermediate computations
into variables.
In x86_64 builds compilers use SSE registers for float and double by default, so that no extended precision is used and this issue doesn't occur.
gcc compiler option -mfpmath controls that.