为什么 gcc -O3 自动矢量化阶乘？许多额外的指令看起来更糟
Why is gcc -O3 auto-vectorizing factorial? That many extra instructions looks worse
这是一个非常简单的阶乘函数。
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 上。
这并没有使情况变得更糟。对于大量数据，它运行得更快。以下是 factorial(1000000000) 的结果：
-O2: 0.78 秒-O3: 0.5 秒
当然，使用这么大的数字是未定义的行为（因为有符号算术溢出）。但是时间与无符号数相同，这不是未定义的行为。
请注意，阶乘的这种用法通常没有意义，因为它不计算 num!，而是计算 num! & UINT_MAX。但是编译器并不知道这一点。
也许对于 PGO，编译器不会向量化这段代码，如果它总是用小数字调用的话。
如果您不喜欢这种行为，但又想使用 -O3，请使用 -fno-tree-loop-vectorize 关闭自动矢量化。
imul r32,r32 在典型的现代 x86 CPU (http://agner.org/optimize/) 上有 3 个周期延迟。所以标量实现可以每 3 个时钟周期执行一次乘法，因为它们是相关的。不过，它是完全流水线化的，因此您的标量循环留下了 2/3 的潜在吞吐量未使用。
在 3 个周期内，Core2 或更高版本中的流水线可以将 12 微指令馈送到内核的乱序部分。对于小输入，最好保持代码小，让无序执行与后面的代码重叠依赖链，特别是如果后面的代码并不完全依赖于阶乘结果。但是编译器并不擅长知道何时针对延迟与吞吐量进行优化，并且如果没有配置文件引导的优化，他们就没有关于 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 次打包乘法（5c 延迟/1c 或 0.5c Intel CPU 上的吞吐量） 在 AMD CPU 上类似；请参阅 Agner Fog 的说明 tables。 每个向量循环迭代执行 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 和更高版本上是 2 微指令（在 Sandybridge 上是单微指令），但会避免所有 gcc 的愚蠢改组。
当然这里也存在延迟瓶颈，尤其是因为 gcc 错过了优化，增加了涉及累加器的循环携带的 dep 链的长度。
使用更多向量累加器展开可以隐藏很多 pmuludq 延迟。
通过良好的矢量化，SIMD 整数乘法器可以管理 2 倍或 4 倍的标量整数乘法单元的吞吐量。 （或者，对于 AVX2，使用 8 个 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 个依赖 uops）。使用多个累加器展开可以为大型输入提供高达 10 倍的吞吐量提升，对于 Skylake 上的 SIMD 整数乘法运算，5c 延迟/0.5c 吞吐量。
但是对于标量，每 5 个周期 4 次乘法仍然比每 3 次 1 次好得多。
Clang 默认使用多个累加器展开，所以应该不错。不过代码比较多，就没有手动分析了。将其插入 IACA 或针对大输入进行基准测试。 (What is IACA and how do I use it?)
处理展开结语的有效策略：
查找 table 阶乘 [0..7] 可能是最好的选择。安排事物，使您的矢量/展开循环执行 n%8 .. n，而不是 1 .. n/8*8，因此每个 n.
的剩余部分始终相同
水平矢量乘积后，再与table查找结果进行一次标量乘法运算。 SIMD 循环已经需要一些向量常量，因此您可能无论如何都会接触到内存，并且 table 查找可以与主循环并行进行。
8! 是 40320，适合 16 位，所以 1..8 查找 table 只需要 8 * 2 字节的存储空间。或者使用 32 位条目，这样您就可以使用 imul 的内存源操作数，而不是单独的 movzx.
这是一个非常简单的阶乘函数。
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 上。
这并没有使情况变得更糟。对于大量数据，它运行得更快。以下是 factorial(1000000000) 的结果：
-O2: 0.78 秒-O3: 0.5 秒
当然，使用这么大的数字是未定义的行为（因为有符号算术溢出）。但是时间与无符号数相同，这不是未定义的行为。
请注意，阶乘的这种用法通常没有意义，因为它不计算 num!，而是计算 num! & UINT_MAX。但是编译器并不知道这一点。
也许对于 PGO，编译器不会向量化这段代码，如果它总是用小数字调用的话。
如果您不喜欢这种行为，但又想使用 -O3，请使用 -fno-tree-loop-vectorize 关闭自动矢量化。
imul r32,r32 在典型的现代 x86 CPU (http://agner.org/optimize/) 上有 3 个周期延迟。所以标量实现可以每 3 个时钟周期执行一次乘法，因为它们是相关的。不过，它是完全流水线化的，因此您的标量循环留下了 2/3 的潜在吞吐量未使用。
在 3 个周期内，Core2 或更高版本中的流水线可以将 12 微指令馈送到内核的乱序部分。对于小输入，最好保持代码小，让无序执行与后面的代码重叠依赖链，特别是如果后面的代码并不完全依赖于阶乘结果。但是编译器并不擅长知道何时针对延迟与吞吐量进行优化，并且如果没有配置文件引导的优化，他们就没有关于 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 次打包乘法（5c 延迟/1c 或 0.5c Intel CPU 上的吞吐量） 在 AMD CPU 上类似；请参阅 Agner Fog 的说明 tables。 每个向量循环迭代执行 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 和更高版本上是 2 微指令（在 Sandybridge 上是单微指令），但会避免所有 gcc 的愚蠢改组。
当然这里也存在延迟瓶颈，尤其是因为 gcc 错过了优化，增加了涉及累加器的循环携带的 dep 链的长度。
使用更多向量累加器展开可以隐藏很多 pmuludq 延迟。
通过良好的矢量化，SIMD 整数乘法器可以管理 2 倍或 4 倍的标量整数乘法单元的吞吐量。 （或者，对于 AVX2，使用 8 个 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 个依赖 uops）。使用多个累加器展开可以为大型输入提供高达 10 倍的吞吐量提升，对于 Skylake 上的 SIMD 整数乘法运算，5c 延迟/0.5c 吞吐量。
但是对于标量，每 5 个周期 4 次乘法仍然比每 3 次 1 次好得多。
Clang 默认使用多个累加器展开，所以应该不错。不过代码比较多，就没有手动分析了。将其插入 IACA 或针对大输入进行基准测试。 (What is IACA and how do I use it?)
处理展开结语的有效策略：
查找 table 阶乘 [0..7] 可能是最好的选择。安排事物，使您的矢量/展开循环执行 n%8 .. n，而不是 1 .. n/8*8，因此每个 n.
水平矢量乘积后，再与table查找结果进行一次标量乘法运算。 SIMD 循环已经需要一些向量常量，因此您可能无论如何都会接触到内存，并且 table 查找可以与主循环并行进行。
8! 是 40320，适合 16 位，所以 1..8 查找 table 只需要 8 * 2 字节的存储空间。或者使用 32 位条目，这样您就可以使用 imul 的内存源操作数，而不是单独的 movzx.