gcc compiler optimization influences result of floating point comparison

Question

Problem

In using automated CI tests I found some code, which breaks if optimization of gcc is set to -O2. The code should increment a counter if a double value crosses a threshold in either direction.

Workaround

Going down to -O1 or using -ffloat-store option works around the problem.

Example

Here is a small example which shows the same problem. The update() function should return true whenever a sequence of *pNextState * 1e-6 crosses a threshold of 0.03. I used call by reference because the values are part of a large struct in the full code.

The idea behind using < and >= is that if a sequence hits the value exactly, the function should return 1 this time and return 0 the next cycle.

test.h:

extern int update(double * pState, double * pNextState);

test.c:

#include "test.h"

int update(double * pState, double * pNextState_scaled) {
    static double threshold = 0.03;
    double oldState = *pState;
    *pState = *pNextState_scaled * 1e-6;

    return oldState < threshold && *pState >= threshold;
}

main.c:

#include <stdio.h>
#include <stdlib.h>

#include "test.h"

int main(void) {

    double state = 0.01;
    double nextState1 = 20000.0;
    double nextState2 = 30000.0;
    double nextState3 = 40000.0;

    printf("%d\n", update(&state, &nextState1));
    printf("%d\n", update(&state, &nextState2));
    printf("%d\n", update(&state, &nextState3));

    return EXIT_SUCCESS;
}

Using gcc with at least -O2 the output is:

0
0
0

Using gcc with -O1, -O0 or -ffloat-store produces the desired output

0
1
0

As i understand the problem from debugging a problem arises if the compiler optimizes out the local variable oldstate on stack and instead compares against an intermediate result in an floting point register with higher precision (80bit) and the value *pState is a tiny bit smaller than the threshold. If the value for comparison is stored in 64bit precision, the logic can't miss crossing the threshold. Because of the multiplication by 1e-6 the result is probably stored in an floating point register.

Would you consider this a gcc bug? clang does not show the problem.

I am using gcc version 9.2.0 on an Intel Core i5, Windows and msys2.

Update

It is clear to me that floating point comparison is not exact and i would consider the following result as valid:

0
0
1

The idea was that, if (*pState >= threshold) == false in one cycle then comparing the same value (oldstate = *pState) against the same threshold in a subsequent call (*pState < threshold) has to be true.

Answer 1

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.

Answer 2

[Disclaimer: This is a generic, shoot-from-the-hip answer. Floating-point issues can be subtle, and I have not analyzed this one carefully. Once in a while, suspicious-looking code like this can be made to work portably and reliably after all, and per the accepted answer, that appears to be the case here. Nevertheless, the general-purpose answer stands in the, well, general case.]

I would consider this a bug in the test case, not in gcc. This sounds like a classic example of code that's unnecessarily brittle with respect to exact floating-point equality.

I would recommend either:

• rewriting the test case, or perhaps
• removing the test case.

I would not recommend:

• working around it by switching compilers
• working around it by using different optimization levels or other compiler options
• submitting compiler bug reports [although in this case, it appears there was a compiler bug, although it needn't be submitted as it has been submitted already]