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

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 次减法（我希望得到两个零），然后将它们平方（再次两个零）并将它们加在一起（零） .

事实证明，问题的根本原因是使用了 fused-multiply-add 操作，这使得结果不准确（从我的角度来看）。一般来说，我不反对这种优化，因为它承诺提供更准确的结果，但在这种情况下，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()) 中有很多，所以这将是一个“部分解决方法”，而不是一个解决方案......我也想避免像“不要使用浮点”这样的解决方案（;

Answer 1

通常不会：这正是您使用-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

Answer 2

有趣的是，多亏了 fma，浮点数 mx 和 my 为您提供了将 r 和 cos 相乘时产生的舍入误差。

fma( r,cos, -r*cos) = theoretical(r*cos) - float(r*cos)

因此，由于浮点数相乘（但不考虑 cos 和 sin 计算中的舍入误差），您得到的结果以某种方式表明计算出的 (sx,sy) 与理论 (sx,sy) 相差多远。

所以问题是您的程序如何依赖与浮点舍入相关的不确定区间内的差异（ex-sx,ey-sy）？

Answer 3

我可以看到这个问题已经存在了一段时间，但如果其他人在寻找答案时遇到它，我想我会提到几点..

首先，如果不分析生成的汇编代码，很难准确判断，但我怀疑 FMA 给出的结果远远超出预期的原因不仅仅是 FMA 本身，还在于您假设所有计算是按照您指定的顺序进行的，但在优化 C/C++ 编译器时，情况通常并非如此。这也可能是取消注释打印语句会改变结果的原因。

如果 mx 和 my 按照 cmets 的建议进行计算，那么即使最终的 mx*mx + my*my 是使用 FMA 完成的，它仍然会导致预期的 0 结果。问题是，由于没有任何 sx/sy/ex/ey/mx/my 变量被其他任何东西使用，因此编译器很可能从未实际评估它们作为独立变量，只需将所有数学运算组合成大量的乘法、加法和减法，即可一步计算出distance，然后可以用机器代码以多种不同的方式表示（在任何顺序，可能有多个 FMA 等），但它认为它会为这一大计算获得最佳性能。

但是，如果其他内容（如打印语句）引用 mx 和 my，则编译器更有可能在第二步计算 distance 之前单独计算它们。在这种情况下，数学确实按照 cmets 建议的方式进行计算，即使是最终 distance 计算中的 FMA 也不会改变结果（因为输入全都是 0）。

答案

但这实际上并不能回答真正的问题。为了回答这个问题，一般来说，避免此类问题的最稳健（并且通常推荐）的方法是：永远不要假设浮点运算会永远产生一个精确的数字，即使该数字是 0。这意味着，一般来说，使用== 比较浮点数是个坏主意。相反，您应该选择一个较小的数字（通常称为 epsilon），它大于任何可能/可能的累积误差，但仍小于任何显着的结果（例如，如果您知道您关心的距离只是真的显着到小数点后几位，那么您可以选择EPSILON = 0.01，这意味着“任何小于 0.01 的差异我们都将视为与零相同”）。然后，而不是说：

assert(distance == 0.f);

你会说：

assert(distance < EPSILON);

（您的 epsilon 的确切值可能取决于应用程序，当然，对于不同类型的计算甚至可能会有所不同）

同样，对于浮点数，不要说if (a == b) 之类的东西，而是说if (abs(a - b) < EPSILON) 之类的东西。

减少（但不一定消除）此问题的另一种方法是在您的应用程序中实现“快速失败”逻辑。例如，在上面的代码中，与其一路计算distance，然后看看它最后是否为0，不如在你得到之前通过测试if (mx < EPSILON && my < EPSILON)“短路”一些数学。到计算distance 并在它们都为零时跳过其余部分（因为您知道在这种情况下结果将为零）。越快掌握情况，累积错误的机会就越少（有时您也可以避免在不需要的情况下进行一些成本更高的计算）。