比较Fortran和C++汇编程序的int=floor(sqrt(…))

html5 • • 问答

我在 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              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

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

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

和

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++ 代码以实现匹配性能？

回答

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

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

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

回答

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

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

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

通过省略floor(). 似乎都没有意识到 sqrt 的任何结果都是非负或 NaN，因此floor(sqrt) = trunc(sqrt)

（相关：如果要将 double 舍入到最近的 int，而不是 floor 或 trunc，请使用lrint(x)or (int)nearbyint(x)，它可以内联到使用当前 FP 舍入模式的指令，该模式（除非您更改它）将舍入到-nearest （即使作为 tiebreak）。有关它如何为 x86-64 和 AArch64 编译，请参阅此答案。）

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 也选择使用地板转换到 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

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

gcc可以为 C 这样做，因为如果double->int 转换结果不适合 int ，则行为未定义。（有趣的事实，负双 -> unsigned 也是 UB 出于这个原因；模块化范围缩减规则仅适用于宽整数类型 -> unsigned。）所以这是一个 gcc 错过了 C 没有 SSE4 的优化错误。

当 SSE4 可用时，最好roundsd在转换前使用（在不能直接优化地板的一般情况下）。在 Skylake 上往返转换为整数和返回是12 个周期的延迟，所以额外的可能是 7 次，而不是一次转换，加上 ucomisd + sbb 延迟。与 8 个周期相比roundsd。而且roundsd总 uops ( 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 uops 并且有 8 个周期延迟（就像两个添加的链）https://uops.info/，因此避免它从 x -> retval 和前面的总延迟中取出相当一部分 -最终吞吐量成本。

（你的吞吐量将基准瓶颈的背部末端吞吐量的sqrtsd，每4.5个周期上SKYLAKE微架构OOO EXEC皮的延迟例如一，和，所以你当前的测试设置可能会无法衡量的差异。）

以上是比较Fortran和C++汇编程序的int=floor(sqrt(…))的全部内容。
THE END
分享
二维码
< <上一篇
下一篇>>