python - 为什么以不同的方式生成数组会导致代码不同部分的大幅加速？

Question

我遇到过一个非常奇怪的案例，即在更改一个看似微不足道的细节后，代码会大幅提速。我有以下代码，它是在 SageMath 中实现的用于计算阶乘的 Borwein 算法的实现(但除了一些次要的东西，比如 ^ 用于求幂，它应该在纯 Python 2.7 中工作)

def sieve(n): #My implementation of the sieve of Eratosthenes
    T = [1]*n
    for i in xrange(2,n):
        if T[i]==1:
            for j in xrange(2,ceil(n/i)):
                T[i*j]=0
    return [i for i in xrange(2,n) if T[i]]
def expp(n,p): #Exponent of p in the factorization of n!
    k = p
    s = 0
    while k<=n:
        s += n//k
        k = k*p
    return s
def quick_prod(T): #Computing product of the elements of an array using binary splitting
    if len(T)==1:
        return T[0]
    if len(T)==2:
        return T[0]*T[1]
    if len(T)>2:
        s = len(T)//2
        return quick_prod(T[0:s])*quick_prod(T[s:len(T)])

n = 10^6
P = sieve(n) #Array of primes up to n
exps = [expp(n,p) for p in P] #exponents of all primes in P
l = len(bin(abs(n)))-2
nums = [quick_prod([P[j] for j in xrange(len(P)) if (exps[j] >> i)%2])^(2^i) for i in range(l)] #Array of numbers appearing in Borwein's algorithm, whose product is n!
quick_prod(nums)

(请原谅我糟糕的命名约定(可能还有其他糟糕的编码习惯)，我是一个业余爱好者，只真正编写“快速而肮脏”的东西)

没想到这段代码效率特别高，所以看到运行了10分钟我并不意外。但是当我开始修改代码以尝试改进它时，我注意到将行 P = sieve(n) 替换为 P = prime_range(n)(生成相同的数组，除了它使用 SageMath 中内置的函数)将运行时间减少到 3.5 秒。

现在，当我看到这个时，我的第一个想法是解释很明显 - 我对筛子的实现一定很糟糕，它花了很长时间，而 prime_range 这样做效率更高。但结果让我大吃一惊——sieve(10^6) 花了 4 秒，而 prime_range(10^6) 花了 2 秒。这甚至无法解释 10 分钟的差异!

我和我 friend 的一些想法可能可以解释它:

• 这两个数组可能不同，例如它们的顺序可能不同。情况并非如此，因为 sieve(10^6)==prime_range(10^6) 返回 True
• 尽管相等，但数组可以存储为不同类型。 type(...) 为两者返回 list。
• 中间结果被缓存。可能不是这种情况，因为无论编译顺序如何，结果都大致相同，即使在重新启动内核之后也是如此。

如此巨大的加速(或减速，取决于你如何看待)的唯一方法是，如果原始代码以某种方式返回到 之后 生成 P 的方式计算它。什么可以解释这种行为？

Answer 1

SageMath(我不熟悉)可能像 numpy 一样工作。也就是说，它使用(和返回)的数组结构和内部数据类型比标准的 python 列表更高效。这可能适用于列表理解和您之后使用它执行的其他计算。

这是这种现象的一个例子(基于 numpy)。

import numpy as np
def sieve2(n):
    s       = np.ones(n+1)
    s[4::2] = 0
    s[:2]   = 0
    p = 3
    while p*p<=n:
        if s[p]:s[p*p::p] = 0
        p+=2
    return np.arange(n+1)[s==1]

此函数在 0.013 秒内返回最多 10^6 个素数，而你的函数在我的计算机上需要 0.27 秒(大约快 20 倍)。基于 numpy 的函数返回一个 numpy 数组，它有自己的基本函数实现，例如加法、乘法、求幂等。sageMath 也可能是这种情况，它可能会加速程序的其他部分。

请注意，使用 numpy 的巨大时间差异来自其向量化计算和使用 GPU 并行执行多个操作的能力。 SageMath 可能会对其大整数计算使用相同的技巧(Python 可能不会这样做)