Different Floating Point Result With Optimization Enabled - Compiler Bug

Different floating point result with optimization enabled - compiler bug?

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.

gcc compiler optimization influences result of floating point comparison

I've analysed your code and my take is that it is sound according the standard but you are hit by gcc bug 323 about which you may find more accessible information in gcc FAQ.

A way to modify your function and render it robust in presence of gcc bug is to store the fact that the previous state was below the threshold instead (or in addition) of storing that state. Something like this:

int update(int* pWasBelow, double* pNextState_scaled) {
    static double const threshold = 0.03;
    double const nextState = *pNextState_scaled * 1e-6;
    int const wasBelow = *pWasBelow;
    *pWasBelow = nextState < threshold;
    return wasBelow && !*pWasBelow;
}

Note that this does not guarantee reproductibility. You may get 0 1 0 in one set-up and 0 0 1 in another, but you'll detect the transition sooner or later.

Why can the compiler not optimize floating point addition with 0?

IEEE 754 floating-point numbers have two zero values, one negative, one positive. When added together, the result is the positive one.

So id1(-0.f) is 0.f, not -0.f.

Note that id1(-0.f) == -0.f because 0.f == -0.f.

Demo

Also, note that compiling with -ffast-math in GCC does make the optimization and changes the result.

floating point rounding giving different results for 80-bit register and 64-bit double: ill-formed code or gcc/clang bug?

Without -ffloat-store, GCC targeting x87 does violate the standard: it keeps values un-rounded even across statements. (-mfpmath=387 is the default for -m32). Assignment like double x = y; is supposed to round to actual double in ISO C++, and probably also passing a function arg.

So I think your code is safe for ISO C++ rules, even with the FLT_EVAL_METHOD == 2 that GCC claims to be doing. (https://en.cppreference.com/w/cpp/types/climits/FLT_EVAL_METHOD)

See also https://randomascii.wordpress.com/2012/03/21/intermediate-floating-point-precision/ for more about the real-world issues, with actual compilers for x86.

https://gcc.gnu.org/wiki/x87note doesn't really mention the difference between when GCC rounds vs. when ISO C++ requires rounding, just describes GCC's actual behaviour.

How does -march native affect floating point accuracy?

The use of FMA can both decrease and increase error, both of those may result in a testcase failing, depending on how the test works. FMA improves error "locally", but the effect may be the opposite when put in a wider context.

For example, a * c - b * d (determinant of a 2x2 matrix) famously gives some (usually minor) trouble when FMA-contracted. Without FMA, the subtraction has the potential to eliminate the rounding error, if it is the same on both sides. That does not always happen, but it can happen when a * c = b * d which is of special interest because that means the determinant should be zero. Without FMA the result would actually be zero, with FMA it won't be.

#include <math.h>
#include <stdio.h>

double determinant(double a, double b, double c, double d)
{
    return a * c - b * d;
}

int main()
{
    volatile double a = M_PI;
    double x = determinant(a, a, a, a);
    printf("%E\n", x);
    return 0;
}

This program, compiled by GCC 11.2 with optimizations enabled and FMA allowed, does not print zero, but something on the order of 1E-16.

Some variants of an "is this result close enough"-test in a unit-test would conclude that this result is, relative to zero, extremely wrong. Another way to look at it though, is that if one of the inputs changed by just 1 ULP, that would have introduced an error on the order of 1E-15, which is even worse.

Most special/new instructions either don't affect accuracy or are by default restricted. For example addsubpd and haddpd (from SSE3) are just equivalents of what would have cost more code before, and roundpd (from SSE4.1) is by default only used in ways that don't affect results (using roundpd for floor and ceil is safe, ironically using it for round itself is non-trivial due to the different halfway rounding).

-O1 alters floating point math

With -O1, the floating computation happens at compile time, using the GNU MPFR library. MPFR is expected to give a correctly rounded result even for functions such as sin and cos. Your math library likely has different accuracy goals for these functions, which is why run-time computation (at the -O0 optimization level) sometimes gives different results. For example, the GNU C library has a general accuracy goal of a few ulp.

Reportedly, IEEE 754 only has accuracy requirements for a subset of the math library functions (sqrt, apparently), which enables math libraries to choose different trade-offs between speed and accuracy for the transcendental functions. (I do not have access to IEEE 754 because IEEE is opposed to the open dissemination of knowledge unfortunately.)

---

Related Topics