gcc - 为什么gcc -O3自动矢量分解？许多额外的指令看起来更糟

Question

这是一个非常简单的阶乘函数。

int factorial(int num) {
    if (num == 0)
        return 1;
    return num*factorial(num-1);
}

GCC在-O2上为此功能进行的组装是合理的。

factorial(int):
        mov     eax, 1
        test    edi, edi
        je      .L1
.L2:
        imul    eax, edi
        sub     edi, 1
        jne     .L2
.L1:
        ret

但是，在-O3或-Ofast上，它决定使事情变得更复杂（几乎100行！）：

factorial(int):
        test    edi, edi
        je      .L28
        lea     edx, [rdi-1]
        mov     ecx, edi
        cmp     edx, 6
        jbe     .L8
        mov     DWORD PTR [rsp-12], edi
        movd    xmm5, DWORD PTR [rsp-12]
        mov     edx, edi
        xor     eax, eax
        movdqa  xmm0, XMMWORD PTR .LC0[rip]
        movdqa  xmm4, XMMWORD PTR .LC2[rip]
        shr     edx, 2
        pshufd  xmm2, xmm5, 0
        paddd   xmm2, XMMWORD PTR .LC1[rip]
.L5:
        movdqa  xmm3, xmm2
        movdqa  xmm1, xmm2
        paddd   xmm2, xmm4
        add     eax, 1
        pmuludq xmm3, xmm0
        psrlq   xmm1, 32
        psrlq   xmm0, 32
        pmuludq xmm1, xmm0
        pshufd  xmm0, xmm3, 8
        pshufd  xmm1, xmm1, 8
        punpckldq       xmm0, xmm1
        cmp     eax, edx
        jne     .L5
        movdqa  xmm2, xmm0
        movdqa  xmm1, xmm0
        mov     edx, edi
        psrldq  xmm2, 8
        psrlq   xmm0, 32
        and     edx, -4
        pmuludq xmm1, xmm2
        psrlq   xmm2, 32
        sub     edi, edx
        pmuludq xmm0, xmm2
        pshufd  xmm1, xmm1, 8
        pshufd  xmm0, xmm0, 8
        punpckldq       xmm1, xmm0
        movdqa  xmm0, xmm1
        psrldq  xmm1, 4
        pmuludq xmm0, xmm1
        movd    eax, xmm0
        cmp     ecx, edx
        je      .L1
        lea     edx, [rdi-1]
.L3:
        imul    eax, edi
        test    edx, edx
        je      .L1
        imul    eax, edx
        mov     edx, edi
        sub     edx, 2
        je      .L1
        imul    eax, edx
        mov     edx, edi
        sub     edx, 3
        je      .L1
        imul    eax, edx
        mov     edx, edi
        sub     edx, 4
        je      .L1
        imul    eax, edx
        mov     edx, edi
        sub     edx, 5
        je      .L1
        imul    eax, edx
        sub     edi, 6
        je      .L1
        imul    eax, edi
.L1:
        ret
.L28:
        mov     eax, 1
        ret
.L8:
        mov     eax, 1
        jmp     .L3
.LC0:
        .long   1
        .long   1
        .long   1
        .long   1
.LC1:
        .long   0
        .long   -1
        .long   -2
        .long   -3
.LC2:
        .long   -4
        .long   -4
        .long   -4
        .long   -4

我使用Compiler Explorer获得了这些结果，因此在实际的用例中应该是相同的。

那是怎么回事？在任何情况下这会更快吗？ Clang似乎也做了类似的事情，但是在-O2上。

Answer 1

在典型的现代x86 CPU（http://agner.org/optimize/）上，imul r32,r32具有3个周期的延迟。因此标量实现可以每3个时钟周期执行一次乘法运算，因为它们是相关的。但是，它已完全流水线化，因此标量循环未使用潜在吞吐量的2 / 3rd。

在3个周期中，Core2或更高版本中的管道可以将12 uops馈入内核的乱序部分。对于较小的输入，最好使代码保持较小，并让乱序执行将依赖链与更高版本的代码重叠，尤其是如果更高版本的代码并不完全取决于阶乘结果的话。但是，编译器并不善于知道何时针对延迟与吞吐量进行优化，并且如果没有配置文件引导的优化，他们就没有通常n的大小的数据。

我怀疑gcc的自动向量化器没有考虑大型n溢出的速度。



一个有用的标量优化应该是使用多个累加器展开的，例如利用乘法是关联的这一事实，并在循环中并行执行这些操作：prod(n*3/4 .. n) * prod(n/2 .. n*3/4) * prod(n/4 .. n/2) * prod(1..n/4)（当然，具有不重叠的范围）。即使自动换行，乘法也具有关联性。乘积位仅取决于该位置的低位，而不取决于（丢弃的）高位。

或更简单地，执行f0 *= i; f1 *= i+1; f2 *= i+2; f3 *= i+3; i+=4;。然后在循环之外，return (f0*f1) * (f2*f3);。这也将是标量代码的胜利。当然，展开时还必须考虑n % 4 != 0。



gcc选择要做的基本上是后者，使用pmuludq用一条指令进行2次打包乘法（在Intel CPU上为5c延迟/ 1c或0.5c吞吐量）。请参阅Agner Fog的说明表。每个向量循环迭代都会在您的C源代码中执行阶乘循环的4次迭代，并且一次迭代内有明显的指令级并行性

内部循环只有12微秒长（cmp / jcc宏融合为1），因此它可以每3个周期发出1次迭代，与标量版本中的延迟瓶颈具有相同的吞吐量，但每次迭代完成4倍的工作量。

.L5:
    movdqa  xmm3, xmm2         ; copy the old i vector
    movdqa  xmm1, xmm2
    paddd   xmm2, xmm4         ; [ i0,  i1 |  i2,  i3 ]  += 4
    add     eax, 1
    pmuludq xmm3, xmm0         ; [ f0      |  f2  ] *= [ i0   |  i2  ]

    psrlq   xmm1, 32           ; bring odd 32 bit elements down to even: [ i1  | i3 ]
    psrlq   xmm0, 32
    pmuludq xmm1, xmm0         ; [ f1  | f3 ] *= [ i1  | i3 ]

    pshufd  xmm0, xmm3, 8
    pshufd  xmm1, xmm1, 8
    punpckldq       xmm0, xmm1   ; merge back into [ f0  f1  f2  f3 ]
    cmp     eax, edx
    jne     .L5

因此，在使用 pmuludq时，gcc浪费了大量精力来模拟一个打包的32位乘法，而不是将两个单独的向量累加器分开。我也看了clang6.0。我认为它陷入了同样的陷阱。 （ Source+asm on the Godbolt compiler explorer）

您没有使用 -march=native或任何东西，因此仅SSE2（x86-64的基线）可用，因此只有扩展的32x32 => 64位SIMD乘法（如 pmuludq）可用于32位输入元素。 SSE4.1 pmulld在Haswell及之后版本（Sandybridge）上为2 oups，但可以避免gcc的所有愚蠢改组。

当然，这里也存在一个延迟瓶颈，尤其是由于gcc错过了优化，从而增加了累加器所带来的循环承载的dep链的长度。

展开更多的向量累加器可能会隐藏很多 pmuludq延迟。

通过良好的矢量化，SIMD整数乘法器可以管理标量整数乘法单元2倍或4倍的吞吐量。 （或者，对于AVX2，使用8x 32位整数的向量将吞吐量提高8倍。）

但是向量越宽，展开越多，所需的清理代码就越多。



gcc -march=haswell

我们得到一个这样的内循环：

.L5:
    inc     eax
    vpmulld ymm1, ymm1, ymm0
    vpaddd  ymm0, ymm0, ymm2
    cmp     eax, edx
    jne     .L5

超级简单，但带有10c延迟循环承载的依赖链：/（ pmulld是Haswell及更高版本上的2个依赖uop）。对于大型输入，使用多个累加器进行展开可以使吞吐量提高10倍，对于Skylake上的SIMD整数乘法，则可以达到5c延迟/ 0.5c吞吐量。

但是对于标量，每5个周期4个乘数仍然比每3个周期1个乘数好得多。

默认情况下，Clang会使用多个累加器展开，因此应该不错。但这是很多代码，所以我没有手工分析它。将其插入IACA或对大型输入进行基准测试。 （ What is IACA and how do I use it?）



处理展开尾声的有效策略：

查找阶乘 [0..7]的表可能是最好的选择。安排事情，使向量/展开循环执行 n%8 .. n而不是 1 .. n/8*8，因此每个 n的剩余部分始终相同。

在水平向量乘积之后，再对表查找结果进行一个标量乘法。 SIMD循环已经需要一些向量常量，因此无论如何您都可能会触摸内存，并且表查找可以与主循环并行进行。

8!是40320，适合16位，因此1..8查找表仅需要8 * 2字节的存储空间。或使用32位条目，以便可以将内存源操作数用于 imul而不是单独的 movzx。