比较 int = floor(sqrt(...)) 的 Fortran 和 C++ 汇编器

Question

我在 Fortran 和 C++ 中分别实现了一个函数：

#include <math.h>

void dbl_sqrt_c(double *x, double *y){
   *y = sqrt(*x - 1.0);
   return;
}

pure subroutine my_dbl_sqrt(x,y) bind(c, name="dbl_sqrt_fort")
   USE, INTRINSIC :: ISO_C_BINDING
   implicit none
   real(kind=c_double), intent(in)  :: x
   real(kind=c_double), intent(out) :: y

   y = sqrt(x - 1d0)
end subroutine my_dbl_sqrt

我在编译器资源管理器中比较了它们：

Fortran：https://godbolt.org/z/froz4rx97
C++：https://godbolt.org/z/45aex99Yz

按照我阅读汇编程序的方式，它们的作用基本相同，但 C++ 会检查 sqrt 的参数是否为负，而 Fortran 不会。我使用 googles benchmark 比较了它们的性能，但它们非常匹配：

--------------------------------------------------------
Benchmark              Time             CPU   Iterations
--------------------------------------------------------
bm_dbl_c/8          2.07 ns         2.07 ns    335965892
bm_dbl_fort/8       2.06 ns         2.06 ns    338643106

这是有趣的部分。如果我把它变成基于整数的函数：

void int_sqrt_c(int *x, int *y){
   *y = floor(sqrt(*x - 1.0));
   return;
}

和

pure subroutine my_int_root(x,y) bind(c, name="int_sqrt_fort")
   USE, INTRINSIC :: ISO_C_BINDING
   implicit none
   integer(kind=c_int), intent(in)  :: x
   integer(kind=c_int), intent(out) :: y

   y = floor(sqrt(x - 1d0))
end subroutine my_int_root

那么这就是他们开始分歧的地方：

--------------------------------------------------------
Benchmark              Time             CPU   Iterations
--------------------------------------------------------
bm_int_c/8          3.05 ns         3.05 ns    229239198
bm_int_fort/8       2.13 ns         2.13 ns    328933185

Fortran 代码似乎并没有因为这种变化而明显变慢，但 C++ 代码却减慢了 50%。这似乎相当大。这些是程序集：

Fortran：https://godbolt.org/z/axqqrc5E1
C++：https://godbolt.org/z/h7K75oKbn

Fortran 程序集看起来非常简单。它只是增加了double 和int 之间的转换，其他不多，但C++ 似乎做得更多，我不完全理解。

为什么 C++ 汇编器要复杂得多？如何改进 C++ 代码以实现匹配性能？

Answer 1

TL;DR：你被糟糕的默认值和与过时机器的兼容性所困扰：糟糕的默认值是 gcc 设置 errno 用于浮点计算（尽管 C 语言不需要这样做），以及与 x86 机器的兼容性没有比 SSE2 更好的 SSE 指令。如果您想生成体面的代码，请将-fno-math-errno -msse4 添加到compiler flags。

包含浮点硬件的现代处理器架构通常将平方根计算作为基本操作（指令）提供，如果平方根指令的操作数超出范围（负）。另一方面，旧的指令集架构可能没有浮点指令，或者没有硬件加速的平方根指令，因此 C 语言允许实现在超出范围的参数上设置 errno，但 @987654326 @ 作为线程本地内存位置实际上可以防止任何健全的架构直接从平方根指令设置 errno。为了获得不错的性能，gcc 通过调用硬件指令 (sqrtsd) 内联平方根计算，但要设置 errno，它会单独检查参数的符号，并仅在以防万一时调用库函数 参数是否定的，所以库函数可以设置errno。是的，这很疯狂，但这反过来又是课程的标准。您可以通过在编译器标志中设置 -fno-math-errno 来避免这种没人需要或想要的脑损伤。

较新的 x86-64 处理器具有比 AMD 首次开发的原始 x86-64 中更多的指令（​​其中仅包括 SSE2 向量/浮点指令）。添加的指令中有浮点/整数转换指令，允许受控舍入/截断，因此不必在软件中实现。您可以通过指定支持这些指令的目标来让 gcc 使用这些新指令，例如使用 -msse4 编译器标志。请注意，如果生成的程序在不支持这些指令的目标上运行，这将导致生成的程序出错，因此生成的程序的可移植性会降低（尽管它不会明显降低源代码的可移植性）。

Answer 2

正如@EOF 解释的那样，始终使用-fno-math-errno，如果可以，请使用-march=native，或者如果您不需要二进制文件在第一代Core 2等旧机器上运行，至少使用-msse4，或者AMD 飞鸿 II。 （最好还有其他好东西，比如popcnt和SSE4.2，比如
-march=nehalem -mtune=haswell）

最好是-O3，以便在可能的情况下启用 SIMD 自动矢量化，尽管有时 不可能 用于浮点，例如像点积或数组总和这样的减少，因为 FP 数学不是很相关。

(-ffast-math 会让编译器假装它是这样，或者您可以使用#pragma omp simd reduction(+:my_var_name) 授予它在一个循环中重新排列操作顺序的许可。但这通常是不安全的，例如，如果您的代码执行 Kahan 之类的操作求和以补偿 FP 舍入误差，并启用其他优化，例如将非规范化视为 0。）

---

您可以通过省略 floor() 从 gfortran / gcc 中获得更好的 asm。似乎都没有意识到 sqrt 的任何结果都是非负数或 NaN，因此 floor(sqrt) = trunc(sqrt)

（相关：如果要将双精度数舍入到最接近的 int，而不是 floor 或 trunc，请使用 lrint(x) 或 (int)nearbyint(x)，它们可以内联到使用当前 FP 舍入模式的指令，其中（除非你改变了它）将是四舍五入到最近的（即使作为决胜局）。有关它如何为 x86-64 和 AArch64 编译，请参阅 this answer。）

---

GCC 的后端没有优化底层（使用任何一种前端语言），-msse4 只是让它足够便宜，足以让sqrtsd 的吞吐量瓶颈隐藏其成本你的 C 基准测试。具体来说，我们得到

# gcc -O2 -fno-math-errno -msse4   with floor() still in the source.
        sqrtsd  xmm0, xmm0
        roundsd xmm0, xmm0, 9              # floor separately, result as a double
        cvttsd2si       eax, xmm0          # convert (with Truncation) to signed int

即使没有 SSE4，GFortran 也选择使用 floor-conversion 到 int（它只适用于 int 的范围，而不适用于该范围之外的双精度数，这是 GCC 的代码在没有 -msse4 的情况下手动执行的操作):

# GFortran -O2 -msse4      # with floor()   # chooses not to use SSE4
        sqrtsd  xmm0, xmm0
        cvttsd2si       eax, xmm0          # eax = int trunc(sqrt result)
        cvtsi2sd        xmm1, eax          # convert back to double
        ucomisd xmm0, xmm1                 # compare, setting CF if truncation rounded *up*
        sbb     eax, 0                     # eax -= CF

有趣的事实：此代码避免在ucomisd FLAGS 结果中检查“无序”(PF=1)。这种情况下 CF = 1，但显然 gfortran 并不关心它会使结果为 0x7FFFFFFF 而不是 0x80000000。

gcc 可以为 C 完成此操作，因为the behaviour is undefined if double->int 转换结果不适合 int。 （有趣的是，负双精度 -> 无符号也是 UB，因为这个原因；模块化范围缩减规则仅适用于宽整数类型 -> 无符号。）所以这是一个 gcc 错过了没有 SSE4 的 C 的优化错误。

当 SSE4 可用时，最好在转换前使用roundsd（在一般情况下，不能只优化地板）。到整数和返回的往返转换是12 cycle latency on Skylake，因此可能额外增加了 7 次转换，而只有一次转换，加上 ucomisd + sbb 延迟。与 roundsd 的 8 个周期相比。而roundsd 的总微指令更少 (https://uops.info)。所以这是对gfortran -msse4 的错过优化，它继续使用其 double->int->double compare / sbb 技巧进行下限转换。

---

优化源码：去掉floor

C（和 Fortran）FP -> int 转换默认截断（向0.0 舍入）。对于非负整数，这相当于floor，所以它是多余的。由于编译器没有意识到和/或没有利用 sqrt 结果为非负（或 NaN）这一事实，我们可以自己删除它：

#include <math.h>
int int_sqrt_c(int x){
   return sqrt(x - 1.0);
   //return floor(sqrt(x - 1.0));
}

我还简化了您的 C 以获取/返回一个 int，而不是通过指针。 https://godbolt.org/z/4zWeaej9T

int_sqrt_c(int):
        pxor    xmm0, xmm0
        cvtsi2sd        xmm0, edi             # int -> double
        subsd   xmm0, QWORD PTR .LC0[rip]
        sqrtsd  xmm0, xmm0
        cvttsd2si       eax, xmm0             # truncating double -> int
        ret

Fortan 代码的区别完全相同，https://godbolt.org/z/EhbTjhGen - 删除 floor 会删除无用的指令，只留下 cvttsd2si。

roundsd 在 Skylake 上花费 2 微秒并且有 8 个周期延迟（就像两个加法链）https://uops.info/，因此避免它从 x -> retval 和前面的总延迟中减少了相当一部分- 最终吞吐量成本。

（您的吞吐量基准测试将成为 sqrtsd 的瓶颈 on the back-end throughput，例如在 Skylake 上每 4.5 个周期一个，并且 OoO exec 隐藏了延迟，因此您当前的测试设置可能无法测量差异。）