4

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

#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())，所以这将是一个“部分解决方法”，而不是一个解决方案......我也想避免像“不要使用浮点”这样的解决方案（;

4

3 回答 3

4

一般不会：这正是您使用的价格（巧合的是，威廉·卡汉（William Kahan）在自动收缩问题中指出的-ffp-contract=fast正是这个例子）

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

#pragma STDC FP_CONTRACT OFF
// non-contracted code
#pragma STDC FP_CONTRACT ON
于 2017-02-09T10:31:38.690 回答
3

有趣的是，多亏了 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）？

于 2017-02-09T17:44:57.010 回答
1

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

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

如果mx并且my按照评论建议的那样计算，那么即使最终mx*mx + my*my使用 FMA 完成，它仍然会导致预期的 0 结果。问题是，由于没有任何其他变量使用sx//乘法、加法和减法以一步计算，然后可以在机器代码中以任意数量的不同方式表示（以任何顺序，可能使用多个 FMA 等），但它认为它将获得最佳性能大计算。syexeymxmydistance

但是，如果其他东西（如打印语句）引用mxand my，那么编译器更有可能在distance作为第二步计算之前单独计算它们。在这种情况下，数学确实按照评论建议的方式计算，即使最终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

于 2019-02-21T00:52:48.790 回答