Accélérer bitstring / opérations de bits en Python?
https://stackoverflow.com/questions/2897297
-
04-10-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
J'ai écrit un générateur de nombres premiers en utilisant d'Eratosthène Sieve et Python 3.1. Le code fonctionne correctement et avec élégance à 0,32 secondes sur ideone.com pour générer des nombres premiers jusqu'à 1.000.000.
# from bitstring import BitString
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
yield 2
sub_limit = int(limit**0.5)
flags = [False, False] + [True] * (limit - 2)
# flags = BitString(limit)
# Step through all the odd numbers
for i in range(3, limit, 2):
if flags[i] is False:
# if flags[i] is True:
continue
yield i
# Exclude further multiples of the current prime number
if i <= sub_limit:
for j in range(i*3, limit, i<<1):
flags[j] = False
# flags[j] = True
Le problème est, je manque de mémoire lorsque je tente de générer des nombres jusqu'à 1000000000.
flags = [False, False] + [True] * (limit - 2)
MemoryError
Comme vous pouvez l'imaginer, allouer 1 milliard de valeurs booléennes ( 1 octet 4 ou 8 octets (voir le commentaire) chacun en Python) est vraiment pas possible, alors j'ai regardé dans bitstring . Je me suis dit, en utilisant 1 bit pour chaque indicateur serait beaucoup plus économe en mémoire. Cependant, la performance du programme a chuté de manière drastique - 24 secondes de l'exécution, pour nombre premier à 1.000.000. Ceci est probablement dû à la mise en œuvre interne de bitstring.
Vous pouvez commenter / décommenter les trois lignes pour voir ce que j'ai changé pour utiliser BitString, comme l'extrait de code ci-dessus.
Ma question est, est-il un moyen d'accélérer mon programme, avec ou sans bitstring?
Edit: S'il vous plaît tester le code vous-même avant de poster. Je ne peux pas accepter les réponses qui fonctionnent plus lentement que mon code existant naturellement.
Modifier à nouveau:
La solution
Il y a quelques petites optimisations pour votre version. En inversant les rôles de True et False, vous pouvez changer « if flags[i] is False: » à « if flags[i]: ». Et la valeur de départ pour la deuxième instruction range peut être i*i au lieu de i*3. Votre version originale prend 0,166 secondes sur mon système. Avec ces changements, la version ci-dessous prend 0,156 secondes sur mon système.
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
yield 2
sub_limit = int(limit**0.5)
flags = [True, True] + [False] * (limit - 2)
# Step through all the odd numbers
for i in range(3, limit, 2):
if flags[i]:
continue
yield i
# Exclude further multiples of the current prime number
if i <= sub_limit:
for j in range(i*i, limit, i<<1):
flags[j] = True
Cela ne permet pas votre problème de mémoire, cependant.
Déplacement dans le monde des extensions C, j'ai utilisé la version de développement de gmpy . (Disclaimer: Je suis l'un des mainteneurs.) La version de développement est appelé gmpy2 et supports entiers mutables appelés xmpz. En utilisant gmpy2 et le code suivant, j'ai une durée de 0,140 secondes. Exécution de temps pour une limite de 1.000.000.000 est de 158 secondes.
import gmpy2
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
yield 2
sub_limit = int(limit**0.5)
# Actual number is 2*bit_position + 1.
oddnums = gmpy2.xmpz(1)
current = 0
while True:
current += 1
current = oddnums.bit_scan0(current)
prime = 2 * current + 1
if prime > limit:
break
yield prime
# Exclude further multiples of the current prime number
if prime <= sub_limit:
for j in range(2*current*(current+1), limit>>1, prime):
oddnums.bit_set(j)
Repousser les optimisations, et pour autant sacrifier la clarté, je suis en cours d'exécution et temps de 0,107 secondes 123 avec le code suivant:
import gmpy2
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
yield 2
sub_limit = int(limit**0.5)
# Actual number is 2*bit_position + 1.
oddnums = gmpy2.xmpz(1)
f_set = oddnums.bit_set
f_scan0 = oddnums.bit_scan0
current = 0
while True:
current += 1
current = f_scan0(current)
prime = 2 * current + 1
if prime > limit:
break
yield prime
# Exclude further multiples of the current prime number
if prime <= sub_limit:
list(map(f_set,range(2*current*(current+1), limit>>1, prime)))
Edit: Sur la base de cet exercice, j'ai modifié gmpy2 d'accepter xmpz.bit_set(iterator). En utilisant le code suivant, le temps d'exécution pour tous les nombres premiers moins 1000000000 est de 56 secondes pour Python 2.7 et 74 secondes pour Python 3.2. (Comme il est indiqué dans les commentaires, xrange est plus rapide que range.)
import gmpy2
try:
range = xrange
except NameError:
pass
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
yield 2
sub_limit = int(limit**0.5)
oddnums = gmpy2.xmpz(1)
f_scan0 = oddnums.bit_scan0
current = 0
while True:
current += 1
current = f_scan0(current)
prime = 2 * current + 1
if prime > limit:
break
yield prime
if prime <= sub_limit:
oddnums.bit_set(iter(range(2*current*(current+1), limit>>1, prime)))
Edit # 2: Encore un essai! Je modifié gmpy2 d'accepter xmpz.bit_set(slice). En utilisant le code suivant, le temps d'exécution pour tous les nombres premiers moins 1000000000 est d'environ 40 secondes à la fois Python 2.7 et Python 3.2.
from __future__ import print_function
import time
import gmpy2
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
yield 2
sub_limit = int(limit**0.5)
flags = gmpy2.xmpz(1)
# pre-allocate the total length
flags.bit_set((limit>>1)+1)
f_scan0 = flags.bit_scan0
current = 0
while True:
current += 1
current = f_scan0(current)
prime = 2 * current + 1
if prime > limit:
break
yield prime
if prime <= sub_limit:
flags.bit_set(slice(2*current*(current+1), limit>>1, prime))
start = time.time()
result = list(prime_numbers(1000000000))
print(time.time() - start)
Edit # 3: J'ai mis à jour gmpy2 à trancher correctement support au niveau binaire d'un xmpz. Aucun changement dans la performance, mais une API beaucoup de bien. Je l'ai fait un peu de peaufinage et j'ai le temps d'arrêt à environ 37 secondes. (Voir Modifier n ° 4 à des changements dans gmpy2 2.0.0b1.)
from __future__ import print_function
import time
import gmpy2
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
sub_limit = int(limit**0.5)
flags = gmpy2.xmpz(1)
flags[(limit>>1)+1] = True
f_scan0 = flags.bit_scan0
current = 0
prime = 2
while prime <= sub_limit:
yield prime
current += 1
current = f_scan0(current)
prime = 2 * current + 1
flags[2*current*(current+1):limit>>1:prime] = True
while prime <= limit:
yield prime
current += 1
current = f_scan0(current)
prime = 2 * current + 1
start = time.time()
result = list(prime_numbers(1000000000))
print(time.time() - start)
Edit # 4: J'ai fait quelques changements dans gmpy2 2.0.0b1 qui rompent l'exemple précédent. gmpy2 ne traite plus vraie en tant que valeur spéciale qui fournit une source infinie de bits 1. -1 devrait être utilisé à la place.
from __future__ import print_function
import time
import gmpy2
def prime_numbers(limit=1000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
sub_limit = int(limit**0.5)
flags = gmpy2.xmpz(1)
flags[(limit>>1)+1] = 1
f_scan0 = flags.bit_scan0
current = 0
prime = 2
while prime <= sub_limit:
yield prime
current += 1
current = f_scan0(current)
prime = 2 * current + 1
flags[2*current*(current+1):limit>>1:prime] = -1
while prime <= limit:
yield prime
current += 1
current = f_scan0(current)
prime = 2 * current + 1
start = time.time()
result = list(prime_numbers(1000000000))
print(time.time() - start)
Edit # 5: Je l'ai fait quelques améliorations à gmpy2 2.0.0b2. Vous pouvez maintenant parcourir tous les bits qui sont actifs ou inactifs. Durée est améliorée de ~ 30%.
from __future__ import print_function
import time
import gmpy2
def sieve(limit=1000000):
'''Returns a generator that yields the prime numbers up to limit.'''
# Increment by 1 to account for the fact that slices do not include
# the last index value but we do want to include the last value for
# calculating a list of primes.
sieve_limit = gmpy2.isqrt(limit) + 1
limit += 1
# Mark bit positions 0 and 1 as not prime.
bitmap = gmpy2.xmpz(3)
# Process 2 separately. This allows us to use p+p for the step size
# when sieving the remaining primes.
bitmap[4 : limit : 2] = -1
# Sieve the remaining primes.
for p in bitmap.iter_clear(3, sieve_limit):
bitmap[p*p : limit : p+p] = -1
return bitmap.iter_clear(2, limit)
if __name__ == "__main__":
start = time.time()
result = list(sieve(1000000000))
print(time.time() - start)
print(len(result))
Autres conseils
D'accord, voici une (presque complète) je complète l'analyse comparative ont fait ce soir pour voir quel code court le plus vite. quelqu'un trouvera un peu de chance cette liste utile. J'omis tout ce qui prend plus de 30 secondes pour terminer sur ma machine.
Je voudrais remercier tout le monde qui a mis dans une entrée. J'ai gagné beaucoup de perspicacité de vos efforts, et j'espère que vous avez trop.
Ma machine: AMD ZM-86, 2,40 Ghz Dual-Core, avec 4 Go de RAM. Ceci est un ordinateur portable HP TouchSmart TX2. Notez que si je suis lié à un pastebin, Je benchmarkée tous les éléments suivants sur ma machine propre.
Je vais ajouter la référence de gmpy2 une fois que je suis en mesure de le construire.
Tous les points de référence sont testés en Python 2.6 x86
De retour une liste des nombres premiers n jusqu'à 1000000: ( Utiliser Python générateurs)
numpy version générateur de Sebastian (mise à jour) - 121 ms @
Mark Sieve + roue - 154 ms
version de Robert Slicing - 159 ms
Ma version améliorée avec découpage en tranches - 205 ms
générateur Numpy avec Énumérer - 249 ms @
Sieve de base Mark - 317 ms
amélioration de casevh sur mon original solution - 343 ms
Ma solution modifiée générateur numpy - 407 ms
Ma méthode originale dans la question - 409 ms
Python pur avec bytearray - 1394 ms @
solution BitString Scott - 6659 ms @
« @ » désigne le présent procédé est capable de générer jusqu'à n <1000000000 sur ma configuration de la machine.
En outre, si vous n'avez pas besoin générateur et que vous voulez juste la liste complète à la fois:
numpy solution de RosettaCode - 32 ms @
(La solution numpy est également capable de générer jusqu'à 1 milliard, ce qui a pris 61.6259 secondes. Je soupçonne que la mémoire a été permuté une fois, d'où le double.)
OK, donc voici ma deuxième réponse, mais que la vitesse est de l'essence, je pensais que je devais parler de la BitArray Module - même si elle est bitstring 's ennemi juré: ). Il est idéalement adapté à ce cas est non seulement une extension C (et donc plus rapide que pur Python a un espoir d'être), mais il prend également en charge les affectations de tranche. Ce n'est pas encore disponible pour Python 3 bien.
J'ai même pas essayé d'optimiser cela, je viens récrit la version bitstring. Sur ma machine, je reçois 0,16 secondes pour les nombres premiers moins d'un million.
Pour un milliard, il fonctionne parfaitement bien et termine en 2 minutes 31 secondes.
import bitarray
def prime_bitarray(limit=1000000):
yield 2
flags = bitarray.bitarray(limit)
flags.setall(False)
sub_limit = int(limit**0.5)
for i in xrange(3, limit, 2):
if not flags[i]:
yield i
if i <= sub_limit:
flags[3*i:limit:i*2] = True
question connexe: moyen le plus rapide à la liste tous les nombres premiers N-dessous en python .
Salut, je cherche aussi un code en Python pour générer des nombres premiers jusqu'à 10 ^ 9 aussi vite que je peux, ce qui est difficile en raison du problème de mémoire. jusqu'à maintenant, je suis venu avec cela pour générer des nombres premiers jusqu'à 10 ^ 6 & 10 ^ 7 (cadencer 0,171s et 1,764s respectivement sur ma vieille machine), qui semble être un peu plus rapide (au moins dans mon ordinateur) que « ma version améliorée avec découpage en tranches » de post précédent, probablement parce que j'utilise n//i-i +1 (voir ci-dessous) au lieu de len(flags[i2::i<<1]) analogue dans l'autre code. S'il vous plait corrigez moi si je me trompe. Toutes les suggestions d'amélioration sont les bienvenus.
def primes(n):
""" Returns a list of primes < n """
sieve = [True] * n
for i in xrange(3,int(n**0.5)+1,2):
if sieve[i]:
sieve[i*i::2*i]=[False]*((n-i*i-1)/(2*i)+1)
return [2] + [i for i in xrange(3,n,2) if sieve[i]]
Dans un de ses codes Xavier utilise flags[i2::i<<1] et len(flags[i2::i<<1]). J'ai calculé la taille explicitement, en utilisant le fait qu'entre i*i..n nous avons des multiples (n-i*i)//2*i de 2*i parce que nous voulons compter i*i aussi nous résumons 1 si len(sieve[i*i::2*i]) est égal à (n-i*i)//(2*i) +1. Cela rend le code plus rapide. Un générateur de base ci-dessus en fonction du code serait:
def primesgen(n):
""" Generates all primes <= n """
sieve = [True] * n
yield 2
for i in xrange(3,int(n**0.5)+1,2):
if sieve[i]:
yield i
sieve[i*i::2*i] = [False]*((n-i*i-1)/(2*i)+1)
for i in xrange(i+2,n,2):
if sieve[i]: yield i
avec un peu d'une modification peut écrire une version légèrement plus lente du code ci-dessus qui commence avec une demi-tamis de la taille sieve = [True] * (n//2) et travaille pour le même n. ne sais pas comment cela va se refléter dans la question de la mémoire. À titre d'exemple de mise en œuvre est ici la
version modifiée du code numpy rosetta (peut-être plus rapide) en commençant par un demi-tamis de la taille.
import numpy
def primesfrom3to(n):
""" Returns a array of primes, 3 <= p < n """
sieve = numpy.ones(n/2, dtype=numpy.bool)
for i in xrange(3,int(n**0.5)+1,2):
if sieve[i/2]: sieve[i*i/2::i] = False
return 2*numpy.nonzero(sieve)[0][1::]+1
Un plus rapide et plus sage générateur mémoire serait:
import numpy
def primesgen1(n):
""" Input n>=6, Generates all primes < n """
sieve1 = numpy.ones(n/6+1, dtype=numpy.bool)
sieve5 = numpy.ones(n/6 , dtype=numpy.bool)
sieve1[0] = False
yield 2; yield 3;
for i in xrange(int(n**0.5)/6+1):
if sieve1[i]:
k=6*i+1; yield k;
sieve1[ ((k*k)/6) ::k] = False
sieve5[(k*k+4*k)/6::k] = False
if sieve5[i]:
k=6*i+5; yield k;
sieve1[ ((k*k)/6) ::k] = False
sieve5[(k*k+2*k)/6::k] = False
for i in xrange(i+1,n/6):
if sieve1[i]: yield 6*i+1
if sieve5[i]: yield 6*i+5
if n%6 > 1:
if sieve1[i+1]: yield 6*(i+1)+1
ou avec un peu plus de code:
import numpy
def primesgen(n):
""" Input n>=30, Generates all primes < n """
size = n/30 + 1
sieve01 = numpy.ones(size, dtype=numpy.bool)
sieve07 = numpy.ones(size, dtype=numpy.bool)
sieve11 = numpy.ones(size, dtype=numpy.bool)
sieve13 = numpy.ones(size, dtype=numpy.bool)
sieve17 = numpy.ones(size, dtype=numpy.bool)
sieve19 = numpy.ones(size, dtype=numpy.bool)
sieve23 = numpy.ones(size, dtype=numpy.bool)
sieve29 = numpy.ones(size, dtype=numpy.bool)
sieve01[0] = False
yield 2; yield 3; yield 5;
for i in xrange(int(n**0.5)/30+1):
if sieve01[i]:
k=30*i+1; yield k;
sieve01[ (k*k)/30::k] = False
sieve07[(k*k+ 6*k)/30::k] = False
sieve11[(k*k+10*k)/30::k] = False
sieve13[(k*k+12*k)/30::k] = False
sieve17[(k*k+16*k)/30::k] = False
sieve19[(k*k+18*k)/30::k] = False
sieve23[(k*k+22*k)/30::k] = False
sieve29[(k*k+28*k)/30::k] = False
if sieve07[i]:
k=30*i+7; yield k;
sieve01[(k*k+ 6*k)/30::k] = False
sieve07[(k*k+24*k)/30::k] = False
sieve11[(k*k+16*k)/30::k] = False
sieve13[(k*k+12*k)/30::k] = False
sieve17[(k*k+ 4*k)/30::k] = False
sieve19[ (k*k)/30::k] = False
sieve23[(k*k+22*k)/30::k] = False
sieve29[(k*k+10*k)/30::k] = False
if sieve11[i]:
k=30*i+11; yield k;
sieve01[ (k*k)/30::k] = False
sieve07[(k*k+ 6*k)/30::k] = False
sieve11[(k*k+20*k)/30::k] = False
sieve13[(k*k+12*k)/30::k] = False
sieve17[(k*k+26*k)/30::k] = False
sieve19[(k*k+18*k)/30::k] = False
sieve23[(k*k+ 2*k)/30::k] = False
sieve29[(k*k+ 8*k)/30::k] = False
if sieve13[i]:
k=30*i+13; yield k;
sieve01[(k*k+24*k)/30::k] = False
sieve07[(k*k+ 6*k)/30::k] = False
sieve11[(k*k+ 4*k)/30::k] = False
sieve13[(k*k+18*k)/30::k] = False
sieve17[(k*k+16*k)/30::k] = False
sieve19[ (k*k)/30::k] = False
sieve23[(k*k+28*k)/30::k] = False
sieve29[(k*k+10*k)/30::k] = False
if sieve17[i]:
k=30*i+17; yield k;
sieve01[(k*k+ 6*k)/30::k] = False
sieve07[(k*k+24*k)/30::k] = False
sieve11[(k*k+26*k)/30::k] = False
sieve13[(k*k+12*k)/30::k] = False
sieve17[(k*k+14*k)/30::k] = False
sieve19[ (k*k)/30::k] = False
sieve23[(k*k+ 2*k)/30::k] = False
sieve29[(k*k+20*k)/30::k] = False
if sieve19[i]:
k=30*i+19; yield k;
sieve01[ (k*k)/30::k] = False
sieve07[(k*k+24*k)/30::k] = False
sieve11[(k*k+10*k)/30::k] = False
sieve13[(k*k+18*k)/30::k] = False
sieve17[(k*k+ 4*k)/30::k] = False
sieve19[(k*k+12*k)/30::k] = False
sieve23[(k*k+28*k)/30::k] = False
sieve29[(k*k+22*k)/30::k] = False
if sieve23[i]:
k=30*i+23; yield k;
sieve01[(k*k+24*k)/30::k] = False
sieve07[(k*k+ 6*k)/30::k] = False
sieve11[(k*k+14*k)/30::k] = False
sieve13[(k*k+18*k)/30::k] = False
sieve17[(k*k+26*k)/30::k] = False
sieve19[ (k*k)/30::k] = False
sieve23[(k*k+ 8*k)/30::k] = False
sieve29[(k*k+20*k)/30::k] = False
if sieve29[i]:
k=30*i+29; yield k;
sieve01[ (k*k)/30::k] = False
sieve07[(k*k+24*k)/30::k] = False
sieve11[(k*k+20*k)/30::k] = False
sieve13[(k*k+18*k)/30::k] = False
sieve17[(k*k+14*k)/30::k] = False
sieve19[(k*k+12*k)/30::k] = False
sieve23[(k*k+ 8*k)/30::k] = False
sieve29[(k*k+ 2*k)/30::k] = False
for i in xrange(i+1,n/30):
if sieve01[i]: yield 30*i+1
if sieve07[i]: yield 30*i+7
if sieve11[i]: yield 30*i+11
if sieve13[i]: yield 30*i+13
if sieve17[i]: yield 30*i+17
if sieve19[i]: yield 30*i+19
if sieve23[i]: yield 30*i+23
if sieve29[i]: yield 30*i+29
i += 1
mod30 = n%30
if mod30 > 1:
if sieve01[i]: yield 30*i+1
if mod30 > 7:
if sieve07[i]: yield 30*i+7
if mod30 > 11:
if sieve11[i]: yield 30*i+11
if mod30 > 13:
if sieve13[i]: yield 30*i+13
if mod30 > 17:
if sieve17[i]: yield 30*i+17
if mod30 > 19:
if sieve19[i]: yield 30*i+19
if mod30 > 23:
if sieve23[i]: yield 30*i+23
Ps. Si vous avez des suggestions sur la façon d'accélérer ce code, tout de changer l'ordre des opérations à quoi que ce soit pré-calcul, s'il vous plaît commentaire
Voici une version que j'ai écrit un dos while; il pourrait être intéressant de comparer avec le vôtre pour la vitesse. Il ne fait rien au sujet des problèmes d'espace, bien que (en fait, ils sont probablement pires que avec votre version).
from math import sqrt
def basicSieve(n):
"""Given a positive integer n, generate the primes < n."""
s = [1]*n
for p in xrange(2, 1+int(sqrt(n-1))):
if s[p]:
a = p*p
s[a::p] = [0] * -((a-n)//p)
for p in xrange(2, n):
if s[p]:
yield p
Je versions plus rapides, en utilisant une roue, mais ils sont beaucoup plus compliquées. Source originale est .
D'accord, voici la version à l'aide d'une roue. wheelSieve est le point d'entrée principal.
from math import sqrt
from bisect import bisect_left
def basicSieve(n):
"""Given a positive integer n, generate the primes < n."""
s = [1]*n
for p in xrange(2, 1+int(sqrt(n-1))):
if s[p]:
a = p*p
s[a::p] = [0] * -((a-n)//p)
for p in xrange(2, n):
if s[p]:
yield p
class Wheel(object):
"""Class representing a wheel.
Attributes:
primelimit -> wheel covers primes < primelimit.
For example, given a primelimit of 6
the wheel primes are 2, 3, and 5.
primes -> list of primes less than primelimit
size -> product of the primes in primes; the modulus of the wheel
units -> list of units modulo size
phi -> number of units
"""
def __init__(self, primelimit):
self.primelimit = primelimit
self.primes = list(basicSieve(primelimit))
# compute the size of the wheel
size = 1
for p in self.primes:
size *= p
self.size = size
# find the units by sieving
units = [1] * self.size
for p in self.primes:
units[::p] = [0]*(self.size//p)
self.units = [i for i in xrange(self.size) if units[i]]
# number of units
self.phi = len(self.units)
def to_index(self, n):
"""Compute alpha(n), where alpha is an order preserving map
from the set of units modulo self.size to the nonnegative integers.
If n is not a unit, the index of the first unit greater than n
is given."""
return bisect_left(self.units, n % self.size) + n // self.size * self.phi
def from_index(self, i):
"""Inverse of to_index."""
return self.units[i % self.phi] + i // self.phi * self.size
def wheelSieveInner(n, wheel):
"""Given a positive integer n and a wheel, find the wheel indices of
all primes that are less than n, and that are units modulo the
wheel modulus.
"""
# renaming to avoid the overhead of attribute lookups
U = wheel.units
wS = wheel.size
# inverse of unit map
UI = dict((u, i) for i, u in enumerate(U))
nU = len(wheel.units)
sqroot = 1+int(sqrt(n-1)) # ceiling of square root of n
# corresponding index (index of next unit up)
sqrti = bisect_left(U, sqroot % wS) + sqroot//wS*nU
upper = bisect_left(U, n % wS) + n//wS*nU
ind2 = bisect_left(U, 2 % wS) + 2//wS*nU
s = [1]*upper
for i in xrange(ind2, sqrti):
if s[i]:
q = i//nU
u = U[i%nU]
p = q*wS+u
u2 = u*u
aq, au = (p+u)*q+u2//wS, u2%wS
wp = p * nU
for v in U:
# eliminate entries corresponding to integers
# congruent to p*v modulo p*wS
uvr = u*v%wS
m = aq + (au > uvr)
bot = (m + (q*v + u*v//wS - m) % p) * nU + UI[uvr]
s[bot::wp] = [0]*-((bot-upper)//wp)
return s
def wheelSieve(n, wheel=Wheel(10)):
"""Given a positive integer n, generate the list of primes <= n."""
n += 1
wS = wheel.size
U = wheel.units
nU = len(wheel.units)
s = wheelSieveInner(n, wheel)
return (wheel.primes[:bisect_left(wheel.primes, n)] +
[p//nU*wS + U[p%nU] for p in xrange(bisect_left(U, 2 % wS)
+ 2//wS*nU, len(s)) if s[p]])
Quant à ce qu'une roue est: bien, vous savez que (à l'exception de 2), tous les nombres premiers sont impairs, donc la plupart des tamis manquent tous les nombres pairs. De même, vous pouvez aller un peu plus loin et que tous les nombres premiers avis (sauf 2 et 3) sont conformes à 1 ou 5 modulo 6 (== 2 * 3), de sorte que vous pouvez vous en sortir avec seulement stocker les entrées pour ces numéros dans votre tamis . La prochaine étape est de constater que tous les nombres premiers (à l'exception de 2, 3 et 5) sont congruentes l'une de 1, 7, 11, 13, 17, 19, 23, 29 (modulo 30) (ici 30 == 2 * 3 * 5), et ainsi de suite.
Une amélioration de la vitesse que vous pouvez faire en utilisant bitstring est d'utiliser la méthode du « jeu » lorsque vous excluez des multiples du nombre actuel.
Ainsi, la section devient vitale
for i in range(3, limit, 2):
if flags[i]:
yield i
if i <= sub_limit:
flags.set(1, range(i*3, limit, i*2))
Sur ma machine, cela va environ 3 fois plus rapide que l'original.
Mon autre pensée était d'utiliser la chaîne binaire pour représenter uniquement les nombres impairs. Vous pouvez ensuite trouver les bits hors service en utilisant flags.findall([0]) qui retourne un générateur. Je ne sais pas si ce serait beaucoup plus rapide (trouver des modèles de bits n'est pas très facile à faire efficacement).
[La divulgation complète: J'ai écrit le module bitstring, donc j'ai une certaine fierté en jeu ici!]
En comparaison, j'ai aussi pris le courage de la méthode bitstring afin qu'il soit le faire de la même manière, mais sans le contrôle d'exécution, etc. Je pense que cela donne une limite inférieure raisonnable pour Python pur qui travaille pour un milliard éléments (sans changer l'algorithme, qui je pense est la triche!)
def prime_pure(limit=1000000):
yield 2
flags = bytearray((limit + 7) // 8)
sub_limit = int(limit**0.5)
for i in xrange(3, limit, 2):
byte, bit = divmod(i, 8)
if not flags[byte] & (128 >> bit):
yield i
if i <= sub_limit:
for j in xrange(i*3, limit, i*2):
byte, bit = divmod(j, 8)
flags[byte] |= (128 >> bit)
Sur ma machine, ce fonctionne en environ 0,62 secondes pour un million d'éléments, ce qui signifie qu'il est d'environ un quart de la vitesse de la réponse BitArray. Je pense que est tout à fait raisonnable pour Python pur.
Les types entiers de Python peut être de taille arbitraire, de sorte que vous ne devriez pas avoir besoin d'une bibliothèque bitstring intelligente pour représenter un ensemble de bits, juste un seul numéro.
Voici le code. Il gère une limite de 1.000.000 avec aisance, et peut même gérer 10.000.000 sans trop se plaindre:
def multiples_of(n, step, limit):
bits = 1 << n
old_bits = bits
max = 1 << limit
while old_bits < max:
old_bits = bits
bits += bits << step
step *= 2
return old_bits
def prime_numbers(limit=10000000):
'''Prime number generator. Yields the series
2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ...
using Sieve of Eratosthenes.
'''
yield 2
sub_limit = int(limit**0.5)
flags = ((1 << (limit - 2)) - 1) << 2
# Step through all the odd numbers
for i in xrange(3, limit, 2):
if not (flags & (1 << i)):
continue
yield i
# Exclude further multiples of the current prime number
if i <= sub_limit:
flags &= ~multiples_of(i * 3, i << 1, limit)
La fonction multiples_of évite le coût de la mise en chaque multiple individuellement. Il définit le bit initial, il se déplace par l'étape initiale (i << 1) et ajoute le résultat à lui-même. Il double alors l'étape, de sorte que les copies de quart de travail suivant les deux bits, et ainsi de suite jusqu'à ce qu'il atteigne la limite. Ceci définit tous les multiples d'un nombre jusqu'à la limite dans O (log (limite)) temps.
L'une des options que vous pouvez regarder est simplement un module compilation C / C ++ de sorte que vous avez un accès direct aux fonctionnalités tripotant bits. Vous pouvez écrire autant que je sache quelque chose de cette nature et ne renvoie les résultats à la fin des calculs effectués en C / C ++. Maintenant que je tape cela, vous pouvez également regarder numpy qui ne stockent des tableaux en C pour la vitesse native. Je ne sais pas si ce sera plus vite que le module bitstring, bien!
Une autre option vraiment stupide, mais qui peut être utile si vous avez vraiment besoin d'une grande liste de nombres premiers nombre disponible très rapidement. Dites, si vous avez besoin d'outil pour les problèmes de projet de réponse d'Euler (si le problème est tout simplement votre code comme l'optimisation d'un jeu de l'esprit, il est hors de propos).
Utiliser une solution lente à générer la liste et l'enregistrer dans un fichier source python, dit primes.py qui contiennent seulement:
primes = [ a list of a million primes numbers here ]
Maintenant, pour les utiliser il vous suffit de faire import primes et vous les obtenez avec une empreinte mémoire minimale à la vitesse du disque IO. La complexité est le nombre de nombres premiers: -)
Même si vous avez utilisé une solution mal optimisé pour générer cette liste, il ne sera fait une fois et il n'a pas beaucoup d'importance.
Vous pourriez probablement le rendre encore plus rapide en utilisant cornichon / unpickle, mais vous voyez l'idée ...
Vous pouvez utiliser un segmentée d'Eratosthène Sieve. La mémoire utilisée pour chaque segment est réutilisé pour le segment suivant.
Voici un code python3 qui utilise moins de mémoire que les BitArray / solutions chaîne binaire affichée à ce jour et environ 1/8 la mémoire de primesgen de Robert William Hanks (), lors de l'exécution plus rapide que primesgen () (légèrement plus rapide à 1.000.000, en utilisant 37KB de mémoire, 3x plus rapide que primesgen () à l'aide 1.000.000.000 sous 34MB). L'augmentation de la taille du morceau (chunk variable dans le code) utilise plus de mémoire, mais accélère le programme, jusqu'à une limite - je pris une valeur de telle sorte que sa contribution à la mémoire est sous environ 10% de son tamis n // 30 octets. Il est pas aussi efficace que mémoire Will générateur infini de Ness (voir aussi https://stackoverflow.com/a/19391111/5439078 ) car il enregistre (et retourne à la fin, sous forme de comprimé tous les nombres premiers) calculées.
Ce correctement devrait fonctionner aussi longtemps que le calcul de la racine carrée sort précise (environ 2 ** 51 si Python utilise double 64 bits). Cependant, vous ne devriez pas utiliser ce programme pour trouver des nombres premiers que les grands!
(je ne l'ai pas recalculent les compensations, juste les a pris de code de Robert William Hanks. Merci Robert!)
import numpy as np
def prime_30_gen(n):
""" Input n, int-like. Yields all primes < n as Python ints"""
wheel = np.array([2,3,5])
yield from wheel[wheel < n].tolist()
powers = 1 << np.arange(8, dtype='u1')
odds = np.array([1, 7, 11, 13, 17, 19, 23, 29], dtype='i8')
offsets=np.array([[0,6,10,12,16,18,22,28],[6,24,16,12,4,0,22,10],
[0,6,20,12,26,18,2,8], [24,6,4,18,16,0,28,10],
[6,24,26,12,14,0,2,20], [0,24,10,18,4,12,28,22],
[24,6,14,18,26,0,8,20], [0,24,20,18,14,12,8,2]],
dtype='i8')
offsets = {d:f for d,f in zip(odds, offsets)}
sieve = 255 * np.ones((n + 28) // 30, dtype='u1')
if n % 30 > 1: sieve[-1] >>= 8 - odds.searchsorted(n % 30)
sieve[0] &= ~1
for i in range((int((n - 1)**0.5) + 29) // 30):
for odd in odds[(sieve[i] & powers).nonzero()[0]]:
k = i * 30 + odd
yield int(k)
for clear_bit, off in zip(~powers, offsets[odd]):
sieve[(k * (k + off)) // 30 :: k] &= clear_bit
chunk = min(1 + (n >> 13), 1<<15)
for i in range(i + 1, len(sieve), chunk):
a = (sieve[i : i + chunk, np.newaxis] & powers)
a = np.flatnonzero(a.astype('bool')) + i * 8
a = (a // 8 * 30 + odds[a & 7]).tolist()
yield from a
return sieve
Note Side: Si vous voulez vitesse réelle et ont la 2 Go de RAM nécessaire pour la liste des nombres premiers à 10 ** 9, vous devez utiliser pyprimesieve (sur https://pypi.python.org/ , en utilisant primesieve http: / /primesieve.org/ ).
