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

5

我一直在为了学术目的而玩弄C++中基本的数学函数实现。今天,我对以下代码进行了平方根的基准测试:

inline float sqrt_new(float n)
{
    __asm {
        fld n
            fsqrt
    }
}

我很惊讶地发现它比标准的sqrt函数更快(执行时间仅为标准函数的85%左右)。

我不太明白原因,希望能更好地理解。以下是我用于分析性能的完整代码(在Visual Studio 2015中编译Release模式并打开所有优化):

#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 ??                   ?? ?? 

4
可能的原因之一是标准的 sqrt 函数来自于一个链接库,因此无法内联,并且容易受到函数调用开销的影响,而您的实现则可以完全内联化。 - Sergey
1
我在我的电脑上得到了相同的时间。使用/fp:fast选项使库版本快了近两倍。 - Ross Ridge
1
我猜如果我真的想知道答案,我会输出汇编代码(/FAs)并检查代码。 - David Wohlferd
1
/fp:fast 通常不会有太大的危险,除非你特别选择了某些操作的顺序,或者你期望遇到 NaNs 时仍然能得到有用的结果。例如,如果你知道你的数组在大数和小数之间交替出现,并且你不希望自动向量化最终将只把小数加在一起,直到它们足够大以至于大数的总和超过了1 ulp。或者,如果你需要在不同编译器版本的不同构建中获得完全可重复的输出。 - Peter Cordes
1
/fp:fast选项并不会对精度产生太大影响,而是更多地关注符合性。实际上,它可以使结果更加精确。由于你的sqrt实现不符合规范(例如,在对负数进行平方根运算时没有设置errno),所以在这种情况下进行比较是公平的。 - Ross Ridge
显示剩余16条评论
2个回答

10
你的循环都非常糟糕,除了sqrt函数调用或FSQRT指令之外,还存在许多瓶颈。至少比最优标量SQRTSS(单精度)代码慢2倍,并且可能比体面的SSE2向量化循环慢8倍。即使不重新排序任何数学运算,也可以击败SQRTSS吞吐量。 https://gcc.gnu.org/wiki/DontUseInlineAsm中的许多原因适用于您的示例。编译器将无法通过函数传播常量,并且如果结果不是NaN,则不知道结果始终为非负数。如果您稍后平方该数字,它也无法将其优化为 fabs()
同时非常重要的是,您可以使用SSE2 SQRTPS(_mm_sqrt_ps())击败自动向量化。使用内部函数的“无错误检查”标量sqrt()函数也会遇到这个问题。我不知道是否有任何方法可以在没有/fp:fast的情况下获得最佳结果,但我对此表示怀疑。(除了使用汇编编写整个循环,或者使用内部函数将整个循环向量化)。
很不错,你的Haswell CPU能够以如此快的速度运行函数调用循环,尽管内联汇编循环可能甚至未饱和FSQRT吞吐量。由于某种原因,您的库函数调用调用了double sqrt(double),而不是C++重载float sqrt(float)。这导致转换为double然后再转回float。可能需要#include <cmath>以获取重载, 或者您可以调用sqrtf()。在Linux上,gcc和clang使用当前代码调用sqrtf()(而不进行double和float的转换),但是可能它们的<random>头文件恰好包括<cmath>,而MSVC则不包括。或者也许还有其他事情发生了。
库函数调用循环将总和保存在内存中(而不是寄存器中)。显然,32位版本的__libm_sse2_sqrt_precise使用的调用约定不保留任何XMM寄存器。Windows x64 ABI确实保留XMM6-XMM15,但维基百科称这是新的,32位ABI没有这样做。我认为如果有任何调用保留的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个周期一次,因此循环传递的依赖关系不是Haswell的限制瓶颈。


The inline-asm version有一个存储/重新加载sqrt结果的往返,但它不是循环依赖链的一部分。MSVC内联汇编语法使避免存储前转发往返困难以将数据输入/输出内联汇编。更糟的是,你在x87堆栈上产生了结果,而编译器想要在XMM寄存器中进行SSE数学运算。
然后MSVC无缘无故地自食恶果,将总和保留在内存中而不是在XMM寄存器中。它会查看内联汇编语句影响哪些寄存器,所以我不知道为什么它看不到你的内联汇编语句不会破坏任何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 + store-forwarding)。Haswell FSQRT每8-17个周期具有一个吞吐量,因此可能仍然是瓶颈。(涉及数组值的所有存储/重新加载都是每次迭代独立的,并且乱序执行可以重叠许多迭代以隐藏该延迟链。但是,它们将阻塞负载/存储执行单元,并有时通过额外的周期延迟关键路径的负载/存储。这称为资源冲突。)
没有使用/fp:fastsqrtf()库函数必须在结果为NaN时设置errno。这就是为什么它不能内联到一个SQRTSS中。
如果你确实想要自己实现一个无检查的标量平方根函数,你可以使用英特尔指令语法来实现:
// 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)。在Godbolt编译器探索器上进行查看:

# 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 不是真正的可交换的):聪明的编译器(或使用内部函数的程序员)将使用 SQRTPS 来获得 4 个结果,其吞吐量与 SQRTSS 的 1 个结果相同。将 SQRT 结果的向量解包成 4 个标量,然后您可以保持完全相同的操作顺序和中间结果的相同舍入。我对 clang 和 gcc 没有这样做感到失望。

然而,gcc 和 clang 确实成功避免了调用库函数。clang3.9(只使用 -O3)甚至在不检查 NaN 的情况下就使用了 SQRTSS。我认为这是合法的,不是编译器错误。也许它看到代码没有使用 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在内联汇编中一样:存在一个存储转发的循环依赖。所有的溢出/加载操作都可以在JBE跳过的块内部执行。也许gcc认为负输入很常见。


即便您使用了/fp:fast或者-ffast-math,像clang这样的聪明编译器也无法将_mm_sqrt_ss重写为SQRTPS。通常情况下,clang并不只是简单地将内部函数映射到指令,如果您错过了合并操作的机会,它会提供更优化的洗牌和混合功能。

因此,启用快速FP数学后,使用_mm_sqrt_ss将会损失很大。clang将sqrt()库函数调用版本编译为RSQRTPS + 牛顿-拉弗森迭代。


请注意,您的微基准测试代码对sqrt_new()实现的延迟不敏感,只对吞吐量敏感。在真正的FP代码中,延迟通常很重要,而不仅仅是吞吐量。但在其他情况下,比如独立地对许多数组元素执行相同的操作时,延迟并不重要,因为乱序执行可以通过从许多循环迭代中有指令在飞行中来隐藏它。
正如我之前提到的,来自额外存储/重新加载数据在进出MSVC样式内联汇编时所需的往返时间的延迟在这里是一个严重的问题。当MSVC内联函数时,fld n并不直接来自数组。
顺便提一下,Skylake的SQRTPS/SS吞吐量为每3个周期1次,但延迟仍为12个周期。SQRTPD/SD吞吐量为每4-6个周期1次,延迟为15-16个周期。因此,在Skylake上,FP平方根的流水线更多,比Haswell更多。这扩大了基准测试FP平方根延迟与吞吐量之间的差异。

3

以发布模式编译并开启所有优化选项

你没有开启全部优化选项,错过了其中一个。在IDE中,可以通过“项目”、“属性”、“C/C++”、“代码生成”、“浮点数模型”来进行设置。你将其保留在默认设置“/fp:precise”上。这会对生成的机器代码产生非常明显的副作用:

  00091197 E8 F2 18 00 00       call        __libm_sse2_sqrt_precise (092A8Eh) 

也许调用CRT中的帮助函数总是比内联指令如FSQRT慢,这一点很直观。
关于/fp的确切语义有很多要说的,MSDN文章并不是非常好。而且它也很难进行逆向工程,因为Microsoft从Intel购买了代码,并且无法获得允许他们重新发布汇编代码的源代码许可证。它最初的目标肯定是处理由Intel的8087FPU设计引起的可怕的浮点一致性问题。但今天这已经不那么相关了,现在所有主流的C和C++编译器都会生成SSE2代码。自VS2012以来,MSVC++就是这样做的。这些英特尔库函数现在主要确保浮点运算仍然产生与旧版本编译器一致的结果。
__libm_sse2_sqrt_precise()实际上做了很多事情。冒着试图记录未记录的函数的风险,我认为这就是它的功能:
- 检查MXCSR寄存器的值,这是SSE的控制寄存器。如果它没有默认值(0x1F80),那么代码假定程序员使用了_control_fp(),并跳转到计算具有FPU语义的传统sqrt()实现。 - 检查FPU控制寄存器的值,以查看浮点异常是否已启用。通常是关闭的,如果有任何异常被启用,那么它会像上面那样跳转到sqrt()。 - 检查参数是否小于零但不是负零,并调用用户提供的_matherr()函数。 - 最后使用SQRTSD指令计算结果。
这实际上与精度无关:) 看到这个85%性能的执行结果相当不错,尽管FSQRT比SQRTSD慢得多。现代处理器对后者进行了更多的硅爱。
如果你关心快速浮点运算,那么将设置更改为/fp:fast即可。它会产生以下结果:
  00D91310  sqrtsd      xmm0,xmm0 

一个内联指令而不是库调用。换句话说,跳过前面列表中的前3个项目。同时也比FSQRT更加高效。

非常有见地,谢谢!为了测试目的,如果我不能使用/fp:fast来进行整个应用程序,那么指定一个自定义的sqrtsd xmm0,xmm0函数的最佳方法是什么?我认为这将是我情况下的理想比较。我尝试过float sqrt_new(float n) { _mm_cvtss_f32(_mm_sqrt_ss(_mm_set_ss(n))); },但实际上比带有/fp:fast的标准sqrt(n)要慢得多。 - Louis15
该设置适用于每个源代码文件。因此,只需移动“需要快速执行”的代码即可。 - Hans Passant
天啊,那个库函数听起来比我想象的还要糟糕。我认为glibc的libm sqrt只需要检查NaN并设置errno,而不是检查MXCSR。你们有什么想法,为什么OP的代码会转换为“double”呢?这只是因为他没有“#include <cmath>”,所以“sqrt(float)”的浮点重载不可用吗?但是必须有某些东西声明“sqrt(double)”,否则它就无法编译,对吧? - Peter Cordes
@Louis15:/fp:fast 让整个过程比标量平方根更快是正常的。它允许整个过程(包括加法)向量化(通过将FP数学视为可结合的)。你的 sqrt_new 仍然只能一次进行标量平方根,而不能像 SQRTPS 一样一次进行4次运算。如果你的代码仍在调用 sqrt(double) 而不是 sqrt(float),那么可能是 SQRTPD - Peter Cordes
这是一个带有1-800支持电话号码的库函数,很难进行比较。我只关注于原始问题中的机器代码,完全没有尝试源代码。 - Hans Passant

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接