Универсальный способ обработки неточностей с плавающей запятой, умножения и сложения

Вчера я отслеживал ошибку в своем проекте, которая — через несколько часов — я сузился до фрагмента кода, который более или менее делал что-то вроде этого:

#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);
}

После компиляции и исполнения:

$ 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)

С моей точки зрения, что-то не так, поскольку я попросил 2 вычитания двух поразрядно-идентичных пар (я ожидал получить два нуля), затем возводил их в квадрат (снова два нуля) и складывал их вместе (ноль).

Оказывается, что основной причиной проблемы является использование операции слияния-умножения-добавления, которая где-то вдоль линии делает результат неточным (с моей точки зрения). Вообще я ничего не имею против этой оптимизации, так как она обещает дать результаты, которые Больше точный, но в этом случае 1.34925e-06 действительно далек от 0, который я ожидал.

Тестовый пример очень «хрупкий» — если вы включаете больше распечаток или больше утверждений, он перестает утверждать, потому что компилятор больше не использует fused-multiply-add. Например, если я раскомментирую все строки:

$ 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

Поскольку я считал, что это ошибка в компиляторе, я сообщил об этом, но он был закрыт объяснением, что это правильное поведение.

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

Поэтому мне интересно — как следует кодировать такие вычисления, чтобы избежать проблемы? Я думал о универсальном решении, но что-то лучше, чем:

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

Я хотел бы исправить или улучшить свой код — если есть что исправить / улучшить — вместо настройки -ffp-contract=off для всего моего проекта, поскольку fused-multiply-add используется внутренне в библиотеках компилятора (я вижу много этого в sinf () и cosf ()), так что это будет «частичный обход», а не решение … Я также хотел бы избежать решений, таких как «не использовать с плавающей точкой» (;