Different floating point result with optimization enabled - compiler bug?

The below code works on Visual Studio 2008 with and without optimization. But it only works on g++ without optimization (O0).

#include <cstdlib>
#include <iostream>
#include <cmath>

double round(double v, double digit)
{
    double pow = std::pow(10.0, digit);
    double t = v * pow;
    //std::cout << "t:" << t << std::endl;
    double r = std::floor(t + 0.5);
    //std::cout << "r:" << r << std::endl;
    return r / pow;
}

int main(int argc, char *argv[])
{
    std::cout << round(4.45, 1) << std::endl;
    std::cout << round(4.55, 1) << std::endl;
}

The output should be:

4.5
4.6

But g++ with optimization (O1 - O3) will output:

4.5
4.5

If I add the volatile keyword before t, it works, so might there be some kind of optimization bug?

Test on g++ 4.1.2, and 4.4.4.

Here is the result on ideone: http://ideone.com/Rz937

And the option I test on g++ is simple:

g++ -O2 round.cpp

The more interesting result, even I turn on /fp:fast option on Visual Studio 2008, the result still is correct.

Further question:

I was wondering, should I always turn on the -ffloat-store option?

Because the g++ version I tested is shipped with CentOS/Red Hat Linux 5 and CentOS/Redhat 6.

I compiled many of my programs under these platforms, and I am worried it will cause unexpected bugs inside my programs. It seems a little difficult to investigate all my C++ code and used libraries whether they have such problems. Any suggestion?

Is anyone interested in why even /fp:fast turned on, Visual Studio 2008 still works? It seems like Visual Studio 2008 is more reliable at this problem than g++?

Solution 1:

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.

Solution 2:

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.