Generic way of handling fused-multiply-add floating-point inaccuracies

Question

Yesterday I was tracking a bug in my project, which - after several hours - I've narrowed down to a piece of code which more or less was doing something like this:

#include <iostream>
#include <cmath>
#include <cassert>

volatile float r = -0.979541123;
volatile float alpha = 0.375402451;

int main()
{
    float sx = r * cosf(alpha); // -0.911326
    float sy = r * sinf(alpha); // -0.359146
    float ex = r * cosf(alpha); // -0.911326
    float ey = r * sinf(alpha); // -0.359146
    float mx = ex - sx;     // should be 0
    float my = ey - sy;     // should be 0
    float distance = sqrtf(mx * mx + my * my) * 57.2958f;   // should be 0, gives 1.34925e-06

//  std::cout << "sv: {" << sx << ", " << sy << "}" << std::endl;
//  std::cout << "ev: {" << ex << ", " << ey << "}" << std::endl;
//  std::cout << "mv: {" << mx << ", " << my << "}" << std::endl;
    std::cout << "distance: " << distance << std::endl;

    assert(distance == 0.f);
//  assert(sx == ex && sy == ey);
//  assert(mx == 0.f && my == 0.f);
} 

After compilation and execution:

$ g++ -Wall -Wextra -Wshadow -march=native -O2 vfma.cpp && ./a.out 
distance: 1.34925e-06
a.out: vfma.cpp:23: int main(): Assertion `distance == 0.f' failed.
Aborted (core dumped)

From my point of view something is wrong, as I've asked for 2 subtractions of two bitwise-identical pairs (I expected to get two zeroes), then squaring them (two zeroes again) and adding them together (zero).

It turns out that the root cause of problem is the use of fused-multiply-add operation, which somewhere along the line makes the result inexact (from my point of view). Generally I have nothing against this optimization, as it promises to give results which are more exact, but in this case 1.34925e-06 is really far from the 0 that I was expecting.

The test case is very "fragile" - if you enable more prints or more asserts, it stops asserting, because compiler doesn't use fused-multiply-add anymore. For example if I uncomment all lines:

$ g++ -Wall -Wextra -Wshadow -march=native -O2 vfma.cpp && ./a.out 
sv: {-0.911326, -0.359146}
ev: {-0.911326, -0.359146}
mv: {0, 0}
distance: 0

As I've considered this to be a bug in the compiler, I've reported that, but it got closed with the explanation that this is correct behaviour.

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=79436

So I'm wondering - how should one code such calculations to avoid the problem? I was thinking about a generic solution, but something better than:

mx = ex != sx ? ex - sx : 0.f;

I would like to fix or improve my code - if there's anything to fix/improve - instead of setting -ffp-contract=off for my whole project, as fused-multiply-add is used internally in the compiler libraries anyway (I see a lot of that in sinf() and cosf()), so it would be a "partial work-around", not a solution... I would also like to avoid solutions like "don't use floating-point" (;

Answer 1

In general no: this is exactly the price you pay for using -ffp-contract=fast (coincidently, it is precisely this example that William Kahan notes in the problems with automatic contraction)

Theoretically, if you were using C (not C++), and your compiler supported C-1999 pragmas (i.e. not gcc), you could use

#pragma STDC FP_CONTRACT OFF
// non-contracted code
#pragma STDC FP_CONTRACT ON