クックブックのerat2関数はさらに高速化できます (約 20 ～ 25%)。
import itertools as it
def erat2a( ):
D = { }
yield 2
for q in it.islice(it.count(3), 0, None, 2):
p = D.pop(q, None)
if p is None:
D[q*q] = q
yield q
else:
# old code here:
# x = p + q
# while x in D or not (x&1):
# x += p
# changed into:
x = q + 2*p
while x in D:
x += 2*p
D[x] = p
チェックはそれが奇数not (x&1)であることを確認します。xただし、 との両方 が奇数であるため、手順の半分を追加することで、奇数のテストとともに回避されます。qp2*p
少し余分な空想を気にしない場合はerat2、次の変更で 35-40% 高速化できます (注意: 関数のために Python 2.7+ または Python 3+ が必要ですitertools.compress):
import itertools as it
def erat3( ):
D = { 9: 3, 25: 5 }
yield 2
yield 3
yield 5
MASK= 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0,
MODULOS= frozenset( (1, 7, 11, 13, 17, 19, 23, 29) )
for q in it.compress(
it.islice(it.count(7), 0, None, 2),
it.cycle(MASK)):
p = D.pop(q, None)
if p is None:
D[q*q] = q
yield q
else:
x = q + 2*p
while x in D or (x%30) not in MODULOS:
x += 2*p
D[x] = p
このerat3関数は、30 を法とするすべての素数 (2、3、5 を除く) が 8 つの数 ( MODULOSfrozenset に含まれる数) になるという事実を利用しています。したがって、最初の 3 つの素数が得られた後、7 から始めて候補のみを処理します。
候補フィルタリングはitertools.compress関数を使用します。「魔法」はMASKシーケンスにあります。15 の要素 (関数MASKによって選択された 30 の数字ごとに 15 の奇数がある) を持ち、 7 から始まるすべての可能な候補に対して。サイクルは関数によって指定されたとおりに繰り返されます。
候補フィルタリングの導入には、チェックという別の変更が必要です。のitertools.islice1itertools.cycleor (x%30) not in MODULOSerat2アルゴリズムはすべての奇数を処理しました。erat3アルゴリズムが r30 候補のみを処理するようになったので、すべてがそのような —false — 候補のみになる可能性があることを確認する必要がありますD.keys()。
Atom 330 Ubuntu 9.10 サーバー、バージョン 2.6.4 および 3.1.1+:
$ testit
up to 8192
==== python2 erat2 ====
100 loops, best of 3: 18.6 msec per loop
==== python2 erat2a ====
100 loops, best of 3: 14.5 msec per loop
==== python2 erat3 ====
Traceback (most recent call last):
…
AttributeError: 'module' object has no attribute 'compress'
==== python3 erat2 ====
100 loops, best of 3: 19.2 msec per loop
==== python3 erat2a ====
100 loops, best of 3: 14.1 msec per loop
==== python3 erat3 ====
100 loops, best of 3: 11.7 msec per loop
AMD Geode LX Gentoo ホーム サーバー、Python 2.6.5 および 3.1.2:
$ testit
up to 8192
==== python2 erat2 ====
10 loops, best of 3: 104 msec per loop
==== python2 erat2a ====
10 loops, best of 3: 81 msec per loop
==== python2 erat3 ====
Traceback (most recent call last):
…
AttributeError: 'module' object has no attribute 'compress'
==== python3 erat2 ====
10 loops, best of 3: 116 msec per loop
==== python3 erat2a ====
10 loops, best of 3: 82 msec per loop
==== python3 erat3 ====
10 loops, best of 3: 66 msec per loop
primegen.pyモジュールには、、および関数がerat2含まれます。テストスクリプトは次のとおりです。erat2aerat3
#!/bin/sh
max_num=${1:-8192}
echo up to $max_num
for python_version in python2 python3
do
for function in erat2 erat2a erat3
do
echo "==== $python_version $function ===="
$python_version -O -m timeit -c \
-s "import itertools as it, functools as ft, operator as op, primegen; cmp= ft.partial(op.ge, $max_num)" \
"next(it.dropwhile(cmp, primegen.$function()))"
done
done
OP は効率的な実装を求めているため、 David Eppstein /Alex Martelli によるアクティブ ステート 2002 コードの大幅な改善を次に示します (ここで彼の回答を参照)。候補。生成される n 個の素数 ( π(sqrt(n log n)) ~ 2 sqrt ( n log n) / log(n log n) ~ 2 sqrt(n / log n) )。その結果、時間の複雑さも改善されます。つまり、実行速度が速くなります。
各基本素数の現在の倍数 (つまり、現在の生産点の sqrt より下) の辞書として「スライディングふるい」を作成し、それらのステップ値を一緒に作成します。
from itertools import count
# ideone.com/aVndFM
def postponed_sieve(): # postponed sieve, by Will Ness
yield 2; yield 3; yield 5; yield 7; # original code David Eppstein,
sieve = {} # Alex Martelli, ActiveState Recipe 2002
ps = postponed_sieve() # a separate base Primes Supply:
p = next(ps) and next(ps) # (3) a Prime to add to dict
q = p*p # (9) its sQuare
for c in count(9,2): # the Candidate
if c in sieve: # c's a multiple of some base prime
s = sieve.pop(c) # i.e. a composite ; or
elif c < q:
yield c # a prime
continue
else: # (c==q): # or the next base prime's square:
s=count(q+2*p,2*p) # (9+6, by 6 : 15,21,27,33,...)
p=next(ps) # (5)
q=p*p # (25)
for m in s: # the next multiple
if m not in sieve: # no duplicates
break
sieve[m] = s # original test entry: ideone.com/WFv4f
(ここの古い元のコードは、以下のTim Petersによる回答に見られるように、変更を組み込むために編集されました)。関連する議論については、これも参照してください。
同様の2-3-5-7 wheelベースのコードは ~ 2.15x 速く実行されます(これは の理論上の改善に非常に近いです3/2 * 5/4 * 7/6 = 2.1875)。
後世のために、Will Ness の美しいアルゴリズムを Python 3 用に書き直したものを次に示します。いくつかの変更が必要です (イテレータにはメソッドがなくなりまし.next()たが、新しいnext()組み込み関数があります)。他の変更は楽しみのためです (新しいものを使用すると、元のyield from <iterable>4 つのステートメントが置き換えyieldられます。さらに読みやすくするためのものです (私は 1 文字の変数名を使いすぎるのが好きではありません ;-))。
オリジナルよりも大幅に高速ですが、アルゴリズム上の理由によるものではありません。高速化は主に、元のadd()関数を削除し、代わりにインラインで行うことによるものです。
def psieve():
import itertools
yield from (2, 3, 5, 7)
D = {}
ps = psieve()
next(ps)
p = next(ps)
assert p == 3
psq = p*p
for i in itertools.count(9, 2):
if i in D: # composite
step = D.pop(i)
elif i < psq: # prime
yield i
continue
else: # composite, = p*p
assert i == psq
step = 2*p
p = next(ps)
psq = p*p
i += step
while i in D:
i += step
D[i] = step
これはもともと私のコードではありませんが、投稿する価値があります。オリジナルはここにあります: http://code.activestate.com/recipes/117119/
def gen_primes():
D = {}
q = 2 # first integer to test for primality.
while True:
if q not in D:
# not marked composite, must be prime
yield q
#first multiple of q not already marked
D[q * q] = [q]
else:
for p in D[q]:
D.setdefault(p + q, []).append(p)
# no longer need D[q], free memory
del D[q]
q += 1
これはジェネレーターなので、他のものと同じように使用してください。
primes = gen_primes()
for p in primes:
print p
デスクトップ上で 100 万個の素数を生成してセットに入れるのに 1.62 秒かかります。
セグメント化されたふるいを実行します。セグメントのサイズは、使用可能なメモリまたはビットセットの最大サイズによって決まります。
各セグメントについて、ある間隔の数値を表します[n; n + segment_size）をビットセットとして、上限の平方根より下のすべての素数でふるいにかけます。
ビットセットを使用すると、高密度の数値セットを操作するため、ハッシュテーブルやツリーデータ構造よりも少ないメモリを使用します。
そして別の答え、ここでの私の答えよりもメモリ効率が良いerat3：
import heapq
def heapprimegen():
hp= []
yield 2
yield 3
cn= 3
nn, inc= 3, 6
while 1:
while cn < nn:
yield cn
heapq.heappush(hp, (3*cn, 2*cn))
cn+= 2
cn= nn+2
nn, inc= heapq.heappushpop(hp, (nn+inc, inc))
辞書ではなく素数のヒープ（リスト）を維持します。明らかに、速度がいくらか失われます。
それを行う別の方法：
import itertools
def primeseq():
prime = [2]
num = 0
yield 2
for i in itertools.count(3, 2):
is_prime = True
for num in prime:
if i % num == 0:
is_prime = False
break
elif num ** 2 > i:
break
if is_prime:
prime.append(i)
yield i
これは非常に高速な無限ジェネレーターで、Python2 で記述されていますが、Python3 に簡単に調整できます。これを使用して 10**9 までの素数を追加するには、次を使用します。
from itertools import takewhile
from functools import partial
from operator import gt
print (sum(takewhile(partial(gt, 10**9), prime_gen_inf())))
これはセグメント化されたふるいで、Will Ness のアルゴリズムよりも高速ですが、明らかにエレガントではありません。
from operator import mul
from functools import reduce
def prod(x): return reduce(mul, x, 1)
def build_sieve(wheel):
w = prod(wheel)
w_phi = prod([p-1 for p in wheel])
rems = [a for a in range(w) if all(a % p for p in wheel)]
assert len(rems) == w_phi
inv = {a:pow(a, w_phi - 1, w) for a in rems}
try:
known_p = wheel + rems[1 : rems.index(rems[1]*rems[1])]
except ValueError:
known_p = wheel + rems[1:]
return wheel, w, w_phi, rems, inv, known_p
#Adjust the chunk variable based on your computer's architecture.
#
#Adjust the line with #! if you don't need "true" infinite. If you don't need
#primes larger than 1<<32, use array('H', []), if 1<<64 use 'L', if 1<<128 (in
#Python3) use 'Q', otherwise use empty list [].
#To save memory, comment out the lines with #*, and uncomment the commented-out
#lines
import itertools
from itertools import islice, count, compress, izip
chain_f = itertools.chain.from_iterable
from array import array
def prime_gen_inf(chunk=250000, sieve_info = build_sieve([2,3,5,7])):
""" Indefinitely yields primes """
wheel, w, w_phi, rems, inv, known_p = sieve_info
for p in known_p: yield p
new_n = 0;
while True:
size = min(chunk, (p * p - new_n) / w)
sieve = bytearray([1]) * size * w_phi
n, new_n = new_n, new_n + size * w
if not n:
zero = bytearray([0])
seen = len(known_p) - len(wheel) + 1
sieve[:seen:1] = zero * seen
p_gen = islice(prime_gen_inf(), len(wheel), None)
new_p = next(p_gen)
ps = [] #! array('H', [])
p_invs = bytearray([]) #*
while new_p * new_p < new_n:
ps.append(new_p)
p_invs.append(inv[new_p % w]) #*
new_p = next(p_gen)
for p, p_inv, modp in izip(ps, p_invs, [-n % p for p in ps]): #*
s = [(modp + p * (p_inv * (r - modp) % w)) / w for r in rems] #*
#for p in ps:
# s = [(-n%p + p * (inv[p%w] * (r - -n%p) % w)) / w for r in rems]
for i, start in enumerate(s):
slice_size = ((size - start - 1) / p + 1)
sieve[i + start * w_phi :: p * w_phi] = zero * slice_size
for p in compress(chain_f(izip(*[count(n+r, w) for r in rems])), sieve):
yield p
これは、辞書の代わりにヒープを使用した単純ですが、それほど遅くはありません。
import heapq
def heap_prime_gen_squares():
yield 2
yield 3
h = [(9, 6)]
n = 5
while True:
a, b = h[0]
while n < a:
yield n
heapq.heappush(h, (n * n, n << 1))
n += 2
heapq.heapreplace(h, (a + b, b)) # Replace h[0], which is still (a, b).
最初の 100 万個の素数のユーザー時間の測定速度 (数値が小さいほど優れています):
そのため、辞書ベースのアプローチが最速のようです。
これは複雑なヒープベースの実装であり、他のヒープベースの実装よりもはるかに高速ではありませんが (私の別の回答の速度比較を参照してください)、使用するメモリははるかに少なくなります。
この実装では、同じ数の要素を含む 2 つのヒープ (tu と wv) を使用します。各要素は int ペアです。q**2(ここでは素数)までのすべての素数を見つけるためにq、各ヒープには最大でもの2*pi(q-1)要素が含まれます。ここpi(x)で、 は より大きくない正の素数の数ですx。したがって、整数の総数は最大で4*pi(floor(sqrt(n)))です。(半分の量のものをヒープにプッシュすることで、メモリの 2 の係数を得ることができますが、それではアルゴリズムが遅くなります。)
上記の他の dict およびヒープベースのアプローチ (たとえば、erat2b、および heap_prime_gen_squares および heapprimegen) は、素数を見つけるたびにヒープまたは dict を拡張するため、約 `2*pi(n)' 整数を格納します。比較として: 1_000_000 素数を見つけるために、この実装は 4141 未満の整数を格納し、他の実装は 1_000_000 を超える整数を格納します。
import heapq
def heap_prime_gen_smallmem():
yield 2
yield 3
f = 5
fmar3 = 2
q = 7
q6 = 7 * 6
qmar3 = 4
tu = [(25, 30), (35, 30)]
vw = [(25, 30), (35, 30)]
while True:
qmar3 += 2
if qmar3 == 6:
qb = q + 4
q6b = q6 + 24
qmar3 = 2
else:
qb = q + 2
q6b = q6 + 12
if q < tu[0][0]:
d = q * q
while f < d:
a, b = vw[0]
if f < a:
yield f
else:
a, b = vw[0]
heapq.heapreplace(vw, (a + b, b))
a, b = vw[0]
while f >= a:
heapq.heapreplace(vw, (a + b, b))
a, b = vw[0]
fmar3 += 2
if fmar3 == 6:
f += 4
fmar3 = 2
else:
f += 2
c = q * qb
heapq.heappush(tu, (d, q6))
heapq.heappush(tu, (c, q6))
heapq.heappush(vw, (d, q6))
heapq.heappush(vw, (c, q6))
else:
a, b = tu[0]
heapq.heapreplace(tu, (a + b, b))
a, b = tu[0]
while q >= a:
heapq.heapreplace(tu, (a + b, b))
a, b = tu[0]
q = qb
q6 = q6b
Haskell で行われる方法に少し忠実なジェネレーターを次に示します。既知の素数の合成に対してフィルター処理を行い、残りの素数をリストに追加します。
def gen_primes():
primes = []
i = 2
while True:
prime = True
for p in primes:
if not (i % p):
prime = False
break
if prime:
yield i
primes.append(i)
i += 1