Comment implémenter un générateur infini efficace de nombres premiers en Python?
ce n'est pas un devoir, je suis juste curieux.
INFINI est le mot clé ici.
je souhaite l'utiliser comme pour P in primes(). Je crois qu'il s'agit d'une fonction intégrée à Haskell.
ainsi, la réponse ne peut pas être aussi naïve que "juste faire un tamis".
tout d'Abord, vous ne savez pas combien de nombres premiers consécutifs seront consommés. Eh bien, supposons que vous puissiez en concocter 100 à la fois. Voulez-vous utiliser l' même approche de tamis ainsi que la fréquence des nombres premiers formule?
je préfère une approche non-concurrente.
je vous Remercie pour la lecture (et l'écriture ;) )!
14 réponses
"Si j'ai vu plus loin..."
la fonction erat2 du livre de recettes peut être accélérée (d'environ 20-25%):
erat2a
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
le contrôle not (x&1) vérifie que x est impair. Toutefois, puisque les deux q et p sont impairs, en ajoutant 2*p , la moitié des pas sont évités en même temps que l'essai pour l'oddity.
erat3
si un peu de maladresse ne dérange pas, erat2 peut être accéléré de 35-40% avec les changements suivants (NB: nécessite Python 2.7+ ou Python 3+ à cause de la fonction 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
la fonction erat3 tire avantage du fait que tous les nombres premiers (sauf 2, 3, 5) modulo 30 ne donnent que huit nombres: ceux inclus dans le MODULOS frozenset. Ainsi, après rendement initial trois nombres premiers, nous commençons à partir de 7 et travaillons seulement avec les candidats.
Le filtrage des candidats utilise la fonction itertools.compress ; la fonction "magic" est dans la séquence MASK ; la fonction MASK comporte 15 éléments (il y a 15 nombres impairs par 30 nombres, comme choisi par la fonction itertools.islice ) avec un 1 pour chaque candidat possible, à partir de 7. Le cycle se répète comme spécifié par la fonction itertools.cycle .
L'introduction du filtrage des candidats nécessite une autre modification: la vérification or (x%30) not in MODULOS . L'algorithme erat2 a traité tous les nombres impairs; maintenant que l'algorithme erat3 ne traite que les candidats r30, nous devons nous assurer que tous les D.keys() ne peuvent être que de tels-faux - candidats.
Repères
résultats
sur un serveur Atom 330 Ubuntu 9.10, versions 2.6.4 et 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
sur un serveur AMD Geode LX Gentoo home, Python 2.6.5 et 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
code de référence
Un primegen.py module contient les erat2 , erat2a et erat3 des fonctions. Voici le script de test:
#!/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
puisque L'OP demande une mise en œuvre efficace , voici une amélioration significative au actif état 2002 code par David Eppstein/Alex Martelli (voir ici dans sa réponse ): ne pas enregistrer une info de prime dans le dictionnaire jusqu'à ce que sa place est vue parmi les candidats . Réduit la complexité de l'espace au-dessous de O (sqrt (n)) au lieu de O(n) , pour n nombres premiers produits ( π (N log n)) ~ 2 (n log n) / log (n log n) ~ 2 sqrt(n / log n) ). Par conséquent, la complexité du temps est également améliorée, i.e. il court plus vite .
crée un "tamis coulissant" comme un dictionnaire des multiples de courant de chaque prime de base (c.-à-d. sous le sqrt du point de production courant), avec leur étape valeurs":
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
(l'ancien code original ici a été modifié pour incorporer les changements comme vu dans la réponse par Tim Peters , ci-dessous). voir aussi ce pour une discussion.
Similaires 2-3-5-7 roue -code pistes ~ 2.15 x plus vite (ce qui est très proche de la théorie l'amélioration de la 3/2 * 5/4 * 7/6 = 2.1875 ).
pour la postérité, voici une réécriture du bel algorithme de Will Ness pour Python 3. Certains changements sont nécessaires (les itérateurs n'ont plus les méthodes .next() , mais il y a une nouvelle fonction next() ). D'autres changements sont pour le plaisir (en utilisant le nouveau yield from <iterable> remplace quatre yield énoncés dans l'original. Plus sont pour la lisibilité (Je ne suis pas un fan de la surutilisation ;-) des noms de variables d'une lettre).
c'est beaucoup plus rapide que l'original, mais pas pour des raisons algorithmiques. L'accélération est principalement due à la suppression de la fonction add() de l'original, faisant cela en ligne à la place.
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
ce n'est pas à l'origine mon code, cependant, il vaut la peine de poster. L'original peut être trouvé ici: 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
c'est un générateur, alors utilisez-le comme n'importe quel autre.
primes = gen_primes()
for p in primes:
print p
il faut 1,62 s pour générer et mettre dans un ensemble, 1 million de nombres premiers, sur mon bureau.
Faire un segmenté d'un tamis, d'où la taille d'un segment est déterminée par la mémoire disponible ou la taille maximum d'un bitset.
pour chaque segment représentent les nombres dans un certain intervalle [n; n + segment_size) comme un ensemble de bits et tamis avec tous les nombres premiers au-dessous de la racine carrée de la limite supérieure.
en utilisant un ensemble de bits utilise moins de mémoire qu'une table de hachage ou une structure de données d'arbre, parce que vous travaillez avec des ensembles denses de nombre.
c'est similaire mais environ 5% plus rapide que erat2a.
import itertools as it
def erat2b( ):
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:
# Only this part is replaced.
q += p
while q in D:
q += p
D[q] = p
une autre façon de le faire:
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
et une autre réponse, plus efficace en mémoire que mon erat3 réponse ici:
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))
il maintient un tas (une liste) de multiples premiers plutôt qu'un dictionnaire. Il perd un peu de vitesse, évidemment.
Voici une implémentation en tas compliquée, qui n'est pas beaucoup plus rapide que d'autres implémentations en tas (voir la comparaison de vitesse dans une autre de mes réponses), mais qui utilise beaucoup moins de mémoire.
cette implémentation utilise deux tas (tu et wv), qui contiennent le même nombre d'éléments. Chaque élément est une paire int. Afin de trouver tous les nombres premiers jusqu'à q**2 (où q est un nombre premier), chaque segment contient au plus 2*pi(q-1) éléments, où pi(x) est le nombre de nombres premiers positifs ne dépassant pas x . Donc le nombre total d'entiers est au plus 4*pi(floor(sqrt(n))) . (Nous pourrions obtenir un facteur sur 2 sur la mémoire en poussant la moitié des choses au tas, mais cela rendrait l'algorithme plus lent.)
d'autres approches basées sur les dicts et les tas (par exemple erat2b, et heap_prime_gen_squares et heaprimegen) au-dessus stockent environ `2*pi(n)' entiers, parce qu'ils étendent leur tas ou dict chaque fois qu'ils trouvent un premier. En tant que comparaison: pour trouver les nombres premiers 1_000_000, cette implémentation stocke moins de 4141 entiers, les autres implémentations stockent plus de 1_000_000 entiers.
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
voici un générateur qui est un peu plus fidèle à la façon dont il est fait à Haskell: filtrage contre les composites de nombres premiers connus, puis l'ajout des nombres premiers restants à la liste.
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
j'ai écrit il y a quelques temps un article sur un générateur de nombres premiers infinis:
http://stacktrace.it/2008/01/progetto-eulero-problema-3 /
c'est en italien mais vous pouvez avoir une traduction embêtante en utilisant Google: http://tinyurl.com/yzpyeom
voici un simple mais pas terriblement lent utilisant un tas au lieu d'un dict:
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).
Mes mesures de vitesse de temps de l'utilisateur pour la première de 1 million de nombres premiers (les plus petits nombres sont de mieux en mieux):
- postponed_sieve( basé sur dict): 8.553 s
- erat2b (basé sur dict): 9.513 s
- erat2a (basé sur dict): 10.313 s
- heap_prime_gen_smallmem (heap): 23.935 s
- heap_prime_gen_squares (heap): 27.302 s
- heaprimegen (dict-based): 145.029 s
les approches basées sur la dict semblent donc être les plus rapides.
voici un générateur infini assez rapide, écrit en Python2 mais facilement ajusté à Python3. De l'utiliser pour ajouter les nombres premiers jusqu'à 10**9, utilisez la commande suivante:
from itertools import takewhile
from functools import partial
from operator import gt
print (sum(takewhile(partial(gt, 10**9), prime_gen_inf())))
c'est un tamis segmenté, plus rapide mais évidemment moins élégant que L'algorithme de 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
je sais que le post est vieux, mais je suis venu par moi-même à travers cette question... Le code suivant est basé sur une idée très simple: un tamis croissant D'Ératosthènes. Cette solution est vraiment plus lente que les meilleures ici, mais elle est facile à saisir et conçue pour être lisible...
j'ai utilisé des entiers pour stocker les résultats du tamis.
En format binaire, un entier est une liste de 0 s et 1 s, 0 à la position i si i est pas un prime, 1 si c'est un prime.
L'infini requis est le résultat du fait que Python 3 entiers sont illimités.
def primes():
container, size = 1 << 2, 3 # we start with 0b100 (from right to left: 0 and 1 are not primes, 2 is
last_prime = 1
while True:
prime = next((j for j in range(last_prime+1, size) if container & 1 << j), None) # find the next prime
while not prime:
container, size = expand(container, size, 2**16) # add 65536 cells and sieve the container
prime = next((j for j in range(last_prime+1, size) if container & 1 << j), None)
yield prime
last_prime = prime
comment agrandir le conteneur? Il suffit d'ajouter un tas de 1 s à la gauche du conteneur (en format binaire) et de les tamiser. C'est
identique au tamis, avec une légère différence. Dans la norme tamis, si nous trouvons un premier i , nous commençons à traverser les cellules i*i , avec une étape de i .
ici, cela peut avoir été fait pour la première partie du conteneur. Nous avons juste besoin de commencer au début de la nouvelle partie du conteneur s'il est plus loin que
i*i .
def expand(container, size, n):
new_size = size + n
container += (1 << (new_size + 1) - 1) - (1 << size) # add n 1's
for i in range(2, new_size):
if container & (1 << i): # i is a prime
t = sum(1 << j for j in range(max(i, size // i)*i, new_size, i)) # set 1 for all mutiple
container &= ~t # cross the cells
return container, new_size
Test pour un million de nombres premiers:
import itertools
assert 78498 == len(list(itertools.takewhile(lambda p: p<1000000, primes())))