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

问题描述 投票：4回答：3

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

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

从我的角度来看，有些事情是错误的，因为我要求两个按位相同的对减去两个（我希望得到两个零），然后将它们平方（再两个零）并将它们加在一起（零）。

事实证明，问题的根本原因是使用了融合乘法 - 加法运算，沿线的某处使得结果不精确（从我的观点来看）。一般来说，我没有反对这种优化，因为它承诺给出更精确的结果，但在这种情况下，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，因为在编译器库中内部使用了fusion-multiply-add（我在sinf中看到了很多内容） （）和cosf（）），所以它将是一个“部分解决方案”，而不是解决方案......我也想避免像“不要使用浮点”这样的解决方案（;

c++ floating-point precision floating-accuracy fma
3个回答
3
投票

一般来说没有：这正是你使用-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
3
投票

有趣的是，多亏了fma，浮点数mx和my给出了乘以r和cos时产生的舍入误差。

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

因此，你得到的结果表明，由于浮点数的乘法（但不考虑cos和sin的计算中的舍入误差），计算出的（sx，sy）离理论值（sx，sy）有多远（sx，sy）。

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

1
投票

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

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

如果按照评论的建议计算mx和my，那么即使最终的mx*mx + my*my是用FMA完成的，它仍然会产生预期的0结果。问题在于，由于其他任何变量都没有使用sx / sy / ex / ey / mx / my变量，因此编译器很可能永远不会将它们作为自变量进行实际评估，而只是简单地计算所有数学一起进行大量的乘法，加法和减法，一步计算distance，然后可以用机器码中的任意数量的不同方式表示（以任何顺序，可能有多个FMA等）但是它表示它会为这一大计算获得最佳性能。

但是，如果其他东西（如print语句）引用mx和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并跳过其余部分（因为你知道在这种情况下结果将为零）。抓住这种情况越快，错误积累的机会越少（有时你也可以避免在不需要的情况下进行更昂贵的计算）。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.