Lun*_*din 9 c floating-point optimization clang
我正在执行这个问题的代码：为什么以下代码的输出不为零？
#include <stdio.h>
int main (void)
{
double A = 373737.0;
double B;
B = A*A*A + 0.37/A - A*A*A - 0.37/A;
printf("The value of B is %f.\n", B);
}
每个主流 x86 编译器的每个优化设置都会给出输出-0.000001。当我使用当前的 clang 15.0.0 时，我-O0也得到了这一点。
但是，使用 14.0.0 版本以上的 clang 和-O1to进行编译-O3会给出输出-1.000001。为什么会发生这种情况？这是一个已知的错误？
Godbolt 为了您的方便：https ://godbolt.org/z/M5j3fGhWf
这是一个相当深的兔子洞，我不知道我是否已经探索过它所有的曲折。但这是答案的初稿；欢迎提出改进建议。
\n从本质上讲，罪魁祸首是所谓的“融合乘加”（或者在本例中为融合乘减）。融合乘加是一条a*b+c单步计算的指令。这可以显着加速某些计算（例如使用霍纳规则的点积和多项式）。大约在2013年被添加到Intel的x86指令集中（Haswell）；一年前，AMD 芯片中添加了类似的指令。但这个想法并不新鲜。至少从 1990 年起，高端处理器就包含了此类指令（使用 IBM 的 POWER1 处理器）。
由于融合运算的结果仅一次（而不是在乘法后舍入两次并在加法后再次舍入），因此通常会产生更准确的结果。不幸的是，在某些情况下它会产生不太准确的结果，这就是其中之一；它是由计算a*b-cwherea*b和c非常相似而触发的，并且c之前已四舍五入。[注 1] 要查看问题的实际情况，将代码减少到最少是很有用的，其结果至少令人惊讶：
#include <stdio.h>\nint main (void) {\n double A = 373737.0;\n printf("A*A*A - A*A*A is %f.\\n", A*A*A - A*A*A);\n return 0;\n}\n从 v14.0.0 开始使用 clang，会打印出 1.000000。[注 2] 结果为 1（而不是 -1），因为该表达式被转换为、和A*A*A - A*A*A的融合乘减。现在，373737\xc2\xb3 正好是 52203339425426553，一个 56 位数字。由于在 x86 平台上只允许 53 个有效位，因此需要四舍五入到最接近的可表示值，即 52203339425426552。在融合操作中，精确计算 373737\xc2\xb2 * 373737，然后四舍五入的值 373737\ xc2\xb3 被减去，剩下 1。A*AAA*A*Adouble
在原始程序中，计算结果为（大约）373737\xc2\xb3 + 1e-6 - 373737\xc2\xb3 - 1e-6。在此计算中，首先计算（使用 FMA）并舍入 373737\xc2\xb3 + 1e-6，这又是 52203339425426552；添加 1e-6 对舍入总和没有影响。然后执行融合求反乘加，将 52203339425426552 与 373737\xc2\xb2 和 373737 的精确求反乘积相加 (-52203339425426553)；结果正是-1. 最后减去 1e-6，得到 的观测结果-1.000001。
这就是戈德堡所说的“灾难性取消”的本质（如果您还没有读过，请参阅注释 1）；两个非常相似的值相减就抵消了所有的意义。
\n（另一方面，您可以小心地利用融合运算中的乘法未舍入的事实，以产生更准确的最终结果，使用加拿大数学家 William Kahan的算法，他是 IEEE 的主要架构师-754 标准。例如，请参阅@njuffa 的启发性答案，了解当 b\xc2\xb2 接近 4ac 时如何准确计算二次根。）
\n那么 Clang v14.0.0 发生了什么变化？Clang 和 GCC 都有一个控制是否使用 FMA 的选项：-ffp-contract。（在 C 标准中，FMA 是“契约操作”的示例之一，该选项控制所有此类操作。）该选项具有三个可能的值：off、on和fast。off始终意味着编译器在编译表达式时不会融合乘法和加法。（如果该操作码在目标机器上可用，它仍然会将fma函数编译为 FMA 操作码。）直到 v13.0.0，off这是 Clang 的默认设置；在 v14.0.0 中，默认值更改为on，这允许在同一表达式中融合乘法和加法。从那时起，如果目标架构实现了 FMA 指令，Clang 将默认发出它们。与这个问题更相关的是，它还将模拟 FMA 以在编译时执行常量计算。
虽然GCC有相同的选项，但语义有些不同。据我所知，GCC 不会模拟 FMA 进行编译时计算。此外，GCC 解释-ffp-contract=on为与 (!) 相同-ffp-contract=off，其默认值为-ffp-contract=fast。该fast设置不仅允许在表达式内（标准 C 允许）进行收缩操作，还允许在跨不同表达式的计算中进行收缩操作。然而，对于这个特定的计算，GCC 的优化器更喜欢保存和重用公共子表达式的值A*A*A，而不是发出 FMA。[注3]
Clang 还允许-ffp-contract=fast，其语义与 GCC 大致相同，但指定该选项的结果是常量文件夹无法模拟 FMA。[注4]
C 标准实际上定义了一种可移植机制来控制契约操作的使用：#pragma STDC FP_CONTRACT, 以及可能的值ON,OFF和DEFAULT。OFF要求抑制FMA作业的排放，但标准没有其他限制；默认值可以是ON和OFF，并且ON不需要做任何特别的事情。然而，GCC 没有实现这个编译指示（从 GCC v12 开始），因此它不像人们希望的那样可移植。（不过，Clang 确实实现了。）
尽管正如这个问题所示，使用融合乘加可能会产生令人惊讶的结果，并且很容易陷入假设此类结果是编译器错误的陷阱，但很明显该标准确实打算编译器可以自由使用 FMA 和其他契约操作，只要有一种方法可以关闭该功能，如 \xc2\xa76.5 第 8 段所示，其措辞自 C99 以来没有改变：
\n\n\n浮动表达式可以被收缩，即，像单个操作一样进行计算，从而省略源代码和表达式计算方法隐含的舍入误差。
\nFP_CONTRACTin 的编译指示提供<math.h>了一种禁止缩写表达式的方法。否则，表达式是否以及如何收缩是由实现定义的。
该条款附有以下脚注：
\n\n\n该许可证专门用于允许实现利用组合多个 C 运算符的快速机器指令。由于收缩可能会破坏可预测性，甚至可能降低包含表达式的准确性，因此需要对其使用进行明确定义和清晰记录。
\n
有人认为附录 F 中关于 IEC-559 合规性的要求（通常描述为 IEEE-754/854）优先于上面明确提到的许可证，但我认为这个论点没有说服力。首先，如上所述，\xc2\xa76.5 非常清楚。其次，附录 F 还考虑了 \xc2\xa7F.7 中的收缩表达式：
\n\n\n收缩表达式被正确舍入（一次），并以与 IEC 60559 涵盖的基本算术运算一致的方式处理无穷大、NaN、有符号零、次正规和舍入方向。
\n
第三，IEEE-754（2008，注释 5）明确允许实现实现契约操作，只要它们提供一种将其关闭的方法：
\n\n\n语言标准应要求默认情况下，当未启用优化且未启用备用异常处理时，语言实现保留源代码的字面含义。\n\xe2\x80\xa6\n语言标准还应定义，并且要求实现提供允许和禁止对块单独或共同进行值更改优化的属性。这些优化可能包括但不限于：
\n\n
\n- 应用结合律或分配律。
\n- \n
fusedMultiplyAdd operation由乘法和加法合成 a 。\n\xe2\x80\xa6
我说这一切时带着一定的痛苦，因为我也很确定这种行为是有问题的。FMA应用的不可预测性似乎不太理想。另一方面，该标准定义了该fma函数，该函数应该（并且通常确实）被内联编译成适当的机器指令，并且有一些机制要求编译器不要发出收缩表达式，除非明确要求，我\'我肯定会考虑更加一致地使用。
这是大卫·戈德堡（David Goldberg）在《每个计算机科学家应该了解浮点运算的知识》一文中描述的“灾难性取消”的场景，任何有关浮点怪癖的讨论都不可避免地会引用该文章。戈德堡所说的“取消”是指通过减法取消有效数字，可能只留下误差范围内的数字。
\n至少，如果您指定了正确的编译器选项。使用默认编译器选项，您将得到 0。
\n正如OP中所述，默认编译器设置不会出现奇怪的结果。那是因为默认情况下没有优化。启用任何优化后，Clang 将在编译时折叠常量表达式，并且常量文件夹模拟融合乘加。如果不进行优化，计算将在运行时完成，并且默认情况下，Clang 不会发出 FMA 指令，因为它们并非在所有受支持的 x86 芯片上可用。您需要指定-mfma（或其他类似的目标选择器）来指示目标架构包含 FMA 指令集，以便在编译的二进制文件中查看 FMA 指令。
我不知道GCC的常量文件夹是否模拟FMA；如果我稍后弄清楚，我会编辑这一段。
\nLLVM 提交者 Andy Kaylor 在bug 54927 的评论-ffp-contract=fast中解释了在常量文件夹中抑制 FMA的原因。
我没有更高版本的副本，但我怀疑本质没有改变。
\n