clang/gcc 仅使用 -ffast-math 生成 fma；为什么？答案

【问题标题】：clang/gcc only generates fma with -ffast-math; why?clang/gcc 仅使用 -ffast-math 生成 fma；为什么？
【发布时间】：2019-09-22 06:04:27
【问题描述】：

在 icc 19 上，点积编译为 fma 指令上的循环。在 clang 和 gcc 上，fma 仅使用-ffast-math 生成。

但是，-ffast-math 违反了 IEEE 合规性，但 fma 完全符合 IEEE-754 2008，所以如果我必须使用-ffast-math 进行编译，则会导致其他问题。

为什么没有-ffast-math，gcc 和 clang 不生成 fma 指令？

Godbolt;编译器标志是-O3 -march=skylake-avx512, +- -ffast-math。

【问题讨论】：

您为什么认为 FMA 符合 IEEE-754？如果程序员写 a*b+c 表示执行 IEEE-754 乘法，然后执行 IEEE-754 加法，这与 fma(a, b, c) 不同，后者执行乘法后加法，就好像乘法具有无限精度一样。在某些情况下，这两件事会产生不同的结果。
@EricPostpischil：传统上，我将“打破 IEEE 合规性”与偷工减料联系在一起，以提高速度同时降低准确性。由于 FMA 是 IEEE-754 2008 中的标准化操作，它通过乘法和加法提高了精度，所以我想象它完全符合将 a*b+c 编译为 fma 的要求。但是你有没有说明这种解释是错误的参考资料？
IEEE-754 没有定义从编程语言到 IEEE-754 操作的绑定。这取决于语言。 C 的附录 F（IEC 60559 浮点算法，IEC 60559 实际上是 IEEE 754）将 * 和 + 绑定到 IEC 60559/IEEE 754 乘法和加法。所以a*b+c 必须是两个操作，一个乘法和一个单独的加法。它不能是融合乘加。不采用附件 F 的 C 实现倾向于遵循这一点。 GCC 和 Clang 似乎正在这样做。
@EricPostpischil：谢谢，这就是我想要的；想把它变成答案吗？
看看这个Godbolt。请注意，gcc 将标量函数编译为 fma，但未在内联它的地方使用 fma。我不确定答案（怀疑优化器的决定有问题），我只知道接受的答案是错误的。使用fma 是合法的取决于FP_CONTRACT。

【解决方案1】：

编译器是否适合使用融合乘法/加法将点积写实现为dot({a,c}, {b,d}) := a*b + c*d，给出 fl(?⋅? + fl(?⋅?)) 就好像它是写fma(a,b, c*d)？一般不会！

假设我们要评估平方差 ?² − ?²。

这可以写成点积dot({x,y}, {-x,y}) = x*x - y*y。当 ? ≈ ? 时，这种幼稚的公式会遭受灾难性的取消，但至少当 ? = ? 时，它可靠地返回零，因为 fl(fl(?²) − fl(?²)) = fl(fl(?²) − fl(?²) ) = fl(0) = 0。

这可以用 FMA 作为fma(x,x, -y*y) 来计算。但如果 ? = ? = fl(1.234) = 0x1.3be76c8b43958p+0，那么结果是 -1.3532e7b3d8ep−55 ≈ −3.352 × 10⁻¹⁷，在 IEEE 754 binary64 算术中，而不是我们希望的零。

它不仅非零，而且是负，所以如果你尝试在下游取平方根，即使你可以保证上游的? ≥ ?，你也会遇到 NaN。

（当然，分解(x + y)*(x - y) 可以更好地避免中间的灾难性取消，但这个问题是关于在没有额外假设的情况下评估点积。）
假设我们要在直角坐标中评估复数乘积 (? + ??)⋅(? + ??) = (?? − ??) + (?? + ??)?。

它的虚部可以写成点积dot({a,d}, {b,c}) = a*d + b*c。它可以用 FMA 作为fma(a,d, b*c) 来计算。

您可能期望复数 ? + ?? 与其复共轭 ? - ?? 的乘积是实数，虚部为零 - 如果使用 a*d + b*c 计算，它会如此，但如果使用 fma(a,d, b*c) 计算则不会。例如，如果? = fl(1.234) = 1.3be76c8b43958p+0 和 ? = fl(5.678) = 1.6b645a1cac083p+2，那么 fl(?⋅(−?) + fl(?⋅?)) = −1.6f6512a94ffp− 55 ≈ 3.983 × 10⁻¹⁷。

因此，在这些场景中使用 FMA 的编译器将是一种糟糕的形式，而您可以通过使用 <math.h> 中的 fma 函数编写 fma(a,b, c*d) 或添加 #pragma STDC FP_CONTRACT ON 来授权此类恶作剧。

也就是说…… 简单地通过-O2 -march=haswell 来persuade GCC 10.2 to abuse vfmadd231sd 似乎并不难，即使是显式的#pragma STDC FP_CONTRACT OFF 和same for ICC 21.1.9。这在我看来就像一个有缺陷的优化器！相比之下，Clang 11.0.1 uses vfmadd231sd with #pragma STDC FP_CONTRACT ON，但没有省略 pragma 或将其设置为 OFF。

【讨论】：

请注意，该问题专门针对点积。
@user14717 我改写了点产品；是不是更清楚了？
不，您已将自己限制为长度为 2 个向量。考虑（按照传统方式）长度为 n 的向量。
编译器不知道您永远不会使用您的点积来实现这些，并且您不关心 FMA 会违反的属性。如果您将x*y + a 重写为fma(x,y,a)，语义会发生实质性变化——这样做会导致出人意料的麻烦，这就是fma 以及FP_CONTRACT 和-ffast-math 等选项存在的原因。