performance - 了解为什么ASM fsqrt实现比标准sqrt函数更快

Question

为了学术目的，我在玩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 ??                   ?? ?? 

Answer 1

您的两个循环都非常可怕，除了sqrt函数调用或FSQRT指令之外，还有许多瓶颈。并且比最佳标量SQRTSS（单精度）代码至少慢2倍。这可能比像样的SSE2矢量化循环所能达到的速度慢8倍。即使不重新排序任何数学运算，您也可以击败SQRTSS吞吐量。

https://gcc.gnu.org/wiki/DontUseInlineAsm中的许多原因都适用于您的示例。编译器将无法通过函数传播常量，并且它将不知道结果总是非负的（如果不是NaN）。如果您以后将数字平方，它也将无法将其优化为fabs()。

同样非常重要的是，您要使用SSE2 SQRTPS（_mm_sqrt_ps()）击败自动矢量化。使用内在函数的“无错误检查”标量sqrt（）函数也存在该问题。 IDK是否可以在没有/fp:fast的情况下获得最佳结果，但我对此表示怀疑。 （除了在汇编中编写整个循环，或者自己使用内部函数向量化整个循环）。



您的Haswell CPU设法以最快的速度运行函数调用循环，这令人印象深刻，尽管inline-asm循环甚至可能也不会使FSQRT吞吐量饱和。

由于某种原因，您的库函数调用正在调用double sqrt(double)，而不是C ++重载float sqrt(float)。这导致转换为双精度，然后又返回为浮点型。可能您需要#include <cmath> to get the overloads，或者可以调用sqrtf()。 Linux上的gcc和clang用您当前的代码调用sqrtf（）（无需转换为double和back），但是也许它们的<random>标头恰好包含了<cmath>，而MSVC却没有。也许还有其他事情正在发生。



库函数调用循环将总和保留在内存中（而不是寄存器中）。显然，__libm_sse2_sqrt_precise的32位版本使用的调用约定不会保留任何XMM寄存器。 Windows x64 ABI确实保留XMM6-XMM15 but wikipedia says this is new and the 32-bit ABI didn't do that。我假设如果有任何保留呼叫的XMM寄存器，则MSVC的优化器将利用它们。

无论如何，除了在每个独立的标量浮点数上调用sqrt的吞吐量瓶颈之外，对sqrt1的循环承载依赖性是一个延迟瓶颈，其中包括存储转发往返：

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加载/存储的存储转发延迟= 6个周期。比整数寄存器的存储转发多1个周期。请参见 Agner Fog's instruction tables。）

SQRTSD的吞吐量为每8-14个循环1个，因此循环携带的依赖性不是Haswell的限制瓶颈。



带有inline-asm版本的sqrt结果具有存储/重载往返行程，但它不是循环承载的依赖链的一部分。 MSVC inline-asm syntax makes it hard to avoid store-forwarding round trips使数据进入/退出内联汇编。但更糟糕的是，您在x87堆栈上生成结果，并且编译器希望在XMM寄存器中进行SSE数学运算。

然后MSVC无缘无故地将自己砸死，将总和保留在内存中，而不是XMM寄存器中。它在inline-asm语句中查找以查看它们影响哪些寄存器，因此IDK为什么看不到inline-asm语句不会破坏任何XMM规则。

因此，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

因此，它具有与函数调用循环相同的循环承载依赖链（ADDSS +存储转发）。 Haswell FSQRT每8-17个周期就有一个吞吐量，因此可能仍然是瓶颈。 （所有涉及数组值的存储/重载对于每次迭代都是独立的，乱序执行可能会使许多迭代重叠以隐藏该延迟链。但是，它们会阻塞加载/存储执行单元，有时会延迟关键-path会以额外的周期加载/存储。这称为资源冲突。）



如果没有 /fp:fast，则如果结果为NaN，则 sqrtf()库函数必须设置 errno。这就是为什么它不能仅内联到SQRTSS。

如果您确实想自己实现无检查标量sqrt函数，则可以使用Intel内部函数语法来实现：

// 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）的近似最佳标量循环。在 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快）上成为瓶颈，并且不会出现资源冲突。但是，与即使不对FP添加重新排序（因为 FP add/mul aren't truly associative）也可以执行的操作相比，它仍然是垃圾：智能编译器（或使用内在函数的程序员）将使用SQRTPS获得4个结果，而吞吐量与1个结果相同SQRTSS。将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中所做的一样：存在存储转发往返作为循环携带的依赖项。所有溢出/重新装料可能在JBE跳过的块内。也许gcc负输入会很常见。



更糟糕的是，如果您确实使用 /fp:fast或 -ffast-math，即使像clang这样的聪明编译器也无法将 _mm_sqrt_ss重写为SQRTPS。 Clang通常不仅在将内在函数映射到指令1：1方面非常擅长，而且如果您错过了将事物组合的机会，它还会提供更优化的混洗和混合。

因此，启用快速FP数学运算后，使用 _mm_sqrt_ss是一个很大的损失。 clang将 sqrt()库函数调用版本编译为RSQRTPS +牛顿拉普森迭代。



还要注意，您的微基准代码对 sqrt_new()实现的延迟不敏感，而对吞吐量不敏感。延迟通常在真实的FP代码中很重要，而不仅仅是吞吐量。但是在其他情况下，就像独立地对许多数组元素执行相同的操作一样，延迟并不重要，因为无序执行可以通过让许多循环迭代中的指令在运行中将其隐藏得足够好。

如前所述， extra store/reload round trip your data takes on its way in/out of MSVC-style inline-asm的延迟在这里是一个严重的问题。当MSVC内联函数时， fld n并非直接来自数组。



顺便说一句，Skylake的SQRTPS / SS吞吐量为每3个周期1个，但仍为12个周期的延迟。 SQRTPD / SD吞吐量=每4-6个周期1个，延迟= 15-16个周期。因此，相对于Haswell，FP平方根在Skylake上的流水线更多。这放大了基准测试FP sqrt延迟与吞吐量之间的差异。