协慌网

登录 贡献 社区

用于测试 Collatz 猜想的 C ++ 代码比手写程序集更快 - 为什么?

我为Project Euler Q14编写了这两个解决方案,在汇编和 C ++ 中。它们是用于测试Collatz 猜想的相同蛮力方法。装配解决方案与组装

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>

using namespace std;

int sequence(long n) {
    int count = 1;
    while (n != 1) {
        if (n % 2 == 0)
            n /= 2;
        else
            n = n*3 + 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;
        }
    }

    cout << maxi << endl;
}

我知道编译器优化以提高速度和一切,但我没有看到很多方法来进一步优化我的程序集解决方案(以编程方式而不是数学方式)。

C ++ 代码在每个术语和每个偶数项的除法中都有模数,其中汇编每个偶数项只有一个除法。

但是组装平均比 C ++ 解决方案长 1 秒。为什么是这样?我主要是好奇地问。

执行时间

我的系统:1.4 GHz Linux 1.4 GHz Intel Celeron 2955U(Haswell 微体系结构)。

答案

如果您认为 64 位 DIV 指令是除以 2 的好方法,那么难怪编译器的 asm 输出会超过您的手写代码,即使使用-O0 (编译速度快,无需额外优化,并存储 / 重新加载到内存中)在每个 C 语句之后 / 之前,因此调试器可以修改变量)。

请参阅Agner Fog 的优化装配指南 ,了解如何编写高效的 asm。他还有针对特定 CPU 的具体细节的指令表和微指南指南。有关更多 perf 链接,另请参阅标记 wiki。

另请参阅关于使用手写 asm 击败编译器的这个更一般的问题: 内联汇编语言是否比本机 C ++ 代码慢? 。 TL:DR:是的,如果你做错了(就像这个问题)。

通常你很好让编译器做它的事情,特别是如果你试图编写可以有效编译的 C ++ 。还看到汇编比编译语言快吗? 。其中一个答案链接到这些简洁的幻灯片,显示各种 C 编译器如何使用很酷的技巧优化一些非常简单的功能。


even:
    mov rbx, 2
    xor rdx, rdx
    div rbx

在 Intel Haswell 上, div r64为 36 div r64延迟为 32-96 个周期 ,吞吐量为每 21-74 个周期一个。 (加上 2 uop 来设置 RBX 和零 RDX,但是乱序执行可以提前运行)。 像 DIV 这样的高 uop 计数指令是微编码的,这也可能导致前端瓶颈。在这种情况下,延迟是最相关的因素,因为它是循环携带的依赖链的一部分。

shr rax, 1执行相同的无符号除法:它是 1 uop,具有 1c 延迟 ,并且每个时钟周期可以运行 2。

相比之下,32 位除法速度更快,但对于换档仍然很可怕。 idiv r32是 9 uops,22-29c 延迟,Haswell 每 8-11c 吞吐量一个。


正如您从 gcc 的-O0 asm 输出( Godbolt 编译器浏览器 )中看到的那样 ,它只使用移位指令 。 clang -O0确实像你想象的那样天真地编译,甚至使用 64 位 IDIV 两次。 (当进行优化时,编译器会在源进行除法时使用 IDIV 的两个输出,并使用相同的操作数使用模数,如果它们完全使用 IDIV 的话)

海湾合作委员会没有完全天真的模式; 它总是通过 GIMPLE 进行转换,这意味着无法禁用某些 “优化” 。这包括识别div_by_13除法和使用移位(2 的幂)或固定点乘法逆 (2 的非幂)以避免 IDIV(参见上面的 godbolt 链接中的 div_by_13)。

gcc -Os (针对大小进行优化) 确实使用 IDIV 进行非 2 次幂除法,遗憾的是,即使在乘法逆码仅略大但速度更快的情况下也是如此。


帮助编译器

(本案例摘要:使用uint64_t n

首先,看一下优化的编译器输出是很有意思的。 ( -O3 )。 -O0速度基本上没有意义。

看看你的 asm 输出(在 Godbolt 上,或者看看如何从 GCC / clang 组件输出中删除 “噪音”? )。当编译器首先没有制作最佳代码时: 以引导编译器制作更好代码的方式编写 C / C ++ 源代码通常是最好的方法 。你必须知道 asm,并知道什么是有效的,但你间接地应用这些知识。编译器也是一个很好的想法来源:有时 clang 会做一些很酷的事情,你可以手动执行 gcc 做同样的事情:看看这个答案以及我在 @ Veedrac 的代码中对非展开循环所做的事情。)

这种方法是可移植的,并且在未来 20 年内,一些未来的编译器可以将其编译为在未来硬件(x86 或不是)上有效的任何东西,可能使用新的 ISA 扩展或自动向量化。从 15 年前手写的 x86-64 asm 通常不会为 Skylake 进行最佳调整。例如,比较和分支宏观融合当时不存在。 对于一个微体系结构,手工制作的 asm 现在最理想的可能不是其他当前和未来 CPU 的最佳选择。 关于 @ johnfound 的回答的评论讨论了 AMD Bulldozer 和 Intel Haswell 之间的主要差异,这对这段代码产生了很大的影响。但理论上, g++ -O3 -march=bdver3g++ -O3 -march=skylake会做正确的事情。 (或者-march=native 。)或者-mtune=...只是调整,而不使用其他 CPU 可能不支持的指令。

我的感觉是指导编译器 asm 对你关心的当前 CPU 有好处,对于未来的编译器来说应该不是问题。他们希望在寻找转换代码的方法方面比当前的编译器更好,并且可以找到适用于未来 CPU 的方法。无论如何,未来的 x86 在当前 x86 上的任何好处都可能不会很糟糕,并且未来的编译器将避免任何特定于 asm 的陷阱,同时实现类似于来自 C 源的数据移动,如果它没有看到更好的东西。

手写 asm 是优化器的黑盒子,因此当内联使输入成为编译时常量时,常量传播不起作用。其他优化也会受到影响。在使用 asm 之前,请阅读https://gcc.gnu.org/wiki/DontUseInlineAsm 。 (并避免使用 MSVC 样式的内联 asm:输入 / 输出必须通过内存,这会增加开销 。)

在这种情况下 :你的n有一个带符号的类型,gcc 使用 SAR / SHR / ADD 序列给出正确的舍入。 (对于负输入,IDIV 和算术移位 “舍入” 不同,请参阅SAR insn set ref 手册输入 )。 (IDK 如果 gcc 尝试并且未能证明n不能为负,或者是什么。签名溢出是未定义的行为,所以应该能够。)

你应该使用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 个周期,延迟瓶颈 。乱序执行与此并行处理其他所有事情(理论上:我没有使用 perf 计数器进行测试,看它是否真的以 5c / iter 运行)。

cmov的 FLAGS 输入(由 TEST 生成)比 RAX 输入(来自 LEA-> MOV)生成得更快,因此它不在关键路径上。

类似地,产生 CMOV 的 RDI 输入的 MOV-> SHR 不在关键路径上,因为它也比 LEA 快。 IvyBridge 上的 MOV 以及之后的延迟为零(在寄存器重命名时间处理)。 (它仍然需要一个 uop 和一个插槽,因此它不是免费的,只是零延迟)。 LEA dep 链中的额外 MOV 是其他 CPU 瓶颈的一部分。

cmp / jne 也不是关键路径的一部分:它不是循环携带的,因为控制依赖性是通过分支预测 + 推测执行来处理的,这与关键路径上的数据依赖性不同。


打败编译器

海湾合作委员会在这里做得很好。它可以通过使用inc edx而不是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)时才为零。 (有趣的事实: 如果你读取标志结果,那么 Nehalem 或更早的 SHR 会计数!= 1 会导致失速 。这就是他们如何使它成为单 uop。但是,按 1 转换特殊编码很好。)

避免 MOV 对 Haswell 的延迟没有帮助( Can x86 的 MOV 真的可以 “免费” 吗?为什么我根本不能重现这个? )。它对英特尔 pre-IvB 和 AMD Bulldozer 系列等 CPU 有很大帮助,其中 MOV 不是零延迟。编译器浪费的 MOV 指令会影响关键路径。 BD 的复杂 LEA 和 CMOV 都具有较低的延迟(分别为 2c 和 1c),因此它是延迟的一小部分。此外,吞吐量瓶颈成为一个问题,因为它只有两个整数 ALU 管道。 请参阅 @ johnfound 的答案 ,他有 AMD CPU 的计时结果。

即使在 Haswell 上,这个版本可能会有所帮助,避免一些偶然的延迟,其中非关键 uop 从关键路径上的一个执行端口窃取执行端口,将执行延迟 1 个周期。 (这称为资源冲突)。它还保存了一个寄存器,这可能有助于在交错循环中并行执行多个n值(见下文)。

LEA 的延迟取决于 Intel SnB 系列 CPU 上的寻址模式 。 3c 为 3 个组件( [base+idx+const] ,需要两个单独的添加),但只有 1c 有 2 个或更少的组件(一个添加)。有些 CPU(如 Core2)在一个周期内甚至可以执行 3 分量 LEA,但 SnB 系列则不行。更糟糕的是, 英特尔 SnB 系列标准化延迟,因此没有 2c uops ,否则 3 组件 LEA 将只有 2c 像 Bulldozer。 (3 组件 LEA 在 AMD 上也较慢,但没有那么多)。

所以lea rcx, [rax + rax*2] / inc rcx在像 Haswell 这样的 Intel SnB 系列 CPU 上只有 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 。这创建了更多的指令级并行性

诀窍是决定是否要等到所有n值都达到1才能获得另一对起始n值,或者是否突破并为只有达到结束条件的一个获得新的起点,而不触及寄存器其他顺序。可能最好让每个链都处理有用的数据,否则你必须有条件地增加它的计数器。


您甚至可以使用 SSE 打包比较的东西来有条件地增加n尚未达到1向量元素的计数器。然后为了隐藏 SIMD 条件增量实现的更长延迟,您需要在空中保留更多n值向量。也许只值得 256b 矢量(4x uint64_t )。

我认为检测1 “粘性” 的最佳策略是屏蔽你添加的所有 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?

    vprllq    ymm3, ymm0, 63                # shift bit 1 to the sign bit

    vpsrlq    ymm0, ymm0, 1                 # n /= 2

    # There may be a better way to do this blend, avoiding the bypass delay for an FP blend between integer insns, not sure.  Probably worth it
    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 s。

所以循环可能如下所示:

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 的 Intel SnB 系列 CPU 上,可变计数移位速度很慢。 3 次 uops,2c 延迟。 (它们对 FLAGS 有输入依赖性,因为 count = 0 表示标志未经修改。它们将此处理为数据依赖性,并且因为 uop 只能有 2 个输入(无论如何都是 HSW / BDW),因此需要多次 uop)。这是人们抱怨 x86 的疯狂 CISC 设计所指的那种。它使得 x86 CPU 比现在从头开始设计 ISA 时的速度慢,即使是以大致相似的方式。 (即这是成本速度 / 功率的 “x86 税” 的一部分。)SHRX / SHLX / SARX(BMI2)是一个巨大的胜利(1 uop / 1c 延迟)。

它还在关键路径上放置 tzcnt(Haswell 及更高版本的 3c),因此它显着延长了循环携带依赖链的总延迟。但它确实消除了对 CMOV 的任何需要,或者准备了一个持有n>>1的寄存器。 @ Veedrac 的答案通过延迟 tzcnt / shift 进行多次迭代来克服所有这些,这非常有效(见下文)。

我们可以安全地交替使用BSFTZCNT ,因为n在此时永远不会为零。 TZCNT 的机器代码在不支持 BMI1 的 CPU 上解码为 BSF。 (忽略无意义的前缀,因此 REP BSF 作为 BSF 运行)。

TZCNT 在支持它的 AMD CPU 上比 BSF 表现要好得多,所以使用REP BSF也是个好主意,即使输入为零而不是输出也不关心设置 ZF。当你使用__builtin_ctzll甚至使用-mno-bmi时,一些编译器会这样做。

它们在 Intel CPU 上执行相同的操作,因此只需保存字节,如果这一切都很重要。英特尔(pre-Skylake)上的 TZCNT 仍然假定依赖于所谓的只写输出操作数,就像 BSF 一样,以支持输入 = 0 的 BSF 未修改的 BSF 的未记录行为。因此,除非仅针对 Skylake 进行优化,否则您需要解决这个问题,因此额外的 REP 字节无法获得任何好处。 (英特尔经常超出 x86 ISA 手册所要求的范围,以避免破坏广泛使用的代码,这些代码依赖于它不应该的东西,或者是追溯不允许的。例如, Windows 9x 假定没有推测预取 TLB 条目 ,这是安全的代码编写时, 在英特尔更新 TLB 管理规则之前 。)

无论如何,Haswell 的 LZCNT / TZCNT 与 POPCNT 具有相同的假设:请参阅此问答 。这就是为什么在 @ Veedrac 的代码的 gcc 的 asm 输出中,当你不使用 dst = src 时,你会看到它在将要用作 TZCNT 的目标的寄存器上用 xor- zero 打破 dep 链 。由于 TZCNT / LZCNT / POPCNT 永远不会使其目标未定义或未修改,因此对 Intel CPU 上的输出的这种错误依赖纯粹是性能错误 / 限制。据推测,值得一些晶体管 / 功率使它们像其他 uops 一样运行到同一个执行单元。唯一的软件可见上行是与另一个微体系结构限制的交互: 它们可以在 Haswell 上使用索引寻址模式内存操作数进行微融合 ,但在 Skylake 上,英特尔取消了对 LZCNT / TZCNT 的错误依赖,他们 “取消层压” 索引寻址模式,而 POPCNT 仍然可以微融合任何地址模式。


从其他答案改进想法 / 代码:

@ hidefromkgb 的答案有一个很好的观察,你保证能够在 3n + 1 之后做一次右移。您可以更有效地计算它,而不仅仅是在步骤之间省略检查。然而,该答案中的 asm 实现被破坏了(它取决于 OF,在 SHRD 之后未定义,计数 > 1),而且慢: ROR rdi,2SHRD rdi,rdi,2快,并且使用两个 CMOV 指令在关键路径上比可以并行运行的额外 TEST 慢。

我整理 / 改进了 C(它指导编译器产生更好的 asm),并在 Godbolt 上测试 + 工作更快的 asm(在 C 下面的注释中):请参阅@ hidefromkgb 的答案中的链接。 (这个答案达到了大型 Godbolt 网址的 30k char 限制,但是短链接可能会腐烂 ,而且对于 goo.gl 来说太长了。)

还改进了输出打印以转换为字符串并进行一次write()而不是一次编写一个 char。使用perf stat ./collatz (记录性能计数器)可以最大限度地减少对整个程序计时的影响,并且我对一些非关键的 asm 进行了去混淆。


@ Veedrac 的代码

我从正确的移位中得到了一个非常小的加速,我们知道需要做的事情,并检查继续循环。从 CoresDuo(Merom)的 7.5s for limit = 1e8 下降到 7.275s,展开因子为 16。

关于 Godbolt 的代码 + 评论。请勿将此版本与 clang 一起使用; 它使用延迟循环做一些愚蠢的事情。使用 tmp 计数器k ,然后将其添加到count以后更改 clang 的作用,但这会轻微伤害 gcc。

请参阅评论中的讨论:Veedrac 的代码在具有 BMI1(即不是赛扬 / 奔腾)的 CPU 上非常出色

声称 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 与 1900 ms 编译时),速度提高了两倍多( - 430 ms 与 830 ms)使用-O3编译 C ++ 代码时。

两个程序的输出相同:i = 837799 时最大序列 = 525。

为了获得更好的性能:一个简单的变化是观察到在 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 次完全取决于这 8 位。例如,如果最后 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,少于 8 个除以 2,并且一次做 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 高速缓存时,我表示速度减慢。

PPS。如果你需要操作次数:在每次迭代中,我们完成八个除以 2 和一个可变数量的(3n + 1)操作,因此计算操作的一个明显方法是另一个数组。但我们实际上可以计算步数(基于循环的迭代次数)。

我们可以稍微重新定义问题:如果奇数将 n 替换为(3n + 1)/ 2,如果是偶数则用 n / 2 替换 n。然后每次迭代都会完成 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 真的很大。