了解为什么 ASM fsqrt 实现比标准 sqrt 函数更快
Understanding why an ASM fsqrt implementation is faster than the standard sqrt function
出于学术目的，我一直在研究 C++ 中的基本数学函数实现。今天，我对以下代码进行了平方根基准测试：
inline float sqrt_new(float n)
{
__asm {
fld n
fsqrt
}
}
我很惊讶地发现它始终比标准 sqrt 函数快（它大约占标准函数执行时间的 85%）。
我不太明白为什么，很想更好地理解它。下面我展示了我用来分析的完整代码（在 Visual Studio 2015 年，在发布模式下编译并打开所有优化）：
#include <iostream>
#include <random>
#include <chrono>
#define M 1000000
float ranfloats[M];
using namespace std;
inline float sqrt_new(float n)
{
__asm {
fld n
fsqrt
}
}
int main()
{
default_random_engine randomGenerator(time(0));
uniform_real_distribution<float> diceroll(0.0f , 1.0f);
chrono::high_resolution_clock::time_point start1, start2;
chrono::high_resolution_clock::time_point end1, end2;
float sqrt1 = 0;
float sqrt2 = 0;
for (int i = 0; i<M; i++) ranfloats[i] = diceroll(randomGenerator);
start1 = std::chrono::high_resolution_clock::now();
for (int i = 0; i<M; i++) sqrt1 += sqrt(ranfloats[i]);
end1 = std::chrono::high_resolution_clock::now();
start2 = std::chrono::high_resolution_clock::now();
for (int i = 0; i<M; i++) sqrt2 += sqrt_new(ranfloats[i]);
end2 = std::chrono::high_resolution_clock::now();
auto time1 = std::chrono::duration_cast<std::chrono::milliseconds>(end1 - start1).count();
auto time2 = std::chrono::duration_cast<std::chrono::milliseconds>(end2 - start2).count();
cout << "Time elapsed for SQRT1: " << time1 << " seconds" << endl;
cout << "Time elapsed for SQRT2: " << time2 << " seconds" << endl;
cout << "Average of Time for SQRT2 / Time for SQRT1: " << time2 / time1 << endl;
cout << "Equal to standard sqrt? " << (sqrt1 == sqrt2) << endl;
system("pause");
return 0;
}
编辑：我正在编辑问题以包括计算平方根的两个循环的反汇编代码，因为它们在 Visual Studio 2015 年出现。
首先，for (int i = 0; i<M; i++) sqrt1 += sqrt(ranfloats[i]);的反汇编：
00091194 0F 5A C0 cvtps2pd xmm0,xmm0
00091197 E8 F2 18 00 00 call __libm_sse2_sqrt_precise (092A8Eh)
0009119C F2 0F 5A C0 cvtsd2ss xmm0,xmm0
000911A0 83 C6 04 add esi,4
000911A3 F3 0F 58 44 24 4C addss xmm0,dword ptr [esp+4Ch]
000911A9 F3 0F 11 44 24 4C movss dword ptr [esp+4Ch],xmm0
000911AF 81 FE 90 5C 46 00 cmp esi,offset __dyn_tls_dtor_callback (0465C90h)
000911B5 7C D9 jl main+190h (091190h)
接下来，for (int i = 0; i<M; i++) sqrt2 += sqrt_new(ranfloats[i]);的反汇编：
00091290 F3 0F 10 00 movss xmm0,dword ptr [eax]
00091294 F3 0F 11 44 24 6C movss dword ptr [esp+6Ch],xmm0
0009129A D9 44 24 6C fld dword ptr [esp+6Ch]
0009129E D9 FA fsqrt
000912A0 D9 5C 24 6C fstp dword ptr [esp+6Ch]
000912A4 F3 0F 10 44 24 6C movss xmm0,dword ptr [esp+6Ch]
000912AA 83 C0 04 add eax,4
000912AD F3 0F 58 44 24 54 addss xmm0,dword ptr [esp+54h]
000912B3 F3 0F 11 44 24 54 movss dword ptr [esp+54h],xmm0
000912B9 ?? ?? ??
000912BA ?? ?? ??
000912BB ?? ?? ??
000912BC ?? ?? ??
000912BD ?? ?? ??
000912BE ?? ?? ??
000912BF ?? ?? ??
000912C0 ?? ?? ??
000912C1 ?? ?? ??
000912C2 ?? ?? ??
000912C3 ?? ?? ??
000912C4 ?? ?? ??
000912C5 ?? ?? ??
000912C6 ?? ?? ??
000912C7 ?? ?? ??
000912C8 ?? ?? ??
000912C9 ?? ?? ??
000912CA ?? ?? ??
000912CB ?? ?? ??
000912CC ?? ?? ??
000912CD ?? ?? ??
000912CE ?? ?? ??
000912CF ?? ?? ??
000912D0 ?? ?? ??
000912D1 ?? ?? ??
000912D2 ?? ?? ??
000912D3 ?? ?? ??
000912D4 ?? ?? ??
000912D5 ?? ?? ??
000912D6 ?? ?? ??
000912D7 ?? ?? ??
000912D8 ?? ?? ??
000912D9 ?? ?? ??
000912DA ?? ?? ??
000912DB ?? ?? ??
000912DC ?? ?? ??
000912DD ?? ?? ??
000912DE ?? ?? ??
compiling in Release mode and with all optimizations turned on
它们没有全部打开，你错过了一个。在 IDE 中是项目 > 属性 > C/C++ > 代码生成 > 浮点模型。您将其保留为默认设置 /fp:precise。这对生成的机器代码有非常明显的副作用：
00091197 E8 F2 18 00 00 call __libm_sse2_sqrt_precise (092A8Eh)
也许很直观，在 CRT 中调用辅助函数总是比像 FSQRT 这样的内联指令慢。
关于/fp 的确切语义有很多说法，关于它的MSDN 文章不是很好。逆向工程也很困难，微软从英特尔购买了代码，无法获得允许他们重新发布汇编代码的源代码许可。它最初的目标当然是处理由 Intel 的 8087 FPU 设计引起的可怕的浮点一致性问题。这在今天已经不那么重要了，所有主流的 C 和 C++ 编译器现在都发出 SSE2 代码。 MSVC++ 从 VS2012 开始就这样做了。这些英特尔库函数现在主要确保浮点运算仍然产生与旧版本编译器一致的结果。
__libm_sse2_sqrt_precise() 做的相当多。尝试记录未记录的函数的风险相当大，我想我看到了：
- 检查MXCSR寄存器的值，SSE的控制寄存器。如果它没有默认值 (0x1F80)，则代码假定程序员已使用 _control_fp() 并跳转到使用 FPU 语义计算的遗留 sqrt() 实现。
- 检查 FPU 控制寄存器的值以查看是否启用了浮点异常。通常关闭，如果有任何启用，则它会跳转到上面的 sqrt()。
- 检查参数是否小于零，但不为负 0，调用用户提供的
_matherr() 函数。
- 最终用SQRTSD指令计算结果。
None 实际上与精度有任何关系 :) 看到这个以 85% 的性能执行是一个很好的结果，但是 FSQRT 比 SQRTSD 慢得多。后者在现代处理器中得到了更多的硅爱。
如果您关心快速浮点运算，请将设置更改为 /fp:fast。产生：
00D91310 sqrtsd xmm0,xmm0
内联指令而不是库调用。换句话说，跳过前一个列表中的前 3 个项目符号。也轻松击败 FSQRT。
你的两个循环都非常糟糕，除了 sqrt 函数调用或 FSQRT 指令之外还有许多瓶颈。并且至少比最佳标量 SQRTSS (single-precision) 代码慢 2 倍。这可能比像样的 SSE2 矢量化循环所能达到的速度慢 8 倍。即使不重新排序任何数学运算，您也可以击败 SQRTSS 吞吐量。
https://gcc.gnu.org/wiki/DontUseInlineAsm 中的许多原因都适用于您的示例。编译器将无法通过您的函数传播常量，并且它不会知道结果总是 non-negative（如果它不是 NaN）。如果您稍后对数字进行平方，它也将无法将其优化为 fabs()。
同样非常重要的是，您使用 SSE2 SQRTPS (_mm_sqrt_ps()) 击败了 auto-vectorization。使用内部函数的 "no-error-checking" 标量 sqrt() 函数也存在该问题。 IDK 如果有任何方法可以在没有 /fp:fast 的情况下获得最佳结果，但我对此表示怀疑。 （除了在汇编中编写一个完整的循环，或者自己用内部函数向量化整个循环）。
令人印象深刻的是，您的 Haswell CPU 设法 运行 和 function-call 循环一样快，尽管 inline-asm 循环甚至可能不会饱和 FSQRT吞吐量。
出于某种原因，您的库函数调用正在调用 double sqrt(double)，而不是 C++ 重载 float sqrt(float)。这导致转换为 double 并返回 float。可能你需要 #include <cmath> to get the overloads，或者你可以调用 sqrtf()。 Linux 上的 gcc 和 clang 使用您当前的代码调用 sqrtf()（无需转换为 double 和返回），但也许它们的 <random> header 恰好包含 <cmath>，而 MSVC 的没有。或者可能还有其他事情发生。
库 function-call 循环 将总和保存在内存中（而不是寄存器）。显然，32 位版本的 __libm_sse2_sqrt_precise 使用的调用约定不保留任何 XMM 寄存器。 Windows x64 ABI 确实保留了 XMM6-XMM15，but wikipedia says this is new and the 32-bit ABI didn't do that。我假设如果有任何 call-preserved XMM 寄存器，MSVC 的优化器会利用它们。
无论如何，除了在每个独立标量浮点数上调用 sqrt 的吞吐量瓶颈外，loop-carried 对 sqrt1 的依赖是一个延迟瓶颈，其中包括 store-forwarding 往返：
000911A3 F3 0F 58 44 24 4C addss xmm0,dword ptr [esp+4Ch]
000911A9 F3 0F 11 44 24 4C movss dword ptr [esp+4Ch],xmm0
乱序执行让每次迭代的其余代码重叠，因此您只是吞吐量的瓶颈，但无论库 sqrt 函数的效率如何，这种延迟瓶颈将循环限制为每 6 + 3 次迭代一次= 9 个周期。 （Haswell ADDSS 延迟 = 3，XMM 的 store-forwarding 延迟 load/store = 6 个周期。整数寄存器比 store-forwarding 多 1 个周期。参见 Agner Fog's instruction tables。）
SQRTSD 的吞吐量为每 8-14 个周期一个，因此 loop-carried 依赖性不是 Haswell 的限制瓶颈。
inline-asm 版本 对于 sqrt 结果有一个 store/reload 往返，但它不是 loop-carried 依赖链的一部分. MSVC inline-asm syntax makes it hard to avoid store-forwarding round trips 将数据输入/输出内联汇编。但更糟糕的是，您在 x87 堆栈上生成结果，而编译器想要在 XMM 寄存器中进行 SSE 数学运算。
然后 MSVC 无缘无故搬起石头砸自己的脚，将总和保存在内存中而不是 XMM 寄存器中。它在 inline-asm 语句中查看它们影响了哪些寄存器，所以 IDK 为什么它看不到你的 inline-asm 语句没有破坏任何 XMM regs。
所以 MSVC 在这里做了很多比必要的更糟糕的工作：
00091290 movss xmm0,dword ptr [eax] # load from the array
00091294 movss dword ptr [esp+6Ch],xmm0 # store to the stack
0009129A fld dword ptr [esp+6Ch] # x87 load from stack
0009129E fsqrt
000912A0 fstp dword ptr [esp+6Ch] # x87 store to the stack
000912A4 movss xmm0,dword ptr [esp+6Ch] # SSE load from the stack (of sqrt(array[i]))
000912AA add eax,4
000912AD addss xmm0,dword ptr [esp+54h] # SSE load+add of the sum
000912B3 movss dword ptr [esp+54h],xmm0 # SSE store of the sum
因此它具有与 function-call 循环相同的 loop-carried 依赖链 (ADDSS + store-forwarding)。 Haswell FSQRT 每 8-17 个周期有一个吞吐量，所以它可能仍然是瓶颈。 （所有涉及数组值的 stores/reloads 对于每次迭代都是独立的，并且 out-of-order 执行可以重叠许多迭代以隐藏延迟链。但是，它们会阻塞 load/store 执行单元和有时会延迟 critical-path loads/stores 一个额外的周期。这称为资源冲突。）
如果没有 /fp:fast，如果结果为 NaN，sqrtf() 库函数必须设置 errno。这就是为什么它不能内联到一个 SQRTSS。
如果您确实想自己实现一个 no-checks 标量 sqrt 函数，您可以使用英特尔内在语法来实现：
// DON'T USE THIS, it defeats auto-vectorization
static inline
float sqrt_scalar(float x) {
__m128 xvec = _mm_set_ss(x);
xvec = _mm_cvtss_f32(_mm_sqrt_ss(xvec));
}
使用 gcc 和 clang（没有 -ffast-math）编译成 near-optimal 标量循环。在 the Godbolt compiler explorer 上查看：
# gcc6.2 -O3 for the sqrt_new loop using _mm_sqrt_ss. good scalar code, but don't optimize further.
.L10:
movss xmm0, DWORD PTR [r12]
add r12, 4
sqrtss xmm0, xmm0
addss xmm1, xmm0
cmp r12, rbx
jne .L10
这个循环应该只在 SQRTSS 吞吐量上出现瓶颈（在 Haswell 上每 7 个时钟一个，明显快于 SQRTSD 或 FSQRT），并且没有资源冲突。然而，即使没有 re-ordering FP 添加 （因为 FP add/mul aren't truly associative），与你可以做的相比，它仍然是垃圾：一个聪明的编译器（或使用内在函数的程序员）会使用 SQRTPS 获得 4 个结果，其吞吐量与来自 SQRTSS 的 1 个结果相同。将 SQRT 结果的向量解压缩为 4 个标量，然后您可以保持完全相同的运算顺序，并对中间结果进行相同的舍入。我很失望 clang 和 gcc 没有这样做
然而，gcc 和 clang 确实设法避免调用库函数。 clang3.9（只有 -O3）使用 SQRTSS，甚至不检查 NaN。我认为这是合法的，而不是编译器错误。也许它看到代码没有使用 errno？
另一方面，gcc6.2 推测内联 sqrtf()，使用 SQRTSS 和检查输入以查看是否需要调用库函数。
# gcc6.2's sqrt() loop, without -ffast-math.
# speculative inlining of SQRTSS with a check + fallback
# spills/reloads a lot of stuff to memory even when it skips the call :(
# xmm1 = 0.0 (gcc -fverbose-asm says it's holding sqrt2, which is zero-initialized, so I guess gcc decides to reuse that zero)
.L9:
movss xmm0, DWORD PTR [rbx]
sqrtss xmm5, xmm0
ucomiss xmm1, xmm0 # compare input against 0.0
movss DWORD PTR [rsp+8], xmm5
jbe .L8 # if(0.0 <= SQRTSS input || unordered(0.0, input)) { skip the function call; }
movss DWORD PTR [rsp+12], xmm1 # silly gcc, this store isn't needed. ucomiss doesn't modify xmm1
call sqrtf # called for negative inputs, but not for NaN.
movss xmm1, DWORD PTR [rsp+12]
.L8:
movss xmm4, DWORD PTR [rsp+4] # silly gcc always stores/reloads both, instead of putting the stores/reloads inside the block that the jbe skips
addss xmm4, DWORD PTR [rsp+8]
add rbx, 4
movss DWORD PTR [rsp+4], xmm4
cmp rbp, rbx
jne .L9
gcc 不幸地搬起石头砸自己的脚，就像 MSVC 处理 inline-asm 一样：有一个 store-forwarding 往返作为 loop-carried 依赖项。所有 spill/reloads 都可能在 JBE 跳过的块内。也许 gcc 东西负输入会很常见。
更糟糕的是，如果您使用 /fp:fast 或 -ffast-math，即使像 clang 这样聪明的编译器也无法将您的 _mm_sqrt_ss 重写为 SQRTPS。 Clang 通常不仅擅长将内在函数映射到指令 1:1，而且如果您错过组合事物的机会，它会提出更优化的洗牌和混合。
所以在启用快速 FP 数学的情况下，使用 _mm_sqrt_ss 是一个很大的损失 。 clang 将 sqrt() 库函数调用版本编译成 RSQRTPS + newton-raphson 迭代。
另请注意，您的微基准测试代码对 sqrt_new() 实施的延迟不敏感，仅对吞吐量敏感。延迟在实际 FP 代码中通常很重要，而不仅仅是吞吐量。但在其他情况下，比如对许多数组元素独立地做同样的事情，延迟并不重要，因为 out-of-order 执行可以通过从许多循环迭代中运行指令来很好地隐藏它。
正如我之前提到的，extra store/reload round trip your data takes on its way in/out of MSVC-style inline-asm 的延迟在这里是一个严重的问题。当 MSVC 内联函数时，fld n 不直接来自数组。
顺便说一句，Skylake 的吞吐量为每 3 个周期一个 SQRTPS/SS，但仍有 12 个周期的延迟。 SQRTPD/SD 吞吐量 = 每 4-6 个周期一个，延迟 = 15-16 个周期。所以 FP 平方根在 Skylake 上比在 Haswell 上更流水线化。这放大了基准 FP sqrt 延迟与吞吐量之间的差异。
出于学术目的，我一直在研究 C++ 中的基本数学函数实现。今天，我对以下代码进行了平方根基准测试：
inline float sqrt_new(float n)
{
__asm {
fld n
fsqrt
}
}
我很惊讶地发现它始终比标准 sqrt 函数快（它大约占标准函数执行时间的 85%）。
我不太明白为什么，很想更好地理解它。下面我展示了我用来分析的完整代码（在 Visual Studio 2015 年，在发布模式下编译并打开所有优化）：
#include <iostream>
#include <random>
#include <chrono>
#define M 1000000
float ranfloats[M];
using namespace std;
inline float sqrt_new(float n)
{
__asm {
fld n
fsqrt
}
}
int main()
{
default_random_engine randomGenerator(time(0));
uniform_real_distribution<float> diceroll(0.0f , 1.0f);
chrono::high_resolution_clock::time_point start1, start2;
chrono::high_resolution_clock::time_point end1, end2;
float sqrt1 = 0;
float sqrt2 = 0;
for (int i = 0; i<M; i++) ranfloats[i] = diceroll(randomGenerator);
start1 = std::chrono::high_resolution_clock::now();
for (int i = 0; i<M; i++) sqrt1 += sqrt(ranfloats[i]);
end1 = std::chrono::high_resolution_clock::now();
start2 = std::chrono::high_resolution_clock::now();
for (int i = 0; i<M; i++) sqrt2 += sqrt_new(ranfloats[i]);
end2 = std::chrono::high_resolution_clock::now();
auto time1 = std::chrono::duration_cast<std::chrono::milliseconds>(end1 - start1).count();
auto time2 = std::chrono::duration_cast<std::chrono::milliseconds>(end2 - start2).count();
cout << "Time elapsed for SQRT1: " << time1 << " seconds" << endl;
cout << "Time elapsed for SQRT2: " << time2 << " seconds" << endl;
cout << "Average of Time for SQRT2 / Time for SQRT1: " << time2 / time1 << endl;
cout << "Equal to standard sqrt? " << (sqrt1 == sqrt2) << endl;
system("pause");
return 0;
}
编辑：我正在编辑问题以包括计算平方根的两个循环的反汇编代码，因为它们在 Visual Studio 2015 年出现。
首先，for (int i = 0; i<M; i++) sqrt1 += sqrt(ranfloats[i]);的反汇编：
00091194 0F 5A C0 cvtps2pd xmm0,xmm0
00091197 E8 F2 18 00 00 call __libm_sse2_sqrt_precise (092A8Eh)
0009119C F2 0F 5A C0 cvtsd2ss xmm0,xmm0
000911A0 83 C6 04 add esi,4
000911A3 F3 0F 58 44 24 4C addss xmm0,dword ptr [esp+4Ch]
000911A9 F3 0F 11 44 24 4C movss dword ptr [esp+4Ch],xmm0
000911AF 81 FE 90 5C 46 00 cmp esi,offset __dyn_tls_dtor_callback (0465C90h)
000911B5 7C D9 jl main+190h (091190h)
接下来，for (int i = 0; i<M; i++) sqrt2 += sqrt_new(ranfloats[i]);的反汇编：
00091290 F3 0F 10 00 movss xmm0,dword ptr [eax]
00091294 F3 0F 11 44 24 6C movss dword ptr [esp+6Ch],xmm0
0009129A D9 44 24 6C fld dword ptr [esp+6Ch]
0009129E D9 FA fsqrt
000912A0 D9 5C 24 6C fstp dword ptr [esp+6Ch]
000912A4 F3 0F 10 44 24 6C movss xmm0,dword ptr [esp+6Ch]
000912AA 83 C0 04 add eax,4
000912AD F3 0F 58 44 24 54 addss xmm0,dword ptr [esp+54h]
000912B3 F3 0F 11 44 24 54 movss dword ptr [esp+54h],xmm0
000912B9 ?? ?? ??
000912BA ?? ?? ??
000912BB ?? ?? ??
000912BC ?? ?? ??
000912BD ?? ?? ??
000912BE ?? ?? ??
000912BF ?? ?? ??
000912C0 ?? ?? ??
000912C1 ?? ?? ??
000912C2 ?? ?? ??
000912C3 ?? ?? ??
000912C4 ?? ?? ??
000912C5 ?? ?? ??
000912C6 ?? ?? ??
000912C7 ?? ?? ??
000912C8 ?? ?? ??
000912C9 ?? ?? ??
000912CA ?? ?? ??
000912CB ?? ?? ??
000912CC ?? ?? ??
000912CD ?? ?? ??
000912CE ?? ?? ??
000912CF ?? ?? ??
000912D0 ?? ?? ??
000912D1 ?? ?? ??
000912D2 ?? ?? ??
000912D3 ?? ?? ??
000912D4 ?? ?? ??
000912D5 ?? ?? ??
000912D6 ?? ?? ??
000912D7 ?? ?? ??
000912D8 ?? ?? ??
000912D9 ?? ?? ??
000912DA ?? ?? ??
000912DB ?? ?? ??
000912DC ?? ?? ??
000912DD ?? ?? ??
000912DE ?? ?? ??
compiling in Release mode and with all optimizations turned on
它们没有全部打开，你错过了一个。在 IDE 中是项目 > 属性 > C/C++ > 代码生成 > 浮点模型。您将其保留为默认设置 /fp:precise。这对生成的机器代码有非常明显的副作用：
00091197 E8 F2 18 00 00 call __libm_sse2_sqrt_precise (092A8Eh)
也许很直观，在 CRT 中调用辅助函数总是比像 FSQRT 这样的内联指令慢。
关于/fp 的确切语义有很多说法，关于它的MSDN 文章不是很好。逆向工程也很困难，微软从英特尔购买了代码，无法获得允许他们重新发布汇编代码的源代码许可。它最初的目标当然是处理由 Intel 的 8087 FPU 设计引起的可怕的浮点一致性问题。这在今天已经不那么重要了，所有主流的 C 和 C++ 编译器现在都发出 SSE2 代码。 MSVC++ 从 VS2012 开始就这样做了。这些英特尔库函数现在主要确保浮点运算仍然产生与旧版本编译器一致的结果。
__libm_sse2_sqrt_precise() 做的相当多。尝试记录未记录的函数的风险相当大，我想我看到了：
- 检查MXCSR寄存器的值，SSE的控制寄存器。如果它没有默认值 (0x1F80)，则代码假定程序员已使用 _control_fp() 并跳转到使用 FPU 语义计算的遗留 sqrt() 实现。
- 检查 FPU 控制寄存器的值以查看是否启用了浮点异常。通常关闭，如果有任何启用，则它会跳转到上面的 sqrt()。
- 检查参数是否小于零，但不为负 0，调用用户提供的
_matherr()函数。 - 最终用SQRTSD指令计算结果。
None 实际上与精度有任何关系 :) 看到这个以 85% 的性能执行是一个很好的结果，但是 FSQRT 比 SQRTSD 慢得多。后者在现代处理器中得到了更多的硅爱。
如果您关心快速浮点运算，请将设置更改为 /fp:fast。产生：
00D91310 sqrtsd xmm0,xmm0
内联指令而不是库调用。换句话说，跳过前一个列表中的前 3 个项目符号。也轻松击败 FSQRT。
你的两个循环都非常糟糕，除了 sqrt 函数调用或 FSQRT 指令之外还有许多瓶颈。并且至少比最佳标量 SQRTSS (single-precision) 代码慢 2 倍。这可能比像样的 SSE2 矢量化循环所能达到的速度慢 8 倍。即使不重新排序任何数学运算，您也可以击败 SQRTSS 吞吐量。
https://gcc.gnu.org/wiki/DontUseInlineAsm 中的许多原因都适用于您的示例。编译器将无法通过您的函数传播常量，并且它不会知道结果总是 non-negative（如果它不是 NaN）。如果您稍后对数字进行平方，它也将无法将其优化为 fabs()。
同样非常重要的是，您使用 SSE2 SQRTPS (_mm_sqrt_ps()) 击败了 auto-vectorization。使用内部函数的 "no-error-checking" 标量 sqrt() 函数也存在该问题。 IDK 如果有任何方法可以在没有 /fp:fast 的情况下获得最佳结果，但我对此表示怀疑。 （除了在汇编中编写一个完整的循环，或者自己用内部函数向量化整个循环）。
令人印象深刻的是，您的 Haswell CPU 设法 运行 和 function-call 循环一样快，尽管 inline-asm 循环甚至可能不会饱和 FSQRT吞吐量。
出于某种原因，您的库函数调用正在调用 double sqrt(double)，而不是 C++ 重载 float sqrt(float)。这导致转换为 double 并返回 float。可能你需要 #include <cmath> to get the overloads，或者你可以调用 sqrtf()。 Linux 上的 gcc 和 clang 使用您当前的代码调用 sqrtf()（无需转换为 double 和返回），但也许它们的 <random> header 恰好包含 <cmath>，而 MSVC 的没有。或者可能还有其他事情发生。
库 function-call 循环 将总和保存在内存中（而不是寄存器）。显然，32 位版本的 __libm_sse2_sqrt_precise 使用的调用约定不保留任何 XMM 寄存器。 Windows x64 ABI 确实保留了 XMM6-XMM15，but wikipedia says this is new and the 32-bit ABI didn't do that。我假设如果有任何 call-preserved XMM 寄存器，MSVC 的优化器会利用它们。
无论如何，除了在每个独立标量浮点数上调用 sqrt 的吞吐量瓶颈外，loop-carried 对 sqrt1 的依赖是一个延迟瓶颈，其中包括 store-forwarding 往返：
000911A3 F3 0F 58 44 24 4C addss xmm0,dword ptr [esp+4Ch]
000911A9 F3 0F 11 44 24 4C movss dword ptr [esp+4Ch],xmm0
乱序执行让每次迭代的其余代码重叠，因此您只是吞吐量的瓶颈，但无论库 sqrt 函数的效率如何，这种延迟瓶颈将循环限制为每 6 + 3 次迭代一次= 9 个周期。 （Haswell ADDSS 延迟 = 3，XMM 的 store-forwarding 延迟 load/store = 6 个周期。整数寄存器比 store-forwarding 多 1 个周期。参见 Agner Fog's instruction tables。）
SQRTSD 的吞吐量为每 8-14 个周期一个，因此 loop-carried 依赖性不是 Haswell 的限制瓶颈。
inline-asm 版本 对于 sqrt 结果有一个 store/reload 往返，但它不是 loop-carried 依赖链的一部分. MSVC inline-asm syntax makes it hard to avoid store-forwarding round trips 将数据输入/输出内联汇编。但更糟糕的是，您在 x87 堆栈上生成结果，而编译器想要在 XMM 寄存器中进行 SSE 数学运算。
然后 MSVC 无缘无故搬起石头砸自己的脚，将总和保存在内存中而不是 XMM 寄存器中。它在 inline-asm 语句中查看它们影响了哪些寄存器，所以 IDK 为什么它看不到你的 inline-asm 语句没有破坏任何 XMM regs。
所以 MSVC 在这里做了很多比必要的更糟糕的工作：
00091290 movss xmm0,dword ptr [eax] # load from the array
00091294 movss dword ptr [esp+6Ch],xmm0 # store to the stack
0009129A fld dword ptr [esp+6Ch] # x87 load from stack
0009129E fsqrt
000912A0 fstp dword ptr [esp+6Ch] # x87 store to the stack
000912A4 movss xmm0,dword ptr [esp+6Ch] # SSE load from the stack (of sqrt(array[i]))
000912AA add eax,4
000912AD addss xmm0,dword ptr [esp+54h] # SSE load+add of the sum
000912B3 movss dword ptr [esp+54h],xmm0 # SSE store of the sum
因此它具有与 function-call 循环相同的 loop-carried 依赖链 (ADDSS + store-forwarding)。 Haswell FSQRT 每 8-17 个周期有一个吞吐量，所以它可能仍然是瓶颈。 （所有涉及数组值的 stores/reloads 对于每次迭代都是独立的，并且 out-of-order 执行可以重叠许多迭代以隐藏延迟链。但是，它们会阻塞 load/store 执行单元和有时会延迟 critical-path loads/stores 一个额外的周期。这称为资源冲突。）
如果没有 /fp:fast，如果结果为 NaN，sqrtf() 库函数必须设置 errno。这就是为什么它不能内联到一个 SQRTSS。
如果您确实想自己实现一个 no-checks 标量 sqrt 函数，您可以使用英特尔内在语法来实现：
// DON'T USE THIS, it defeats auto-vectorization
static inline
float sqrt_scalar(float x) {
__m128 xvec = _mm_set_ss(x);
xvec = _mm_cvtss_f32(_mm_sqrt_ss(xvec));
}
使用 gcc 和 clang（没有 -ffast-math）编译成 near-optimal 标量循环。在 the Godbolt compiler explorer 上查看：
# gcc6.2 -O3 for the sqrt_new loop using _mm_sqrt_ss. good scalar code, but don't optimize further.
.L10:
movss xmm0, DWORD PTR [r12]
add r12, 4
sqrtss xmm0, xmm0
addss xmm1, xmm0
cmp r12, rbx
jne .L10
这个循环应该只在 SQRTSS 吞吐量上出现瓶颈（在 Haswell 上每 7 个时钟一个，明显快于 SQRTSD 或 FSQRT），并且没有资源冲突。然而，即使没有 re-ordering FP 添加 （因为 FP add/mul aren't truly associative），与你可以做的相比，它仍然是垃圾：一个聪明的编译器（或使用内在函数的程序员）会使用 SQRTPS 获得 4 个结果，其吞吐量与来自 SQRTSS 的 1 个结果相同。将 SQRT 结果的向量解压缩为 4 个标量，然后您可以保持完全相同的运算顺序，并对中间结果进行相同的舍入。我很失望 clang 和 gcc 没有这样做
然而，gcc 和 clang 确实设法避免调用库函数。 clang3.9（只有 -O3）使用 SQRTSS，甚至不检查 NaN。我认为这是合法的，而不是编译器错误。也许它看到代码没有使用 errno？
另一方面，gcc6.2 推测内联 sqrtf()，使用 SQRTSS 和检查输入以查看是否需要调用库函数。
# gcc6.2's sqrt() loop, without -ffast-math.
# speculative inlining of SQRTSS with a check + fallback
# spills/reloads a lot of stuff to memory even when it skips the call :(
# xmm1 = 0.0 (gcc -fverbose-asm says it's holding sqrt2, which is zero-initialized, so I guess gcc decides to reuse that zero)
.L9:
movss xmm0, DWORD PTR [rbx]
sqrtss xmm5, xmm0
ucomiss xmm1, xmm0 # compare input against 0.0
movss DWORD PTR [rsp+8], xmm5
jbe .L8 # if(0.0 <= SQRTSS input || unordered(0.0, input)) { skip the function call; }
movss DWORD PTR [rsp+12], xmm1 # silly gcc, this store isn't needed. ucomiss doesn't modify xmm1
call sqrtf # called for negative inputs, but not for NaN.
movss xmm1, DWORD PTR [rsp+12]
.L8:
movss xmm4, DWORD PTR [rsp+4] # silly gcc always stores/reloads both, instead of putting the stores/reloads inside the block that the jbe skips
addss xmm4, DWORD PTR [rsp+8]
add rbx, 4
movss DWORD PTR [rsp+4], xmm4
cmp rbp, rbx
jne .L9
gcc 不幸地搬起石头砸自己的脚，就像 MSVC 处理 inline-asm 一样：有一个 store-forwarding 往返作为 loop-carried 依赖项。所有 spill/reloads 都可能在 JBE 跳过的块内。也许 gcc 东西负输入会很常见。
更糟糕的是，如果您使用 /fp:fast 或 -ffast-math，即使像 clang 这样聪明的编译器也无法将您的 _mm_sqrt_ss 重写为 SQRTPS。 Clang 通常不仅擅长将内在函数映射到指令 1:1，而且如果您错过组合事物的机会，它会提出更优化的洗牌和混合。
所以在启用快速 FP 数学的情况下，使用 _mm_sqrt_ss 是一个很大的损失 。 clang 将 sqrt() 库函数调用版本编译成 RSQRTPS + newton-raphson 迭代。
另请注意，您的微基准测试代码对 sqrt_new() 实施的延迟不敏感，仅对吞吐量敏感。延迟在实际 FP 代码中通常很重要，而不仅仅是吞吐量。但在其他情况下，比如对许多数组元素独立地做同样的事情，延迟并不重要，因为 out-of-order 执行可以通过从许多循环迭代中运行指令来很好地隐藏它。
正如我之前提到的，extra store/reload round trip your data takes on its way in/out of MSVC-style inline-asm 的延迟在这里是一个严重的问题。当 MSVC 内联函数时，fld n 不直接来自数组。
顺便说一句，Skylake 的吞吐量为每 3 个周期一个 SQRTPS/SS，但仍有 12 个周期的延迟。 SQRTPD/SD 吞吐量 = 每 4-6 个周期一个，延迟 = 15-16 个周期。所以 FP 平方根在 Skylake 上比在 Haswell 上更流水线化。这放大了基准 FP sqrt 延迟与吞吐量之间的差异。