floating-point - clang/gcc 只用 -ffast-math 生成 fma；为什么？

Question

在 icc 19 上，点积通过 fma 指令编译为循环。在 clang 和 gcc 上，fma 仅由 -ffast-math 生成.

然而，-ffast-math违反 IEEE 合规性，但 fma 完全符合 IEEE-754 2008，所以如果我必须使用 -ffast-math 进行编译，然后我会引起其他问题。

为什么 gcc 和 clang 在没有 -ffast-math 的情况下不生成 fma 指令?

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

Answer 1

编译器实现点积写成 dot({a,c}, {b,d}) := a*b + c*d 是否合适？使用融合乘法/加法，给出 fl(𝑎⋅𝑏 + fl(𝑐⋅𝑑)) 就像它写成 fma(a,b, c*d) ?
一般不会!
这里有几个改编自 lecture notes by W. Kahan on IEEE 754 的例子:

假设我们要计算平方差 𝑥² − 𝑦²。
这可以写成点积 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 中，我们可能希望二进制 6 算术为零
它不仅非零，而且是负数，因此如果您尝试在下游取平方根，即使您可以保证上游 𝑥 ≥ 𝑦，您也会遇到 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 的编译器形式会很糟糕，而您没有通过编写 fma(a,b, c*d) 来明确要求它。与 fma函数来自 <math.h> ，或添加 #pragma STDC FP_CONTRACT ON授权这种恶作剧。

那说…
看来不难 persuade GCC 10.2 to abuse vfmadd231sd只需通过 -O2 -march=haswell ，即使有明确的 #pragma STDC FP_CONTRACT OFF , 和 same for ICC 21.1.9 .
这在我看来就像一个有问题的优化器!
相比之下， Clang 11.0.1 uses vfmadd231sd with #pragma STDC FP_CONTRACT ON ，但没有省略编译指示或设置为 OFF .