c++ - 处理融合乘加浮点误差的通用方法

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

#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 相差甚远。

测试用例非常“脆弱” - 如果您启用更多打印或更多断言，它会停止断言，因为编译器不再使用 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())，所以这将是一个“部分解决方法”，而不是一个解决方案......我也想避免像“不要使用浮点”这样的解决方案(;

最佳答案

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

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

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

关于c++ - 处理融合乘加浮点误差的通用方法，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/42132260/