我为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毫秒。


当前回答

作为一个通用的答案,并不是专门针对这个任务的:在许多情况下,您可以通过在较高的级别上进行改进来显著地加快任何程序。比如计算一次而不是多次,完全避免不必要的工作,以最佳方式使用缓存,等等。用高级语言做这些事情要容易得多。

编写汇编程序代码,可以改进优化编译器所做的工作,但这是一项艰苦的工作。一旦完成,你的代码就很难修改了,所以要添加算法改进就更加困难了。有时处理器具有高级语言无法使用的功能,内联汇编在这些情况下通常很有用,并且仍然允许您使用高级语言。

在欧拉问题中,大多数情况下你成功的方法是建立一个东西,找到它为什么慢,建立一个更好的东西,找到它为什么慢,等等。使用汇编程序是非常非常困难的。一个更好的算法在可能速度的一半时通常会击败一个更差的算法在全速时,在汇编程序中获得全速并不是一件简单的事情。

其他回答

为了获得更好的性能:一个简单的改变是观察到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真的很大。

在一个相当不相关的注意:更多的性能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字符限制!)

即使不考虑汇编,最明显的原因是/= 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;
}

作为一个通用的答案,并不是专门针对这个任务的:在许多情况下,您可以通过在较高的级别上进行改进来显著地加快任何程序。比如计算一次而不是多次,完全避免不必要的工作,以最佳方式使用缓存,等等。用高级语言做这些事情要容易得多。

编写汇编程序代码,可以改进优化编译器所做的工作,但这是一项艰苦的工作。一旦完成,你的代码就很难修改了,所以要添加算法改进就更加困难了。有时处理器具有高级语言无法使用的功能,内联汇编在这些情况下通常很有用,并且仍然允许您使用高级语言。

在欧拉问题中,大多数情况下你成功的方法是建立一个东西,找到它为什么慢,建立一个更好的东西,找到它为什么慢,等等。使用汇编程序是非常非常困难的。一个更好的算法在可能速度的一半时通常会击败一个更差的算法在全速时,在汇编程序中获得全速并不是一件简单的事情。

评论:

但是,这段代码从未停止(因为整数溢出)!?!伊夫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个数字通过截断影响了序列(其中一些可能仍然意外地命中了正确的步数,我太懒了,没有检查)。