昨天我在跟踪我的项目中的一个错误，几个小时后，我已经缩小到一段或多或少在做这样的事情的代码：
#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次减法运算(我希望得到两个零)，然后对它们进行平方运算(再次得到两个0)并将它们相加(零)。
事实证明，问题的根本原因是使用了融合的乘加运算，这使得结果不精确(从我的角度来看)。一般来说，我并不反对这种优化，因为它承诺会给出更精确的结果，但在这种情况下，1.34925e-06与我预期的0相去甚远。
测试用例是非常"脆弱"的——如果您启用更多的print或更多的assert，它就会停止断言，因为编译器不再使用融合乘加。例如，如果我取消注释所有行：
$ 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，因为编译器库内部无论如何都会使用融合乘法-加法(我在sinf()和cosf()中看到了很多)，所以这将是一个"部分解决方案"，而不是一个解决方案。。。我还想避免像"不要使用浮点"(；
一般来说，否：这正是您使用-ffp-contract=fast所付出的代价(巧合的是，William Kahan在自动收缩问题中指出的正是这个例子)
从理论上讲，如果您使用的是C(而不是C++)，并且您的编译器支持C-1999杂注(即不是gcc)，那么您可以使用
#pragma STDC FP_CONTRACT OFF
// non-contracted code
#pragma STDC FP_CONTRACT ON
有趣的是，由于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)？
我可以看到这个问题已经存在了一段时间，但如果其他人在寻找答案时遇到它，我想我会提到几点。。
首先，如果不分析生成的汇编代码，很难准确地判断，但我怀疑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);
(ε的确切值可能取决于应用程序，当然，对于不同类型的计算，甚至可能不同)
同样，对于浮点数，你应该说if (abs(a - b) < EPSILON)等，而不是if (a == b)。
减少(但不一定消除)这个问题的另一种方法是在应用程序中实现"快速故障"逻辑。例如，在上面的代码中，您可以通过在计算distance之前测试if (mx < EPSILON && my < EPSILON)来"短路"一些数学运算，如果它们都为零，则跳过其余的运算(因为您知道在这种情况下结果将为零)，而不是一直计算distance，然后看看它最后是否为0。你越快发现这种情况，错误积累的机会就越少(有时你也可以避免在不需要的情况下进行一些更昂贵的计算)。