我为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>

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上的64位Linux (Haswell微架构)。

g++(未优化):平均1272毫秒。 g++ -O3:平均578毫秒。 Asm (div)(原始):平均2650毫秒。 Asm (shr):平均679 ms。 @johnfound asm(与NASM组装):平均501毫秒。 @hidefromkgb asm:平均200毫秒。 @hidefromkgb asm,由@Peter Cordes优化:平均145毫秒。 @Veedrac c++: -O3平均81毫秒,-O0平均305毫秒。


当前回答

对于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))

其他回答

评论:

但是,这段代码从未停止(因为整数溢出)!?!伊夫Daoust

对于许多数字,它不会溢出。

如果它会溢出——对于那些不幸的初始种子之一,溢出的数字很可能会收敛到1而不会再次溢出。

这仍然提出了一个有趣的问题,是否存在一些溢出循环的种子数?

任何简单的最终收敛级数都以2的幂值开始(够明显了吗?)

2^64会溢出到零,根据算法,这是一个未定义的无限循环(仅以1结束),但由于shr rax产生ZF=1, answer中的最优解将结束。

能得到2^64吗?如果起始数是0x5555555555555555,它是奇数,那么下一个数是3n+1,它是0xFFFFFFFFFFFFFFFF +1 = 0。理论上算法处于未定义状态,但johnfound的优化答案在ZF=1时退出即可恢复。cmp rax, Peter Cordes的1将在无限循环中结束(QED变体1,“cheapo”通过未定义的0数)。

How about some more complex number, which will create cycle without 0? Frankly, I'm not sure, my Math theory is too hazy to get any serious idea, how to deal with it in serious way. But intuitively I would say the series will converge to 1 for every number : 0 < number, as the 3n+1 formula will slowly turn every non-2 prime factor of original number (or intermediate) into some power of 2, sooner or later. So we don't need to worry about infinite loop for original series, only overflow can hamper us.

我只是把一些数字写在表格里,然后看了看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数,1终止条件,当然可以将算法放入无限循环,特定的输入值在截断后以0结束。

但是对于8位的情况,值27溢出有点警告,这看起来像如果你计算到达值1的步数,你会得到错误的结果,对于整个k位整数集的大多数数字。对于8位整数,256个数字中的146个数字通过截断影响了序列(其中一些可能仍然意外地命中了正确的步数,我太懒了,没有检查)。

c++程序在从源代码生成机器码的过程中被转换为汇编程序。说汇编比c++慢实际上是错误的。此外,生成的二进制代码因编译器而异。因此,一个聪明的c++编译器可能会生成比愚蠢的汇编器代码更优、更有效的二进制代码。

但我认为你的分析方法有一定缺陷。以下是概要分析的一般准则:

确保您的系统处于正常/空闲状态。停止所有已启动或大量使用CPU的正在运行的进程(应用程序)(或在网络上轮询)。 您的数据必须更大。 您的测试必须运行5-10秒以上。 不要仅仅依赖于一个样本。执行测试N次。收集结果并计算结果的平均值或中位数。

在一个相当不相关的注意:更多的性能hack !

[the first «conjecture» has been finally debunked by @ShreevatsaR; removed] When traversing the sequence, we can only get 3 possible cases in the 2-neighborhood of the current element N (shown first): [even] [odd] [odd] [even] [even] [even] To leap past these 2 elements means to compute (N >> 1) + N + 1, ((N << 1) + N + 1) >> 1 and N >> 2, respectively. Let`s prove that for both cases (1) and (2) it is possible to use the first formula, (N >> 1) + N + 1. Case (1) is obvious. Case (2) implies (N & 1) == 1, so if we assume (without loss of generality) that N is 2-bit long and its bits are ba from most- to least-significant, then a = 1, and the following holds: (N << 1) + N + 1: (N >> 1) + N + 1: b10 b1 b1 b + 1 + 1 ---- --- bBb0 bBb where B = !b. Right-shifting the first result gives us exactly what we want. Q.E.D.: (N & 1) == 1 ⇒ (N >> 1) + N + 1 == ((N << 1) + N + 1) >> 1. As proven, we can traverse the sequence 2 elements at a time, using a single ternary operation. Another 2× time reduction.

得到的算法如下所示:

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

请看C和由Peter Cordes在Godbolt上编写的asm的改进/bug修复版本。(编者注:抱歉把我的东西放在你的答案中,但我的答案达到了Godbolt链接+文本的30k字符限制!)

为了获得更好的性能:一个简单的改变是观察到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缓存时,我预计会放缓。

pp。如果你需要运算的次数:在每次迭代中,我们做了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真的很大。

答案很简单:

做一个MOV RBX, 3和MUL RBX是昂贵的;只要加上RBX, RBX两次 add1在这里可能比INC更快 MOV 2和DIV非常昂贵;向右平移 64位代码通常明显比32位代码慢,对齐问题更复杂;对于这样的小程序,你必须打包它们,这样你才能进行并行计算,从而有可能比32位代码更快

如果为c++程序生成程序集列表,可以看到它与程序集的区别。