比较 int = floor(sqrt(...)) 的 Fortran 和 C++ 汇编程序
Comparing Fortran & C++ assembler for int = floor(sqrt(...))
我分别用 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 则不会。我使用谷歌基准比较了它们的性能，但它们相当匹配：
--------------------------------------------------------
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++ 代码以实现匹配性能？
TL;DR：你被错误的默认值和与过时机器的兼容性所困扰：错误的默认值是用于浮点计算的 gcc 设置 errno（尽管 C 语言不需要这样做），并且与没有比 SSE2 更好的 SSE 指令的 x86 机器的兼容性。如果您想要体面的代码生成，请将 -fno-math-errno -msse4 添加到 compiler flags.
包含浮点硬件的现代处理器架构通常将平方根计算作为原始操作（指令）提供，如果平方根指令的操作数超出范围，它会在浮点环境中设置错误标志（消极的）。另一方面，旧的指令集体系结构可能没有浮点指令，或者没有硬件加速平方根指令，因此 C 语言允许实现在超出范围的参数上设置 errno，但是errno 作为线程本地内存位置实际上阻止了任何理智的体系结构直接从平方根指令设置 errno 。为了获得不错的性能，gcc 通过调用硬件指令 (sqrtsd) 内联平方根计算，但要设置 errno，它会单独检查参数的符号，并调用库函数 仅在参数为负数的情况下，所以库函数可以设置errno。是的，这很疯狂，但这反过来又是理所当然的。您可以通过在编译器标志中设置 -fno-math-errno 来避免这种没有人需要或想要的脑损伤。
相当新的 x86-64 处理器比最初由 AMD 开发的原始 x86-64 中的指令更多（仅包含 SSE2 vector/floating-point 指令）。添加的指令中有 float/integer 允许受控 rounding/truncation 的转换指令，因此这不必在软件中实现。您可以通过指定支持这些指令的目标来让 gcc 使用这些新指令，例如通过使用 -msse4 编译器标志。请注意，如果在不支持这些指令的目标上 运行 ，这将导致生成的程序出错，因此生成的程序将不那么可移植（尽管它不会明显降低源代码的可移植性).
正如@EOF 所解释的那样，始终使用 -fno-math-errno，如果可以的话使用 -march=native，或者如果不需要 运行 的二进制文件，至少使用 -msse4 ] 在第一代 Core 2 或 AMD Phenom II 等旧机器上。 （最好还有其他好东西，比如 popcnt 和 SSE4.2，比如
-march=nehalem -mtune=haswell)
最好是 -O3 以在可能的情况下启用 SIMD 自动矢量化，尽管有时 not 可能用于浮点数，例如减少点积或数组的总和，因为 FP 数学不是很相关。
(-ffast-math 会让编译器假装它是，或者您可以使用 #pragma omp simd reduction(+:my_var_name) 授予它在一个循环中重新排列操作顺序的许可。但这通常不安全，例如，如果你的代码会做 Kahan 求和之类的事情来补偿 FP 舍入误差，并启用其他优化，比如将非正规数视为 0。）
您可以通过省略 floor() 从 gfortran / gcc 中获得更好的 asm。他们似乎都没有意识到 sqrt 的任何结果都将是非负数或 NaN，因此 floor(sqrt) = trunc(sqrt)
（相关：如果你想将 double 舍入到最接近的 int，而不是 floor 或 t运行c，请使用 lrint(x) 或 (int)nearbyint(x)，它们可以内联到使用当前 FP 舍入模式的指令，该模式（除非您更改它）将舍入到最接近的值（即使是平分）。有关它如何针对 x86-64 和 AArch64 进行编译，请参阅 this answer。）
GCC 的后端不优化客场（使用任何一种前端语言），-msse4 只是让它足够便宜以应对[的吞吐量瓶颈=27=] 在你的 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 也选择使用 flooring-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（在一般情况下不能直接优化 floor）。整数和返回整数的往返转换是 12 cycle latency on Skylake, so an extra maybe 7 vs. just one convert, plus ucomisd + sbb latency. vs. 8 cycles for roundsd. And roundsd is fewer total uops (https://uops.info)。所以这是对 gfortran -msse4 的一个错过的优化，它继续使用其 double->int->double compare / sbb 技巧进行 floor-conversion。
优化源：删除 floor
默认情况下，C（和 Fortran）FP -> int 转换 t运行cates（向 0.0 舍入）。对于非负整数，这相当于 floor，所以它是多余的。由于编译器没有意识到 and/or 没有利用 sqrt 结果是非负（或 NaN）这一事实，我们可以自己删除它：
#include <math.h>
int int_sqrt_c(int x){
return sqrt(x - 1.0);
//return floor(sqrt(x - 1.0));
}
我还将你的 C 简化为 take/return 一个 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 隐藏了延迟，因此您当前的测试设置可能无法衡量差异。）
我分别用 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 则不会。我使用谷歌基准比较了它们的性能，但它们相当匹配：
--------------------------------------------------------
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++ 代码以实现匹配性能？
TL;DR：你被错误的默认值和与过时机器的兼容性所困扰：错误的默认值是用于浮点计算的 gcc 设置 errno（尽管 C 语言不需要这样做），并且与没有比 SSE2 更好的 SSE 指令的 x86 机器的兼容性。如果您想要体面的代码生成，请将 -fno-math-errno -msse4 添加到 compiler flags.
包含浮点硬件的现代处理器架构通常将平方根计算作为原始操作（指令）提供，如果平方根指令的操作数超出范围，它会在浮点环境中设置错误标志（消极的）。另一方面，旧的指令集体系结构可能没有浮点指令，或者没有硬件加速平方根指令，因此 C 语言允许实现在超出范围的参数上设置 errno，但是errno 作为线程本地内存位置实际上阻止了任何理智的体系结构直接从平方根指令设置 errno 。为了获得不错的性能，gcc 通过调用硬件指令 (sqrtsd) 内联平方根计算，但要设置 errno，它会单独检查参数的符号，并调用库函数 仅在参数为负数的情况下，所以库函数可以设置errno。是的，这很疯狂，但这反过来又是理所当然的。您可以通过在编译器标志中设置 -fno-math-errno 来避免这种没有人需要或想要的脑损伤。
相当新的 x86-64 处理器比最初由 AMD 开发的原始 x86-64 中的指令更多（仅包含 SSE2 vector/floating-point 指令）。添加的指令中有 float/integer 允许受控 rounding/truncation 的转换指令，因此这不必在软件中实现。您可以通过指定支持这些指令的目标来让 gcc 使用这些新指令，例如通过使用 -msse4 编译器标志。请注意，如果在不支持这些指令的目标上 运行 ，这将导致生成的程序出错，因此生成的程序将不那么可移植（尽管它不会明显降低源代码的可移植性).
正如@EOF 所解释的那样，始终使用 -fno-math-errno，如果可以的话使用 -march=native，或者如果不需要 运行 的二进制文件，至少使用 -msse4 ] 在第一代 Core 2 或 AMD Phenom II 等旧机器上。 （最好还有其他好东西，比如 popcnt 和 SSE4.2，比如
-march=nehalem -mtune=haswell)
最好是 -O3 以在可能的情况下启用 SIMD 自动矢量化，尽管有时 not 可能用于浮点数，例如减少点积或数组的总和，因为 FP 数学不是很相关。
(-ffast-math 会让编译器假装它是，或者您可以使用 #pragma omp simd reduction(+:my_var_name) 授予它在一个循环中重新排列操作顺序的许可。但这通常不安全，例如，如果你的代码会做 Kahan 求和之类的事情来补偿 FP 舍入误差，并启用其他优化，比如将非正规数视为 0。）
您可以通过省略 floor() 从 gfortran / gcc 中获得更好的 asm。他们似乎都没有意识到 sqrt 的任何结果都将是非负数或 NaN，因此 floor(sqrt) = trunc(sqrt)
（相关：如果你想将 double 舍入到最接近的 int，而不是 floor 或 t运行c，请使用 lrint(x) 或 (int)nearbyint(x)，它们可以内联到使用当前 FP 舍入模式的指令，该模式（除非您更改它）将舍入到最接近的值（即使是平分）。有关它如何针对 x86-64 和 AArch64 进行编译，请参阅 this answer。）
GCC 的后端不优化客场（使用任何一种前端语言），-msse4 只是让它足够便宜以应对[的吞吐量瓶颈=27=] 在你的 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 也选择使用 flooring-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（在一般情况下不能直接优化 floor）。整数和返回整数的往返转换是 12 cycle latency on Skylake, so an extra maybe 7 vs. just one convert, plus ucomisd + sbb latency. vs. 8 cycles for roundsd. And roundsd is fewer total uops (https://uops.info)。所以这是对 gfortran -msse4 的一个错过的优化，它继续使用其 double->int->double compare / sbb 技巧进行 floor-conversion。
优化源：删除 floor
默认情况下，C（和 Fortran）FP -> int 转换 t运行cates（向 0.0 舍入）。对于非负整数，这相当于 floor，所以它是多余的。由于编译器没有意识到 and/or 没有利用 sqrt 结果是非负（或 NaN）这一事实，我们可以自己删除它：
#include <math.h>
int int_sqrt_c(int x){
return sqrt(x - 1.0);
//return floor(sqrt(x - 1.0));
}
我还将你的 C 简化为 take/return 一个 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 隐藏了延迟，因此您当前的测试设置可能无法衡量差异。）