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


当前回答

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

其他回答

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

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

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

声称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毫秒vs. 1900毫秒),当用-O3编译时,快两倍多(430毫秒vs. 830毫秒)。

两个程序的输出是相同的:在i = 837799上,max sequence = 525。

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

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字符限制!)

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