Infinite loop heisenbug: it exits if I add a printout

Question

This is my source code:

#include <iostream>
#include <cmath>

using namespace std;

double up = 19.0 + (61.0/125.0);
double down = -32.0 - (2.0/3.0);
double rectangle = (up - down) * 8.0;

double f(double x) {
    return (pow(x, 4.0)/500.0) - (pow(x, 2.0)/200.0) - 0.012;
}

double g(double x) {
    return -(pow(x, 3.0)/30.0) + (x/20.0) + (1.0/6.0);
}

double area_upper(double x, double step) {
    return (((up - f(x)) + (up - f(x + step))) * step) / 2.0;
}

double area_lower(double x, double step) {
    return (((g(x) - down) + (g(x + step) - down)) * step) / 2.0;
}

double area(double x, double step) {
    return area_upper(x, step) + area_lower(x, step);
}

int main() {
    double current = 0, last = 0, step = 1.0;
    
    do {
        last = current;
        step /= 10.0;
        current = 0;
        
        for(double x = 2.0; x < 10.0; x += step) current += area(x, step);
        
        current = rectangle - current;
        current = round(current * 1000.0) / 1000.0;
        //cout << current << endl;
    } while(current != last);
    
    cout << current << endl;
    return 0;
}

What it does is calculating area between the curves. There is a loop in main() - its purpose is to calculate value as precise as it's possible within 3 decimal places.

It didn't work. For the sake of debugging, I added the line which is the only commented one. I wanted to know what's going on inside the loop.

//cout << current << endl;

When the line is there - when I uncomment it - everything works fine. When it's not - the loop seems to be infinite.

Holy Compiler, why?

It's not the matter of imprecise floating-point numbers, which I am aware of. Everything is finished within 4 repeats of the loop content when I'm outputting current value inside it.

Answer 1

@Skizz's comment gives the likely problem, but to elaborate:

Floating point math is tricky, and, in particular, rounding errors can often arise. A number such as 1/1000.0 (the results of your round call) can't be precisely represented in floating point.

A further complication is that there are tradeoffs between speed on the one hand and consistent, intuitive results on the other. For example, an Intel processor's FPU stores values in an 80-bit extended precision format, while a C/C++ double is typically 64 bits. For performance, a compiler may leave values in the FPU, as 80-bit temporaries, even though this can produce different results than what you'd get if you truncated them to 64 bits.

With your debug statement enabled, current is likely stored to memory, truncating it to 64 bits, which permits a direct comparison with last.

With the debug statement disabled, current is likely an 80-bit value stored in an FPU register, and thus it can never equal last, as long as last is a 64-bit value and they're both trying to store an inexact floating point representation of x/1000.0.

The solution is to use a floating point comparison with some allowed error (because a direct check for equality with floating point is almost never a good idea).

Further notes: I haven't looked at the assembly output to verify that this is the case; you can do this yourself if you want. I'm only able to reproduce the problem if I enable optimizations. You may be able to "fix" the bug by tweaking compiler flags to choose consistency over speed, but the correct solution is to use an inexact comparison instead of a direct check for equality.

Answer 2

Instead of

while(current != last);

use something like:

while(fabs(current - last) > tolerence);

where tolerance could be a small number, such as 1.0e-6.