c++ - 处理融合乘加浮点不准确的通用方法

Question

昨天我在跟踪我的项目中的一个错误，几个小时后，我已经缩小到一段代码，它或多或少地在做这样的事情:

#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 次减法(我希望得到两个零)，然后对它们进行平方(再次是两个零)并将它们加在一起(零) .

事实证明，问题的根本原因是使用了融合乘加运算，这使得结果不准确(从我的角度来看)。一般来说，我并不反对这种优化，因为它 promise 会给出更准确的结果，但在这种情况下，1.34925e-06 与我预期的 0 相去甚远。

测试用例非常“脆弱” - 如果您启用更多打印或更多断言，它会停止断言，因为编译器不再使用融合乘法加法。例如，如果我取消注释所有行:

$ 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() 中看到很多这样的东西)，所以这将是一个“部分解决方法”，而不是解决方案......我也想避免像“不要使用 float "(;

Answer 1

通常不会:这正是您使用 -ffp-contract=fast 所付出的代价(巧合的是，正是这个例子 William Kahan notes in the problems with automatic contraction )

理论上，如果您使用的是 C(不是 C++)，并且您的编译器支持 C-1999 pragma(即不是 gcc)，您可以使用

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