我用汇编和 C++ 为 Project Euler Q14 编写了这两个解决方案。他们实现了相同的蛮力方法来测试 Collatz conjecture。组装解决方案由以下组件组装而成:
nasm -felf64 p14.asm && gcc p14.o -o p14
C++ 编译时使用:
g++ p14.cpp -o p14
大会,p14.asm
:
section .data
fmt db "%d", 10, 0
global main
extern printf
section .text
main:
mov rcx, 1000000
xor rdi, rdi ; max i
xor rsi, rsi ; i
l1:
dec rcx
xor r10, r10 ; count
mov rax, rcx
l2:
test rax, 1
jpe even
mov rbx, 3
mul rbx
inc rax
jmp c1
even:
mov rbx, 2
xor rdx, rdx
div rbx
c1:
inc r10
cmp rax, 1
jne l2
cmp rdi, r10
cmovl rdi, r10
cmovl rsi, rcx
cmp rcx, 2
jne l1
mov rdi, fmt
xor rax, rax
call printf
ret
C++,p14.cpp
:
#include <iostream>
int sequence(long n) {
int count = 1;
while (n != 1) {
if (n % 2 == 0)
n /= 2;
else
n = 3*n + 1;
++count;
}
return count;
}
int main() {
int max = 0, maxi;
for (int i = 999999; i > 0; --i) {
int s = sequence(i);
if (s > max) {
max = s;
maxi = i;
}
}
std::cout << maxi << std::endl;
}
我知道编译器优化以提高速度和一切,但我没有看到很多方法来进一步优化我的汇编解决方案(以编程方式,而不是数学方式)。
C++ 代码每项使用模数,每隔一项使用除法,而汇编代码每隔一项仅使用一次除法。
但该程序集的平均时间比 C++ 解决方案长 1 秒。为什么是这样?我主要是出于好奇而问。
执行时间
我的系统:1.4 GHz Intel Celeron 2955U(Haswell 微架构)上的 64 位 Linux。
g++(未优化):平均 1272 毫秒。
g++ -O3:平均 578 毫秒。
asm(div)(原始):平均 2650 毫秒。
asm (shr):平均 679 毫秒。
@johnfound asm(使用 NASM 组装):平均 501 毫秒。
@hidefromkgb asm:平均 200 毫秒。
@hidefromkgb asm,由 @Peter Cordes 优化:平均 145 毫秒。
@Veedrac C++:-O3 平均 81 毫秒,-O0 平均 305 毫秒。
-S
进行编译以获取编译器生成的程序集。编译器足够聪明,可以意识到模数同时进行除法。
如果您认为 64 位 DIV 指令是除以二的好方法,那么难怪编译器的 asm 输出会击败您的手写代码,即使使用 -O0
(编译速度快,无需额外优化,并且存储/重新加载到每个 C 语句之后/之前的内存,以便调试器可以修改变量)。
请参阅 Agner Fog's Optimizing Assembly guide 了解如何编写高效的 asm。他还提供指令表和微架构指南,了解特定 CPU 的具体细节。另请参阅 x86 标记 wiki 以获取更多性能链接。
另请参阅有关使用手写 asm 击败编译器的更一般的问题:Is inline assembly language slower than native C++ code?。 TL:DR:是的,如果你做错了(比如这个问题)。
通常你可以让编译器完成它的工作,特别是如果你尝试编写可以高效编译的 C++。另见is assembly faster than compiled languages?。其中一个答案链接到 these neat slides,展示了各种 C 编译器如何使用很酷的技巧优化一些非常简单的函数。 Matt Godbolt 的 CppCon2017 演讲“What Has My Compiler Done for Me Lately? Unbolting the Compiler's Lid”与此类似。
even:
mov rbx, 2
xor rdx, rdx
div rbx
在 Intel Haswell 上,div r64
为 36 微秒,延迟为 32-96 个周期,吞吐量为每 21-74 个周期一个。 (加上设置 RBX 和零 RDX 的 2 个微指令,但乱序执行可以提前运行这些)。 High-uop-count instructions like DIV are microcoded, which can also cause front-end bottlenecks. 在这种情况下,延迟是最相关的因素,因为它是循环承载的依赖链的一部分。
shr rax, 1
执行相同的无符号除法:它是 1 uop,具有 1c 延迟,每个时钟周期可以运行 2 次。
相比之下,32 位除法速度更快,但与移位相比仍然很糟糕。 idiv r32
在 Haswell 上是 9 微指令、22-29c 延迟和每 8-11c 吞吐量一个。
从 gcc 的 -O0
asm 输出 (Godbolt compiler explorer) 可以看出,它只使用移位指令。 clang -O0
确实像您想象的那样天真地编译,即使使用 64 位 IDIV 两次。 (在优化时,如果源使用相同的操作数进行除法和取模,编译器确实会使用 IDIV 的两个输出,如果它们完全使用 IDIV)
GCC 没有完全幼稚的模式。 it always transforms through GIMPLE, which means some "optimizations" can't be disabled。这包括识别除以常数和使用移位(2 的幂)或 a fixed-point multiplicative inverse(非 2 的幂)来避免 IDIV(参见上面的 gobolt 链接中的 div_by_13
)。
gcc -Os
(针对大小进行优化)确实使用 IDIV 进行非 2 次幂除法,不幸的是,即使在乘法逆码稍大但速度更快的情况下也是如此。
帮助编译器
(本案例总结:使用 uint64_t n
)
首先,查看优化的编译器输出很有趣。 (-O3
)。
-O0
speed is basically meaningless.
查看您的 asm 输出(在 Godbolt 上,或参见 How to remove "noise" from GCC/clang assembly output?)。当编译器一开始就没有生成最佳代码时:以引导编译器生成更好代码的方式编写 C/C++ 源代码通常是最好的方法。您必须了解 asm,并且知道什么是有效的,但是您间接地应用了这些知识。编译器也是一个很好的想法来源:有时 clang 会做一些很酷的事情,你可以手持 gcc 做同样的事情:参见 this answer 以及我在下面@Veedrac 的代码中对非展开循环所做的事情。)
这种方法是可移植的,并且在 20 年内,一些未来的编译器可以将其编译为在未来硬件(x86 或非 x86 上)上有效的任何东西,可能使用新的 ISA 扩展或自动矢量化。 15 年前手写的 x86-64 asm 通常不会针对 Skylake 进行优化调整。例如,当时不存在比较和分支宏融合。 目前对于一个微架构的手工 asm 来说最佳的做法可能不适用于其他当前和未来的 CPU。 Comments on @johnfound's answer 讨论 AMD Bulldozer 和 Intel Haswell 之间的主要区别,这对这段代码有很大影响.但理论上,g++ -O3 -march=bdver3
和 g++ -O3 -march=skylake
会做正确的事。 (或 -march=native
。)或 -mtune=...
只是调整,而不使用其他 CPU 可能不支持的指令。
我的感觉是,引导编译器编译出对您关心的当前 CPU 有好处的汇编,对于未来的编译器来说应该不是问题。他们希望在寻找转换代码的方法方面比当前的编译器更好,并且可以找到一种适用于未来 CPU 的方法。无论如何,未来的 x86 可能不会在当前 x86 上的任何优点方面变得糟糕,并且未来的编译器将避免任何特定于 asm 的陷阱,同时实现诸如从 C 源代码移动数据之类的东西,如果它没有看到更好的东西。
手写 asm 是优化器的黑盒,因此当内联使输入成为编译时常量时,常量传播不起作用。其他优化也会受到影响。在使用 asm 之前阅读 https://gcc.gnu.org/wiki/DontUseInlineAsm。 (并避免 MSVC 风格的内联汇编:输入/输出必须通过内存 which adds overhead。)
在这种情况下:您的 n
具有带符号的类型,并且 gcc 使用 SAR/SHR/ADD 序列来提供正确的舍入。 (对于负输入,IDIV 和算术移位“舍入”不同,请参阅 SAR insn set ref manual entry)。 (如果 gcc 尝试并未能证明 n
不能为负,则 IDK。签名溢出是未定义的行为,所以它应该能够。)
您应该使用 uint64_t n
,所以它可以只是 SHR。因此它可以移植到 long
仅为 32 位的系统(例如 x86-64 Windows)。
顺便说一句,gcc 的 优化 asm 输出看起来相当不错(使用 unsigned long n
):它内联到 main()
的内部循环这样做:
# from gcc5.4 -O3 plus my comments
# edx= count=1
# rax= uint64_t n
.L9: # do{
lea rcx, [rax+1+rax*2] # rcx = 3*n + 1
mov rdi, rax
shr rdi # rdi = n>>1;
test al, 1 # set flags based on n%2 (aka n&1)
mov rax, rcx
cmove rax, rdi # n= (n%2) ? 3*n+1 : n/2;
add edx, 1 # ++count;
cmp rax, 1
jne .L9 #}while(n!=1)
cmp/branch to update max and maxi, and then do the next n
内循环是无分支的,循环携带的依赖链的关键路径为:
3 组分 LEA(3 个周期)
cmov(Haswell 上 2 个周期,Broadwell 上 1c 或更高版本)。
总计:每次迭代 5 个周期,延迟瓶颈。乱序执行与此并行处理其他所有事情(理论上:我还没有使用性能计数器测试它是否真的以 5c/iter 运行)。
cmov
的 FLAGS 输入(由 TEST 生成)比 RAX 输入(来自 LEA->MOV)的生成速度更快,因此它不在关键路径上。
同样,产生 CMOV 的 RDI 输入的 MOV->SHR 不在关键路径上,因为它也比 LEA 快。 IvyBridge 上的 MOV 和更高版本的延迟为零(在寄存器重命名时处理)。 (它仍然需要一个微指令和一个管道中的插槽,所以它不是免费的,只是零延迟)。 LEA dep 链中的额外 MOV 是其他 CPU 瓶颈的一部分。
cmp/jne 也不是关键路径的一部分:它不是循环携带的,因为控制依赖是通过分支预测 + 推测执行来处理的,这与关键路径上的数据依赖不同。
击败编译器
GCC 在这里做得很好。它可以通过使用 inc edx
instead of add edx, 1
节省一个代码字节,因为没有人关心 P4 及其对部分标志修改指令的错误依赖性。
它还可以保存所有 MOV 指令,并且 TEST: SHR 设置 CF= 移出的位,因此我们可以使用 cmovc
而不是 test
/ cmovz
。
### Hand-optimized version of what gcc does
.L9: #do{
lea rcx, [rax+1+rax*2] # rcx = 3*n + 1
shr rax, 1 # n>>=1; CF = n&1 = n%2
cmovc rax, rcx # n= (n&1) ? 3*n+1 : n/2;
inc edx # ++count;
cmp rax, 1
jne .L9 #}while(n!=1)
请参阅@johnfound 的另一个聪明技巧的答案:通过在 SHR 的标志结果上进行分支以及将其用于 CMOV 来删除 CMP:仅当 n 为 1(或 0)开始时才为零。 (有趣的事实:SHR with count != 1 on Nehalem or earlier causes a stall if you read the flag results。这就是他们制作单微指令的方式。不过,shift-by-1 特殊编码很好。)
避免 MOV 对 Haswell (Can x86's MOV really be "free"? Why can't I reproduce this at all?) 上的延迟没有任何帮助。它确实对 Intel pre-IvB 和 AMD Bulldozer 系列等 CPU 有显着的帮助,其中 MOV 不是零延迟(以及带有更新微码的 Ice Lake)。编译器浪费的 MOV 指令确实会影响关键路径。 BD 的 complex-LEA 和 CMOV 都具有较低的延迟(分别为 2c 和 1c),因此它占延迟的比例更大。此外,吞吐量瓶颈成为一个问题,因为它只有两个整数 ALU 管道。 See @johnfound's answer,他有来自 AMD CPU 的计时结果。
即使在 Haswell 上,此版本也可能会有所帮助,因为它可以避免一些偶然的延迟,即非关键微指令从关键路径上的一个执行端口窃取执行端口,从而将执行延迟 1 个周期。 (这称为资源冲突)。它还保存了一个寄存器,这在交错循环中并行执行多个 n
值时可能会有所帮助(见下文)。
LEA 的延迟取决于寻址模式,在 Intel SnB 系列 CPU 上。 3c 用于 3 个组件([base+idx+const]
,需要两个单独的添加),但只有 1c 具有 2 个或更少的组件(一个添加)。一些 CPU(如 Core2)甚至可以在一个周期内执行 3 分量 LEA,但 SnB 系列却没有。更糟糕的是,Intel SnB-family standardizes latencies so there are no 2c uops,否则 3 分量 LEA 将像 Bulldozer 一样只有 2c。 (AMD 上的 3 组件 LEA 也较慢,只是没那么慢)。
因此,在 Haswell 等 Intel SnB 系列 CPU 上,lea rcx, [rax + rax*2]
/ inc rcx
的延迟只有 2c,比 lea rcx, [rax + rax*2 + 1]
快。 BD 实现收支平衡,Core2 则更糟。它确实需要额外的 uop,这通常不值得节省 1c 延迟,但延迟是这里的主要瓶颈,Haswell 有足够宽的管道来处理额外的 uop 吞吐量。
gcc、icc 和 clang(在 Godbolt 上)都没有使用 SHR 的 CF 输出,总是使用 AND 或 TEST。愚蠢的编译器。 :P 它们是复杂机器的伟大部件,但聪明的人通常可以在小规模问题上击败它们。 (当然,考虑它的时间要长数千到数百万倍!编译器不会使用详尽的算法来搜索所有可能的做事方式,因为在优化大量内联代码时,这将花费太长时间,这就是他们做得最好。他们也没有在目标微架构中对管道进行建模,至少不像 IACA 或其他静态分析工具那样详细;他们只是使用一些启发式方法。)
简单的循环展开无济于事;这个循环瓶颈是循环携带的依赖链的延迟,而不是循环开销/吞吐量。这意味着它可以很好地处理超线程(或任何其他类型的 SMT),因为 CPU 有大量时间来交错来自两个线程的指令。这意味着在 main
中并行化循环,但这很好,因为每个线程都可以检查一系列 n
值并生成一对整数作为结果。
在单个线程中手动交错也可能是可行的。也许并行计算一对数字的序列,因为每个数字只需要几个寄存器,它们都可以更新相同的 max
/ maxi
。这会创建更多 instruction-level parallelism。
诀窍是决定是等到所有 n
值都达到 1
后再获取另一对起始 n
值,还是突破并仅为达到结束条件的一个值获取新起点,无需接触其他序列的寄存器。可能最好让每个链都处理有用的数据,否则你必须有条件地增加它的计数器。
您甚至可以使用 SSE 打包比较内容来执行此操作,以有条件地增加 n
尚未达到 1
的向量元素的计数器。然后为了隐藏 SIMD 条件增量实现的更长延迟,您需要保持更多的 n
值向量悬而未决。可能只值 256b 向量(4x uint64_t
)。
我认为检测 1
“粘性”的最佳策略是掩盖您添加以增加计数器的全一向量。因此,在您看到元素中的 1
之后,增量向量将具有零,并且 +=0 是无操作的。
手动矢量化的未经测试的想法
# starting with YMM0 = [ n_d, n_c, n_b, n_a ] (64-bit elements)
# ymm4 = _mm256_set1_epi64x(1): increment vector
# ymm5 = all-zeros: count vector
.inner_loop:
vpaddq ymm1, ymm0, xmm0
vpaddq ymm1, ymm1, xmm0
vpaddq ymm1, ymm1, set1_epi64(1) # ymm1= 3*n + 1. Maybe could do this more efficiently?
vpsllq ymm3, ymm0, 63 # shift bit 1 to the sign bit
vpsrlq ymm0, ymm0, 1 # n /= 2
# FP blend between integer insns may cost extra bypass latency, but integer blends don't have 1 bit controlling a whole qword.
vpblendvpd ymm0, ymm0, ymm1, ymm3 # variable blend controlled by the sign bit of each 64-bit element. I might have the source operands backwards, I always have to look this up.
# ymm0 = updated n in each element.
vpcmpeqq ymm1, ymm0, set1_epi64(1)
vpandn ymm4, ymm1, ymm4 # zero out elements of ymm4 where the compare was true
vpaddq ymm5, ymm5, ymm4 # count++ in elements where n has never been == 1
vptest ymm4, ymm4
jnz .inner_loop
# Fall through when all the n values have reached 1 at some point, and our increment vector is all-zero
vextracti128 ymm0, ymm5, 1
vpmaxq .... crap this doesn't exist
# Actually just delay doing a horizontal max until the very very end. But you need some way to record max and maxi.
您可以并且应该使用内在函数而不是手写 asm 来实现它。
算法/实现改进:
除了用更高效的 asm 实现相同的逻辑之外,还要寻找简化逻辑或避免冗余工作的方法。例如 memoize 以检测序列的共同结尾。或者更好的是,一次查看 8 个尾随位(gnasher 的回答)
@EOF 指出 tzcnt
(或 bsf
)可用于在一个步骤中执行多个 n/=2
迭代。这可能比 SIMD 矢量化更好;没有 SSE 或 AVX 指令可以做到这一点。不过,它仍然兼容在不同的整数寄存器中并行执行多个标量 n
。
所以循环可能看起来像这样:
goto loop_entry; // C++ structured like the asm, for illustration only
do {
n = n*3 + 1;
loop_entry:
shift = _tzcnt_u64(n);
n >>= shift;
count += shift;
} while(n != 1);
这可能会显着减少迭代次数,但在没有 BMI2 的英特尔 SnB 系列 CPU 上,变量计数的变化很慢。 3 微指令,2c 延迟。 (它们对 FLAGS 有输入依赖性,因为 count=0 表示标志未修改。它们将其作为数据依赖性处理,并采用多个 uop,因为 uop 只能有 2 个输入(无论如何都是预 HSW/BDW))。这就是那些抱怨 x86 疯狂的 CISC 设计的人们所指的那种。它使 x86 CPU 比现在从头开始设计 ISA 的速度要慢,即使是以最相似的方式。 (即这是消耗速度/功率的“x86 税”的一部分。) SHRX/SHLX/SARX (BMI2) 是一个巨大的胜利(1 uop / 1c 延迟)。
它还将 tzcnt(Haswell 上的 3c 及更高版本)置于关键路径上,因此它显着延长了循环承载的依赖链的总延迟。不过,它确实消除了对 CMOV 或准备保存 n>>1
的寄存器的任何需要。 @Veedrac 的回答通过推迟多次迭代的 tzcnt/shift 克服了所有这些问题,这是非常有效的(见下文)。
我们可以安全地交替使用 BSF 或 TZCNT,因为此时 n
永远不会为零。 TZCNT 的机器代码在不支持 BMI1 的 CPU 上解码为 BSF。 (无意义的前缀被忽略,所以 REP BSF 作为 BSF 运行)。
TZCNT 在支持它的 AMD CPU 上的性能比 BSF 好得多,因此使用 REP BSF
可能是个好主意,即使您不关心在输入为零而不是输出时设置 ZF。一些编译器在您使用 __builtin_ctzll
时甚至使用 -mno-bmi
时也会这样做。
它们在 Intel CPU 上的性能相同,所以如果这很重要,只需保存字节即可。 Intel 上的 TZCNT(Skylake 之前)仍然对所谓的只写输出操作数具有错误依赖性,就像 BSF 一样,以支持输入 = 0 的 BSF 未修改其目标的未记录行为。所以你需要解决这个问题,除非只针对 Skylake 进行优化,所以额外的 REP 字节没有任何好处。 (英特尔经常超出 x86 ISA 手册的要求,以避免破坏广泛使用的代码,这些代码依赖于它不应该依赖的东西,或者追溯不允许。例如 Windows 9x's assumes no speculative prefetching of TLB entries,在编写代码时它是安全的, before Intel updated the TLB management rules。)
无论如何,Haswell 上的 LZCNT/TZCNT 与 POPCNT 具有相同的错误深度:参见 this Q&A。这就是为什么在 gcc 的 @Veedrac 代码的 asm 输出中,当它不使用 dst=src 时,您会在即将用作 TZCNT 目标的寄存器上看到它 breaking the dep chain with xor-zeroing。由于 TZCNT/LZCNT/POPCNT 永远不会离开未定义或未修改的目的地,因此这种对 Intel CPU 输出的错误依赖是性能错误/限制。据推测,值得一些晶体管/电源让它们像其他进入同一执行单元的微指令一样工作。唯一的性能优势是与另一个 uarch 限制的交互:they can micro-fuse a memory operand with an indexed addressing mode 在 Haswell 上,但在 Skylake 上,英特尔删除了 LZCNT/TZCNT 的错误 dep,它们“取消层压”索引寻址模式,而 POPCNT 仍然可以微融合任何 addr 模式。
其他答案对想法/代码的改进:
@hidefromkgb 的回答 有一个很好的观察结果,即保证您能够在 3n+1 之后进行一次右移。您可以更有效地计算这一点,而不是仅仅省略步骤之间的检查。但是,该答案中的 asm 实现被破坏了(它取决于 OF,它在 SHRD 之后未定义且计数 > 1),并且速度慢:ROR rdi,2
比 SHRD rdi,rdi,2
快,并且在关键节点上使用两个 CMOV 指令path 比可以并行运行的额外 TEST 慢。
我将经过整理/改进的 C(引导编译器生成更好的 asm),并在 Godbolt 上测试+工作更快的 asm(在 C 下方的注释中):请参阅 @hidefromkgb's answer 中的链接。 (这个答案达到了大型 Godbolt URL 的 30k 字符限制,但是 shortlinks can rot 并且对于 goo.gl 来说太长了。)
还改进了输出打印以转换为字符串并生成一个 write()
而不是一次写入一个字符。这将使用 perf stat ./collatz
(记录性能计数器)对整个程序计时的影响降至最低,并且我对一些非关键 asm 进行了去混淆处理。
@Veedrac 的代码
正如我们所知道的那样,我从右移获得了轻微的加速,并检查以继续循环。在 Core2Duo (Merom) 上,从 limit=1e8 的 7.5 秒降低到 7.275 秒,展开因子为 16。
代码 + 注释 on Godbolt。请勿将此版本与 clang 一起使用;它对延迟循环做了一些愚蠢的事情。使用 tmp 计数器 k
然后将其添加到 count
会改变 clang 的作用,但是 稍微 会伤害 gcc。
请参阅评论中的讨论: Veedrac 的代码在具有 BMI1 的 CPU 上非常出色(即不是 Celeron/Pentium)
声称 C++ 编译器可以生成比有能力的汇编语言程序员更优化的代码是一个非常严重的错误。尤其是在这种情况下。人类总是可以使代码比编译器更好,这种特殊情况很好地说明了这种说法。
您看到的时间差异是因为问题中的汇编代码在内部循环中远非最佳。
(以下代码是 32 位的,但可以很容易地转换为 64 位)
例如,序列函数可以优化为只有 5 条指令:
.seq:
inc esi ; counter
lea edx, [3*eax+1] ; edx = 3*n+1
shr eax, 1 ; eax = n/2
cmovc eax, edx ; if CF eax = edx
jnz .seq ; jmp if n<>1
整个代码如下所示:
include "%lib%/freshlib.inc"
@BinaryType console, compact
options.DebugMode = 1
include "%lib%/freshlib.asm"
start:
InitializeAll
mov ecx, 999999
xor edi, edi ; max
xor ebx, ebx ; max i
.main_loop:
xor esi, esi
mov eax, ecx
.seq:
inc esi ; counter
lea edx, [3*eax+1] ; edx = 3*n+1
shr eax, 1 ; eax = n/2
cmovc eax, edx ; if CF eax = edx
jnz .seq ; jmp if n<>1
cmp edi, esi
cmovb edi, esi
cmovb ebx, ecx
dec ecx
jnz .main_loop
OutputValue "Max sequence: ", edi, 10, -1
OutputValue "Max index: ", ebx, 10, -1
FinalizeAll
stdcall TerminateAll, 0
为了编译此代码,需要 FreshLib。
在我的测试中(1 GHz AMD A4-1200 处理器),上述代码比问题中的 C++ 代码快大约四倍(使用 -O0
编译时:430 ms vs. 1900 ms)等等比使用 -O3
编译 C++ 代码时快两倍(430 ms vs. 830 ms)。
两个程序的输出是相同的:最大序列 = 525 on i = 837799。
-O3
输出时错过了这一点,但我确实发现了您对内部循环所做的所有其他优化。 (但是为什么你使用 LEA 而不是 INC 来进行计数器增量呢?此时可以破坏标志,并导致除 P4 之外的任何东西都减速(对 INC 和 SHR 的旧标志的错误依赖)。LEA 可以t 在尽可能多的端口上运行,并且可能导致资源冲突更频繁地延迟关键路径。)
为了获得更高的性能:一个简单的变化是观察在 n = 3n+1 之后,n 将是偶数,因此您可以立即除以 2。而且 n 不会是 1,所以你不需要测试它。因此,您可以保存一些 if 语句并编写:
while (n % 2 == 0) n /= 2;
if (n > 1) for (;;) {
n = (3*n + 1) / 2;
if (n % 2 == 0) {
do n /= 2; while (n % 2 == 0);
if (n == 1) break;
}
}
这是一个很大的胜利:如果您查看 n 的最低 8 位,则除以 2 八倍之前的所有步骤完全由这 8 位决定。例如,如果最后八位是 0x01,那么在二进制中你的数字是 ???? 0000 0001 然后接下来的步骤是:
3n+1 -> ???? 0000 0100
/ 2 -> ???? ?000 0010
/ 2 -> ???? ??00 0001
3n+1 -> ???? ??00 0100
/ 2 -> ???? ???0 0010
/ 2 -> ???? ???? 0001
3n+1 -> ???? ???? 0100
/ 2 -> ???? ???? ?010
/ 2 -> ???? ???? ??01
3n+1 -> ???? ???? ??00
/ 2 -> ???? ???? ???0
/ 2 -> ???? ???? ????
所以所有这些步骤都可以预测,并且 256k + 1 被 81k + 1 替换。所有组合都会发生类似的事情。所以你可以用一个大的 switch 语句创建一个循环:
k = n / 256;
m = n % 256;
switch (m) {
case 0: n = 1 * k + 0; break;
case 1: n = 81 * k + 1; break;
case 2: n = 81 * k + 1; break;
...
case 155: n = 729 * k + 425; break;
...
}
运行循环直到 n ≤ 128,因为此时 n 可能变为 1,而除以 2 的次数少于 8 步,并且一次执行 8 步或更多步会使您错过第一次达到 1 的点。然后继续“正常”循环 - 或者准备一个表格,告诉您需要多少步才能达到 1。
PS。我强烈怀疑 Peter Cordes 的建议会使其更快。除了一个之外,根本不会有条件分支,并且除非循环实际结束,否则将正确预测一个。所以代码会是这样的
static const unsigned int multipliers [256] = { ... }
static const unsigned int adders [256] = { ... }
while (n > 128) {
size_t lastBits = n % 256;
n = (n >> 8) * multipliers [lastBits] + adders [lastBits];
}
在实践中,您将衡量一次处理 n 的最后 9、10、11、12 位是否会更快。对于每一位,表中的条目数会增加一倍,并且当表不再适合 L1 缓存时,我希望速度会变慢。
聚苯乙烯。如果您需要操作的数量:在每次迭代中,我们准确地进行 8 除以 2,以及可变数量的 (3n + 1) 操作,因此计算操作的明显方法是另一个数组。但我们实际上可以计算步数(基于循环的迭代次数)。
我们可以稍微重新定义这个问题:如果是奇数,则将 n 替换为 (3n + 1) / 2,如果是偶数,则将 n 替换为 n / 2。然后每次迭代将精确执行 8 步,但您可以考虑作弊:-) 所以假设有 r 操作 n <- 3n+1 和 s 操作 n <- n/2。结果将完全是 n' = n * 3^r / 2^s,因为 n <- 3n+1 意味着 n <- 3n * (1 + 1/3n)。取对数,我们发现 r = (s + log2 (n' / n)) / log2 (3)。
如果我们执行循环直到 n ≤ 1,000,000 并且有一个预先计算的表,从任何起点 n ≤ 1,000,000 需要多少次迭代,那么如上所述计算 r,四舍五入到最接近的整数,将给出正确的结果,除非 s 真的很大。
count
,您需要第三个数组,对吗? adders[]
没有告诉您完成了多少次右移。
uint16_t
的零扩展负载非常便宜。在 x86 上,它与从 32 位 unsigned int
到 uint64_t
的零扩展一样便宜。 (英特尔 CPU 上的内存中的 MOVZX 只需要一个加载端口 uop,但 AMD CPU 确实也需要 ALU。)哦,顺便说一句,你为什么使用 size_t
作为 lastBits
?它是带有 -m32
的 32 位类型,甚至是 -mx32
(带有 32 位指针的长模式)。 n
绝对是错误的类型。只需使用 unsigned
。
在一个相当不相关的说明上:更多的性能黑客!
[@ShreevatsaR 终于揭穿了第一个“猜想”;删除]
在遍历序列时,我们只能得到当前元素 N 的 2 邻域中的 3 种可能情况(先显示): [偶数] [奇数] [奇数] [偶数] [偶数] [偶数] 跳过这 2 个元素意味着分别计算 (N >> 1) + N + 1、((N << 1) + N + 1) >> 1 和 N >> 2。让我们证明对于情况 (1) 和 (2) 都可以使用第一个公式,(N >> 1) + N + 1。情况 (1) 是显而易见的。情况 (2) 意味着 (N & 1) == 1,因此如果我们假设(不失一般性)N 是 2 位长并且它的位从最高有效位到最低有效位是 ba,那么 a = 1,并且以下成立: (N << 1) + N + 1: (N >> 1) + N + 1: b10 b1 b1 b + 1 + 1 ---- --- bBb0 bBb 其中 B = !b。右移第一个结果给了我们想要的结果。 QED: (N & 1) == 1 ⇒ (N >> 1) + N + 1 == ((N << 1) + N + 1) >> 1. 正如证明的那样,我们可以遍历序列 2 个元素一次,使用单个三元运算。又减少了 2 倍的时间。
[偶数] [奇数]
[奇偶]
[偶数] [偶数]
生成的算法如下所示:
uint64_t sequence(uint64_t size, uint64_t *path) {
uint64_t n, i, c, maxi = 0, maxc = 0;
for (n = i = (size - 1) | 1; i > 2; n = i -= 2) {
c = 2;
while ((n = ((n & 3)? (n >> 1) + n + 1 : (n >> 2))) > 2)
c += 2;
if (n == 2)
c++;
if (c > maxc) {
maxi = i;
maxc = c;
}
}
*path = maxc;
return maxi;
}
int main() {
uint64_t maxi, maxc;
maxi = sequence(1000000, &maxc);
printf("%llu, %llu\n", maxi, maxc);
return 0;
}
这里我们比较 n > 2
,因为如果序列的总长度是奇数,则过程可能会在 2 而不是 1 处停止。
[编辑:]
让我们把它翻译成汇编!
MOV RCX, 1000000;
DEC RCX;
AND RCX, -2;
XOR RAX, RAX;
MOV RBX, RAX;
@main:
XOR RSI, RSI;
LEA RDI, [RCX + 1];
@loop:
ADD RSI, 2;
LEA RDX, [RDI + RDI*2 + 2];
SHR RDX, 1;
SHRD RDI, RDI, 2; ror rdi,2 would do the same thing
CMOVL RDI, RDX; Note that SHRD leaves OF = undefined with count>1, and this doesn't work on all CPUs.
CMOVS RDI, RDX;
CMP RDI, 2;
JA @loop;
LEA RDX, [RSI + 1];
CMOVE RSI, RDX;
CMP RAX, RSI;
CMOVB RAX, RSI;
CMOVB RBX, RCX;
SUB RCX, 2;
JA @main;
MOV RDI, RCX;
ADD RCX, 10;
PUSH RDI;
PUSH RCX;
@itoa:
XOR RDX, RDX;
DIV RCX;
ADD RDX, '0';
PUSH RDX;
TEST RAX, RAX;
JNE @itoa;
PUSH RCX;
LEA RAX, [RBX + 1];
TEST RBX, RBX;
MOV RBX, RDI;
JNE @itoa;
POP RCX;
INC RDI;
MOV RDX, RDI;
@outp:
MOV RSI, RSP;
MOV RAX, RDI;
SYSCALL;
POP RAX;
TEST RAX, RAX;
JNE @outp;
LEA RAX, [RDI + 59];
DEC RDI;
SYSCALL;
使用这些命令进行编译:
nasm -f elf64 file.asm
ld -o file file.o
请参阅 Peter Cordes on Godbolt 的 C 和 asm 的改进/错误修复版本。 (编者注:很抱歉把我的东西放在你的答案中,但我的答案达到了 Godbolt 链接 + 文本的 30k 字符限制!)
12 = 3Q + 1
这样的整数 Q
。你的第一点是不对的,我想。
mov reg, imm32
,显然是为了节省字节,但随后它到处使用 64 位版本的寄存器,即使是 xor rax, rax
,所以它有很多不必要的 REX 前缀。我们显然只需要在内循环中持有 n
的 reg 上的 REX 以避免溢出。
-O3 -march=core2
:96ms。 gcc5.2:108 毫秒。从我改进的 clang 的 asm 内循环版本:92ms(应该看到对 SnB 系列有更大的改进,其中复杂的 LEA 是 3c 而不是 1c)。从这个 asm 循环的改进 + 工作版本(使用 ROR+TEST,而不是 SHRD):87 毫秒。打印前以 5 次重复测量
在从源代码生成机器代码的过程中,C++ 程序被翻译成汇编程序。说汇编比 C++ 慢实际上是错误的。此外,生成的二进制代码因编译器而异。因此,一个聪明的 C++ 编译器可以生成比愚蠢的汇编程序代码更优化和更高效的二进制代码。
但是我相信你的分析方法有一定的缺陷。以下是概要分析的一般准则:
确保您的系统处于正常/空闲状态。停止所有正在运行的进程(应用程序)您已启动或密集使用 CPU(或通过网络轮询)。您的数据大小必须更大。您的测试必须运行超过 5-10 秒。不要只依赖一个样本。执行您的测试 N 次。收集结果并计算结果的平均值或中位数。
来自评论:
但是,这段代码永远不会停止(因为整数溢出)!?!伊夫·道斯特
对于许多数字,它不会溢出。
如果它会溢出 - 对于那些不幸的初始种子之一,溢出的数字很可能会收敛到 1 而不会再次溢出。
这仍然提出了一个有趣的问题,是否有一些溢出循环种子数?
任何简单的最终收敛级数都以两个值的幂开始(足够明显吗?)。
2^64 将溢出到零,根据算法这是未定义的无限循环(仅以 1 结束),但由于 shr rax
产生 ZF=1,答案中的最佳解决方案将完成。
我们可以生产 2^64 吗?如果起始编号为0x5555555555555555
,则为奇数,下一个编号为3n+1,即0xFFFFFFFFFFFFFFFF + 1
= 0
。理论上处于算法的未定义状态,但 johnfound 的优化答案将通过退出 ZF=1 来恢复。 Peter Cordes 的 cmp rax,1
将在无限循环中结束(QED 变体 1,“便宜”到未定义的 0
数字)。
一些更复杂的数字怎么样,它会在没有 0
的情况下创建循环?坦率地说,我不确定,我的数学理论太模糊了,无法得到任何严肃的想法,如何严肃地处理它。但直觉上我会说这个系列对于每个数字都会收敛到 1:0 <数,因为 3n+1 公式迟早会慢慢地将原始数(或中间数)的每个非 2 素数因子转换为 2 的某个幂。所以我们不需要担心原始系列的无限循环,只有溢出会阻碍我们。
所以我只是将几个数字放入工作表并查看了 8 位截断数字。
有三个值溢出到 0
:227
、170
和 85
(85
直接流向 0
,另外两个流向 85
)。
但是创建循环溢出种子没有任何价值。
有趣的是,我做了一个检查,这是第一个遭受 8 位截断的数字,并且 27
已经受到影响!它确实在适当的非截断序列中达到值 9232
(第 12 步中的第一个截断值是 322
),并且以非截断方式达到的 2-255 个输入数字中的任何一个的最大值是 13120
(对于 255
本身),收敛到 1
的最大步数约为 128
(+-2,不确定是否要计算“1”等...)。
有趣的是(对我来说)数字 9232
是许多其他源数字的最大值,它有什么特别之处? :-O 9232
= 0x2410
...嗯..不知道。
不幸的是,我无法深入了解这个系列,为什么它会收敛以及将它们截断为 k 位的含义是什么,但是使用 cmp number,1
终止条件当然可以将算法进入无限循环,特定输入值在截断后以 0
结尾。
但是 8 位情况下的值 27
溢出是一种警告,这看起来如果你计算达到值 1
的步数,你会从总 k 位集合中得到大多数数字的错误结果整数。对于 8 位整数,256 个中的 146 个数字通过截断影响了系列(其中一些可能仍然意外地达到正确的步数,我懒得检查)。
27
系列看起来像这样:82 41 124 62 31 94 47 142 71 214 107 66(截断)33 100 50 25 76 38 19 58 29 88 44 22 11 34 17 52 26 13 40 20 10 5 16 8 4 2 1 (其余部分没有截断)。我不明白你,对不起。如果截断值等于当前正在进行的系列中先前达到的一些值,它将永远不会停止,并且我找不到任何这样的值与 k 位截断(但我要么无法弄清楚背后的数学理论,为什么这适用于 8/16/32/64 位截断,直觉上我认为它有效)。
2
-255
数字,或者没有截断(到 1
) ,或使用 8 位截断(对于三个数字,为预期的 1
或 0
)。
cmp rax,1 / jna
(即 do{}while(n>1)
)也终止于零。我考虑制作一个记录最大 n
看到的循环的检测版本,以了解我们有多接近溢出。
您没有发布编译器生成的代码,所以这里有一些猜测,但即使没有看到它,也可以这样说:
test rax, 1
jpe even
... 有 50% 的机会错误预测分支,而且代价高昂。
编译器几乎肯定会进行这两种计算(由于 div/mod 的延迟非常长,因此乘加是“免费的”,因此成本几乎可以忽略不计)并跟进 CMOV。当然,错误预测的可能性为零。
对于 Collatz 问题,您可以通过缓存“尾部”来显着提高性能。这是时间/内存的权衡。请参阅:记忆 (https://en.wikipedia.org/wiki/Memoization)。您还可以研究其他时间/内存权衡的动态编程解决方案。
示例python实现:
import sys
inner_loop = 0
def collatz_sequence(N, cache):
global inner_loop
l = [ ]
stop = False
n = N
tails = [ ]
while not stop:
inner_loop += 1
tmp = n
l.append(n)
if n <= 1:
stop = True
elif n in cache:
stop = True
elif n % 2:
n = 3*n + 1
else:
n = n // 2
tails.append((tmp, len(l)))
for key, offset in tails:
if not key in cache:
cache[key] = l[offset:]
return l
def gen_sequence(l, cache):
for elem in l:
yield elem
if elem in cache:
yield from gen_sequence(cache[elem], cache)
raise StopIteration
if __name__ == "__main__":
le_cache = {}
for n in range(1, 4711, 5):
l = collatz_sequence(n, le_cache)
print("{}: {}".format(n, len(list(gen_sequence(l, le_cache)))))
print("inner_loop = {}".format(inner_loop))
0
条目表示尚不存在。我们可以通过仅在表中存储奇数 N 来进一步优化,因此哈希函数为 n>>1
,丢弃 1。编写步骤代码以始终以 n>>tzcnt(n)
或其他内容结尾以确保它是奇数。
作为一个通用答案,并非专门针对此任务:在许多情况下,您可以通过在高水平上进行改进来显着加快任何程序。就像计算数据一次而不是多次,完全避免不必要的工作,以最佳方式使用缓存等等。用高级语言来做这些事情要容易得多。
编写汇编代码,可以改进优化编译器的功能,但这是一项艰苦的工作。而且一旦完成,您的代码就更难修改,因此添加算法改进也更加困难。有时处理器具有您无法从高级语言使用的功能,内联汇编在这些情况下通常很有用,并且仍然允许您使用高级语言。
在欧拉问题中,大多数情况下,你成功的方式是构建一些东西,找出它为什么慢,构建更好的东西,找出它为什么慢,等等。使用汇编程序非常非常困难。以一半可能速度的更好算法通常会以全速击败更差的算法,并且在汇编程序中获得全速并非易事。
gcc -O3
针对该精确算法,在 Haswell 上生成的代码与最优值相差 20% 以内。 (获得这些加速是我回答的主要重点,只是因为这就是问题所问的,并且有一个有趣的答案,不是,因为它是正确的方法。)从编译器的转换中获得了更大的加速极不可能寻找,例如推迟右移或一次执行 2 步。比记忆化/查找表更大的加速。仍然是详尽的测试,但不是纯粹的蛮力。
即使不看汇编,最明显的原因是 /= 2
可能已优化为 >>=1
并且许多处理器具有非常快速的移位操作。但即使处理器没有移位操作,整数除法也比浮点除法快。
编辑:您的里程可能因上述“整数除法比浮点除法快”语句而异。下面的评论表明,现代处理器优先考虑优化 fp 除法而不是整数除法。因此,如果有人正在寻找该线程问题所询问的加速的最可能原因,那么编译器将 /=2
优化为 >>=1
将是最好的第一选择。
在不相关的注释中,如果 n
是奇数,则表达式 n*3+1
将始终是偶数。所以没有必要检查。您可以将该分支更改为
{
n = (n*3+1) >> 1;
count += 2;
}
所以整个声明将是
if (n & 1)
{
n = (n*3 + 1) >> 1;
count += 2;
}
else
{
n >>= 1;
++count;
}
DIV r32
(32 位无符号整数)或 DIV r64
(慢得多的 64 位无符号整数)进行比较。特别是对于吞吐量,FP 划分要快得多(单 uop 而不是微编码,并且部分流水线化),但延迟也更好。
div r64
是 36 uop、32-96c 延迟和每 21-74c 吞吐量一个。 Skylake 具有更快的 FP 除法吞吐量(以每 4c 一个流水线进行,延迟并没有好得多),但整数 div 的速度并没有快多少。 AMD Bulldozer 系列的情况类似:DIVSD 是 1M-op,9-27c 延迟,每 4.5-11c 吞吐量一个。 div r64
是 16M-ops,16-75c 延迟,每 16-75c 吞吐量一个。
double
具有 53 位尾数,但在 Haswell 上仍然比 div r32
慢得多。因此,这绝对只是英特尔/AMD 投入多少硬件来解决问题的问题,因为它们没有为整数和 fp 分频器使用相同的晶体管。整数 1 是标量(没有整数 SIMD 除法),向量 1 处理 128b 个向量(不像其他向量 ALU 那样处理 256b)。最重要的是整数 div 是很多 uops,对周围的代码影响很大。
简单的答案:
做一个 MOV RBX, 3 和 MUL RBX 很贵;只需添加 RBX,RBX 两次
ADD 1 在这里可能比 INC 快
MOV 2 和 DIV 非常昂贵;右移
64 位代码通常明显比 32 位代码慢,而且对齐问题更复杂;对于像这样的小程序,你必须打包它们,这样你就可以进行并行计算,以便有可能比 32 位代码更快
如果您为您的 C++ 程序生成程序集列表,您可以看到它与您的程序集有何不同。
mul rbx
为 2 微指令,延迟为 3c(每时钟吞吐量为 1)。 imul rcx, rbx, 3
仅为 1 uop,具有相同的 3c 延迟。两个 ADD 指令将是 2 微指令,具有 2c 延迟。
ADD RBX, RBX
两次将乘以 4,而不是 3)。到目前为止,最好的方法是 lea rax, [rbx + rbx*2]
。或者,以使其成为 3 分量 LEA 为代价,也可以使用 lea rax, [rbx + rbx*2 + 1]
进行 +1(HSW 上的延迟为 3c 而不是 1,正如我在回答中所解释的那样)我的观点是 64 位乘法不是在最近的 Intel CPU 上非常昂贵,因为它们具有非常快的整数乘法单元(甚至与 AMD 相比,其中相同的 MUL r64
是 6c 延迟,每 4c 吞吐量一个:甚至没有完全流水线化。
不定期副业成功案例分享
tzcnt
在标量代码中做得更好,并且您被锁定在矢量化中的矢量元素中运行时间最长的序列案子)。1
时,而不是当它们 all 有时(很容易用 PCMPEQ 检测到/PMOVMSK)。然后,您使用 PINSRQ 和其他东西来处理终止的一个元素(及其计数器),然后跳回到循环中。当您过于频繁地跳出内部循环时,这很容易变成一种损失,但这确实意味着您总是在每次内部循环迭代中完成 2 或 4 个有用的工作元素。不过,关于记忆的好点。