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


当前回答

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

其他回答

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

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

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

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

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

你没有发布编译器生成的代码,所以这里有一些猜测,但即使没有看到它,人们也可以说:

test rax, 1
jpe even

... 有50%的几率预测错分支,代价会很高昂。

编译器几乎肯定会进行这两种计算(因为div/mod的延迟相当长,所以乘法-加法是“免费的”),并随后进行CMOV。当然,被错误预测的可能性为零。

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