python - 为什么两种寻找素数的算法在速度上相差如此之大，即使它们看起来迭代次数相同？

Question

我在 Python 中有两种寻找素数的算法。每一个的内循环看起来都执行了相同的次数，同样简单。然而，其中一个比另一个花费了 10 倍。我的问题是:

Why? Is this some quirk of Python that can be optimized away (how?), or am I missing something else?

我正在解决的问题本质上来自http://www.spoj.pl/problems/PRIME1/ .在我的例子中，我有 N = 10 ** 9，delta = 10 ** 5，我想找到 N-delta 和 delta 之间的所有素数。我还有 smallprimes，这是一个预制列表，包含所有小于或等于 N 的平方根的素数。第一个算法非常简单:

def range_f1(lo, hi, smallprimes):
  """Finds all primes p with lo <= p <= hi. 

  smallprimes is the sorted list of all primes up to (at least) square root of hi.
  hi & lo might be large, but hi-lo+1 miust fit into a long."""

  primes =[]
  for i in xrange(hi-lo+1):
    n = lo + i

    isprime = True
    for p in smallprimes:
        if n % p == 0:
            isprime = False
            break

    if isprime:
        primes.append(n)

  return primes

调用 range_f1(N-delta,N,smallprimes) 需要很长时间(大约 10 秒)。内循环调用了 195170 次。我也有这个算法的一个版本，它用列表理解替换循环(这是我实际用于分析的那个；见问题的结尾)但这并不快。

第二个版本是埃拉托色尼筛法的(一个丑陋的实现):

def range_sieve(lo, hi, smallprimes):
  """Parameters as before"""

  # So ugly! But SO FAST!! How??

  delta = hi-lo+1
  iamprime = [True] * delta      # iamprime[i] says whether lo + i is prime
  if lo<= 1:
    iamprime[1 - lo] = False

  def sillyfun():      # For profiling & speed comparison
    pass

  for p in smallprimes:
    rem = lo % p
    if rem == 0:
        iamprime[0] = False
    for i in xrange(p - rem, delta, p):
        iamprime[i] = False
        sillyfun()

    if p >= lo and p <= hi:
        iamprime[p - lo] = True

  return [p + lo for (p, ami) in enumerate(iamprime) if ami]

这大约快了 10 倍，用时不到 2 秒。但是，内部循环 (sillyfun()) 被调用了 259982 次，比最后一个案例多。我无法解释为什么这么快。

我想可能是因为第一个算法的内循环包含模运算，而第二个只有一个赋值。然而，以下似乎暗示赋值并不比模运算快:

>>> from timeit import timeit
>>> timeit("10**9 % 2341234")
0.23445186469234613
>>> timeit("a[5000]=False", "a = [True] * 10000")
0.47924750212666822

---

这是我实际使用的第一个算法的(可读性较差)版本:

def range_f2(lo, hi, smallprimes):

  primes =[]
  for i in xrange(hi-lo+1):
    n = lo + i

    try:
        (1 for p in smallprimes if n % p ==0).next()
    except StopIteration:
        primes.append(n)

  return primes

以下是为 range_f2() 调用探查器的结果(注意计算生成表达式的时间数):

>>> from cProfile import run as prof
>>> prof("range_f2(N-delta,N,sp)")
 200005 function calls in 13.866 CPU seconds

 Ordered by: standard name

 ncalls  tottime  percall  cumtime  percall filename:lineno(function)
      1    0.000    0.000   13.866   13.866 <string>:1(<module>)
 195170   12.632    0.000   12.632    0.000 prime1.py:101(<genexpr>)
      1    1.224    1.224   13.865   13.865 prime1.py:90(range_f2)
   4832    0.009    0.000    0.009    0.000 {method 'append' of 'list' objects}
      1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

这是 range_sieve() 的分析器结果:

>>> prof("range_sieve(N-delta,N,sp)")
259985 function calls in 1.370 CPU seconds

Ordered by: standard name
ncalls  tottime  percall  cumtime  percall filename:lineno(function)
     1    0.003    0.003    1.370    1.370 <string>:1(<module>)
     1    0.877    0.877    1.367    1.367 prime1.py:39(range_sieve)
259982    0.490    0.000    0.490    0.000 prime1.py:51(sillyfun)
     1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

最后，这是他生成小素数列表的完整代码(以一种非常愚蠢的方式)，以便您可以检查得到的结果:http://pastebin.com/7sfN4mG4

更新 根据大众需求，第一段代码的分析数据。没有关于内部循环执行了多少次的数据，但很明显它与第三次相同。

>>> prof("range_f1(N-delta,N,sp)")
      4835 function calls in 14.184 CPU seconds
Ordered by: standard name
ncalls  tottime  percall  cumtime  percall filename:lineno(function)
     1    0.000    0.000   14.184   14.184 <string>:1(<module>)
     1   14.162   14.162   14.183   14.183 prime1.py:69(range_f1)
  4832    0.021    0.000    0.021    0.000 {method 'append' of 'list' objects}
     1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

Answer 1

区别在于算法。在第一个版本中，trial division，你针对所有小素数测试每个候选 - 当小素数超过 candidate ** 0.5 时你不会停止对于 range( 10**9 - 10**5 ,10**9) 如果 smallprimes 有一个很好的上限，但如果范围的长度相对于上限大得多。对于复合 Material ，这不会产生太多成本，因为它们中的大多数至少有一个小素因数。但是对于素数，你必须一路走到 N**0.5。那个区间内大概有10**5/log(10**9)个素数，每一个素数都被试分为10**4.5/log(10**4.5 ) 素数，因此对素数进行了大约 1.47*10**7 次试验。

另一方面，使用筛子，您只能划掉组合，每个组合被划掉的次数与它有素数的次数一样多，而素数根本不会划掉(因此素数是免费的)。 n 的质因数的数量受限于(的倍数)log(n)(这是一个粗略的上限，通常大大高估)，因此给出了上限10**5*log(10**9)(乘以一个小常数)交叉的边界，大约 2*10**6。除此之外，交叉可能比除法更少工作(不知道 Python，它用于 C 数组)。所以你用筛子做的工作更少。

编辑:收集了 10**9-10**5 到 10**9 的实际数字。

Ticks: 259987
4832
Divisions: 20353799
4832

筛子只进行了 259987 次交叉，您会看到上面的粗略上限高估了一个很大的因素。试验除法需要超过 2000 万个除法，其中 16433632 个用于素数(x/log x 低估了素数的数量，对于 x = 10**4.5 大约低估了 10 %), 3435183 用于该范围内最小质因数大于n**(1/3)的3297个数。