2

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

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 上。

4

2 回答 2

5

imul r32,r32在典型的现代 x86 CPU ( http://agner.org/optimize/ ) 上有 3 个周期延迟。因此，标量实现可以每 3 个时钟周期进行一次乘法，因为它们是相互依赖的。但是，它是完全流水线的，因此您的标量循环使 2/3 的潜在吞吐量未被使用。

在 3 个周期内，Core2 或更高版本中的流水线可以将 12 uop 馈送到核心的无序部分。对于较小的输入，最好保持代码较小，并让无序执行与后续代码重叠依赖链，尤其是在后续代码不完全依赖于阶乘结果的情况下。但是编译器不善于知道何时优化延迟与吞吐量，并且如果没有配置文件引导的优化，他们没有关于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 次打包乘法（英特尔 CPU 上的延迟为 5c / 1c 或 0.5c 吞吐量） 在 AMD CPU 上类似；请参阅 Agner Fog 的说明表。 每个向量循环迭代在您的 C 源代码中执行 4 次阶乘循环迭代，并且在一次迭代中存在显着的指令级并行性

内部循环只有 12 uops 长（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

所以 gcc 浪费了很多精力来模拟一个压缩的 32 位乘法，而不是在使用pmuludq. 我还看了clang6.0。我认为它落入了同一个陷阱。（Godbolt 编译器资源管理器上的 Source+asm）

您没有使用-march=native或任何东西，因此只有 SSE2（x86-64 的基线）可用，因此只有扩展 32x32 => 64 位 SIMD 乘法pmuludq可用于 32 位输入元素。SSE4.1pmulld在 Haswell 和更高版本上是 2 微指令（在 Sandybridge 上是单微指令），但会避免 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 个依赖微指令）。使用多个累加器展开可以为大输入提供高达 10 倍的吞吐量提升，对于 Skylake 上的 SIMD 整数乘法微指令，延迟为 5c/0.5c。

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

默认情况下，Clang 使用多个累加器展开，所以它应该很好。但是代码很多，所以我没有手动分析。将其插入 IACA 或对大输入进行基准测试。（什么是 IACA 以及如何使用它？）

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

阶乘查找表[0..7]可能是最好的选择。安排事情，让你的向量/展开循环做n%8 .. n，而不是1 .. n/8*8，所以剩下的部分总是相同的n。

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

8!是 40320，适合 16 位，因此 1..8 的查找表只需要 8 * 2 字节的存储空间。或者使用 32 位条目，这样您就可以使用内存源操作数imul代替单独的movzx.

于 2018-07-11T07:31:57.973 回答
3

它不会使情况变得更糟。对于大数字，它运行得更快。以下是结果factorial(1000000000)：

  • -O2: 0.78 秒
  • -O3: 0.5 秒

当然，使用这么大的数字是未定义的行为（因为有符号算术溢出）。但是时间与无符号数相同，对于它来说，它不是未定义的行为。

请注意，阶乘的这种用法通常没有意义，因为它不计算num!, 而是num! & UINT_MAX. 但是编译器不知道这一点。

如果使用 PGO，编译器不会向量化这段代码，如果它总是用小数字调用的话。

如果您不喜欢这种行为，但您想使用-O3，请使用 关闭自动矢量化-fno-tree-loop-vectorize。

于 2018-07-11T06:42:23.957 回答