Le moyen le plus rapide de lister tous les nombres premiers inférieurs à N
https://stackoverflow.com/questions/2068372
-
20-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
C'est le meilleur algorithme que j'ai pu trouver.
def get_primes(n):
numbers = set(range(n, 1, -1))
primes = []
while numbers:
p = numbers.pop()
primes.append(p)
numbers.difference_update(set(range(p*2, n+1, p)))
return primes
>>> timeit.Timer(stmt='get_primes.get_primes(1000000)', setup='import get_primes').timeit(1)
1.1499958793645562
Peut-on le rendre encore plus rapide ?
Ce code a un défaut :Depuis numbers est un ensemble non ordonné, il n'y a aucune garantie que numbers.pop() supprimera le nombre le plus bas de l’ensemble.Néanmoins, cela fonctionne (du moins pour moi) pour certains nombres d'entrée :
>>> sum(get_primes(2000000))
142913828922L
#That's the correct sum of all numbers below 2 million
>>> 529 in get_primes(1000)
False
>>> 529 in get_primes(530)
True
La solution
Attention: Résultats timeit peuvent varier en raison des différences dans le matériel ou
version de Python.
Voici un script qui compare un certain nombre de mises en œuvre:
- ambi_sieve_plain,
- rwh_primes ,
- rwh_primes1 ,
- rwh_primes2 ,
- sieveOfAtkin ,
- Crible d'Ératosthène,
- sundaram3 ,
- sieve_wheel_30 ,
- ambi_sieve (nécessite numpy )
- primesfrom3to (nécessite numpy)
- primesfrom2to (nécessite numpy)
Un grand merci à stephan pour amener sieve_wheel_30 à mon attention. Le crédit va Robert William Hanks pour primesfrom2to, primesfrom3to, rwh_primes, rwh_primes1 et rwh_primes2.
Des méthodes Python simples testées, avec psyco , pour n = 1000000, rwh_primes1 a été le plus rapide testé.
+---------------------+-------+
| Method | ms |
+---------------------+-------+
| rwh_primes1 | 43.0 |
| sieveOfAtkin | 46.4 |
| rwh_primes | 57.4 |
| sieve_wheel_30 | 63.0 |
| rwh_primes2 | 67.8 |
| sieveOfEratosthenes | 147.0 |
| ambi_sieve_plain | 152.0 |
| sundaram3 | 194.0 |
+---------------------+-------+
Des méthodes Python simples testées, sans psyco , pour n = 1000000, rwh_primes2 a été le plus rapide.
+---------------------+-------+
| Method | ms |
+---------------------+-------+
| rwh_primes2 | 68.1 |
| rwh_primes1 | 93.7 |
| rwh_primes | 94.6 |
| sieve_wheel_30 | 97.4 |
| sieveOfEratosthenes | 178.0 |
| ambi_sieve_plain | 286.0 |
| sieveOfAtkin | 314.0 |
| sundaram3 | 416.0 |
+---------------------+-------+
De toutes les méthodes testées, permettant numpy , pour n = 1000000, primesfrom2to a été le plus rapide testé.
+---------------------+-------+
| Method | ms |
+---------------------+-------+
| primesfrom2to | 15.9 |
| primesfrom3to | 18.4 |
| ambi_sieve | 29.3 |
+---------------------+-------+
synchronisations ont été mesurées à l'aide de la commande:
python -mtimeit -s"import primes" "primes.{method}(1000000)"
avec {method} remplacé par chacun des noms de méthodes.
primes.py:
#!/usr/bin/env python
import psyco; psyco.full()
from math import sqrt, ceil
import numpy as np
def rwh_primes(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" 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]]
def rwh_primes1(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Returns a list of primes < n """
sieve = [True] * (n/2)
for i in xrange(3,int(n**0.5)+1,2):
if sieve[i/2]:
sieve[i*i/2::i] = [False] * ((n-i*i-1)/(2*i)+1)
return [2] + [2*i+1 for i in xrange(1,n/2) if sieve[i]]
def rwh_primes2(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Input n>=6, Returns a list of primes, 2 <= p < n """
correction = (n%6>1)
n = {0:n,1:n-1,2:n+4,3:n+3,4:n+2,5:n+1}[n%6]
sieve = [True] * (n/3)
sieve[0] = False
for i in xrange(int(n**0.5)/3+1):
if sieve[i]:
k=3*i+1|1
sieve[ ((k*k)/3) ::2*k]=[False]*((n/6-(k*k)/6-1)/k+1)
sieve[(k*k+4*k-2*k*(i&1))/3::2*k]=[False]*((n/6-(k*k+4*k-2*k*(i&1))/6-1)/k+1)
return [2,3] + [3*i+1|1 for i in xrange(1,n/3-correction) if sieve[i]]
def sieve_wheel_30(N):
# http://zerovolt.com/?p=88
''' Returns a list of primes <= N using wheel criterion 2*3*5 = 30
Copyright 2009 by zerovolt.com
This code is free for non-commercial purposes, in which case you can just leave this comment as a credit for my work.
If you need this code for commercial purposes, please contact me by sending an email to: info [at] zerovolt [dot] com.'''
__smallp = ( 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59,
61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139,
149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227,
229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311,
313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401,
409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491,
499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599,
601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683,
691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797,
809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887,
907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997)
wheel = (2, 3, 5)
const = 30
if N < 2:
return []
if N <= const:
pos = 0
while __smallp[pos] <= N:
pos += 1
return list(__smallp[:pos])
# make the offsets list
offsets = (7, 11, 13, 17, 19, 23, 29, 1)
# prepare the list
p = [2, 3, 5]
dim = 2 + N // const
tk1 = [True] * dim
tk7 = [True] * dim
tk11 = [True] * dim
tk13 = [True] * dim
tk17 = [True] * dim
tk19 = [True] * dim
tk23 = [True] * dim
tk29 = [True] * dim
tk1[0] = False
# help dictionary d
# d[a , b] = c ==> if I want to find the smallest useful multiple of (30*pos)+a
# on tkc, then I need the index given by the product of [(30*pos)+a][(30*pos)+b]
# in general. If b < a, I need [(30*pos)+a][(30*(pos+1))+b]
d = {}
for x in offsets:
for y in offsets:
res = (x*y) % const
if res in offsets:
d[(x, res)] = y
# another help dictionary: gives tkx calling tmptk[x]
tmptk = {1:tk1, 7:tk7, 11:tk11, 13:tk13, 17:tk17, 19:tk19, 23:tk23, 29:tk29}
pos, prime, lastadded, stop = 0, 0, 0, int(ceil(sqrt(N)))
# inner functions definition
def del_mult(tk, start, step):
for k in xrange(start, len(tk), step):
tk[k] = False
# end of inner functions definition
cpos = const * pos
while prime < stop:
# 30k + 7
if tk7[pos]:
prime = cpos + 7
p.append(prime)
lastadded = 7
for off in offsets:
tmp = d[(7, off)]
start = (pos + prime) if off == 7 else (prime * (const * (pos + 1 if tmp < 7 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 11
if tk11[pos]:
prime = cpos + 11
p.append(prime)
lastadded = 11
for off in offsets:
tmp = d[(11, off)]
start = (pos + prime) if off == 11 else (prime * (const * (pos + 1 if tmp < 11 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 13
if tk13[pos]:
prime = cpos + 13
p.append(prime)
lastadded = 13
for off in offsets:
tmp = d[(13, off)]
start = (pos + prime) if off == 13 else (prime * (const * (pos + 1 if tmp < 13 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 17
if tk17[pos]:
prime = cpos + 17
p.append(prime)
lastadded = 17
for off in offsets:
tmp = d[(17, off)]
start = (pos + prime) if off == 17 else (prime * (const * (pos + 1 if tmp < 17 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 19
if tk19[pos]:
prime = cpos + 19
p.append(prime)
lastadded = 19
for off in offsets:
tmp = d[(19, off)]
start = (pos + prime) if off == 19 else (prime * (const * (pos + 1 if tmp < 19 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 23
if tk23[pos]:
prime = cpos + 23
p.append(prime)
lastadded = 23
for off in offsets:
tmp = d[(23, off)]
start = (pos + prime) if off == 23 else (prime * (const * (pos + 1 if tmp < 23 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# 30k + 29
if tk29[pos]:
prime = cpos + 29
p.append(prime)
lastadded = 29
for off in offsets:
tmp = d[(29, off)]
start = (pos + prime) if off == 29 else (prime * (const * (pos + 1 if tmp < 29 else 0) + tmp) )//const
del_mult(tmptk[off], start, prime)
# now we go back to top tk1, so we need to increase pos by 1
pos += 1
cpos = const * pos
# 30k + 1
if tk1[pos]:
prime = cpos + 1
p.append(prime)
lastadded = 1
for off in offsets:
tmp = d[(1, off)]
start = (pos + prime) if off == 1 else (prime * (const * pos + tmp) )//const
del_mult(tmptk[off], start, prime)
# time to add remaining primes
# if lastadded == 1, remove last element and start adding them from tk1
# this way we don't need an "if" within the last while
if lastadded == 1:
p.pop()
# now complete for every other possible prime
while pos < len(tk1):
cpos = const * pos
if tk1[pos]: p.append(cpos + 1)
if tk7[pos]: p.append(cpos + 7)
if tk11[pos]: p.append(cpos + 11)
if tk13[pos]: p.append(cpos + 13)
if tk17[pos]: p.append(cpos + 17)
if tk19[pos]: p.append(cpos + 19)
if tk23[pos]: p.append(cpos + 23)
if tk29[pos]: p.append(cpos + 29)
pos += 1
# remove exceeding if present
pos = len(p) - 1
while p[pos] > N:
pos -= 1
if pos < len(p) - 1:
del p[pos+1:]
# return p list
return p
def sieveOfEratosthenes(n):
"""sieveOfEratosthenes(n): return the list of the primes < n."""
# Code from: <dickinsm@gmail.com>, Nov 30 2006
# http://groups.google.com/group/comp.lang.python/msg/f1f10ced88c68c2d
if n <= 2:
return []
sieve = range(3, n, 2)
top = len(sieve)
for si in sieve:
if si:
bottom = (si*si - 3) // 2
if bottom >= top:
break
sieve[bottom::si] = [0] * -((bottom - top) // si)
return [2] + [el for el in sieve if el]
def sieveOfAtkin(end):
"""sieveOfAtkin(end): return a list of all the prime numbers <end
using the Sieve of Atkin."""
# Code by Steve Krenzel, <Sgk284@gmail.com>, improved
# Code: https://web.archive.org/web/20080324064651/http://krenzel.info/?p=83
# Info: http://en.wikipedia.org/wiki/Sieve_of_Atkin
assert end > 0
lng = ((end-1) // 2)
sieve = [False] * (lng + 1)
x_max, x2, xd = int(sqrt((end-1)/4.0)), 0, 4
for xd in xrange(4, 8*x_max + 2, 8):
x2 += xd
y_max = int(sqrt(end-x2))
n, n_diff = x2 + y_max*y_max, (y_max << 1) - 1
if not (n & 1):
n -= n_diff
n_diff -= 2
for d in xrange((n_diff - 1) << 1, -1, -8):
m = n % 12
if m == 1 or m == 5:
m = n >> 1
sieve[m] = not sieve[m]
n -= d
x_max, x2, xd = int(sqrt((end-1) / 3.0)), 0, 3
for xd in xrange(3, 6 * x_max + 2, 6):
x2 += xd
y_max = int(sqrt(end-x2))
n, n_diff = x2 + y_max*y_max, (y_max << 1) - 1
if not(n & 1):
n -= n_diff
n_diff -= 2
for d in xrange((n_diff - 1) << 1, -1, -8):
if n % 12 == 7:
m = n >> 1
sieve[m] = not sieve[m]
n -= d
x_max, y_min, x2, xd = int((2 + sqrt(4-8*(1-end)))/4), -1, 0, 3
for x in xrange(1, x_max + 1):
x2 += xd
xd += 6
if x2 >= end: y_min = (((int(ceil(sqrt(x2 - end))) - 1) << 1) - 2) << 1
n, n_diff = ((x*x + x) << 1) - 1, (((x-1) << 1) - 2) << 1
for d in xrange(n_diff, y_min, -8):
if n % 12 == 11:
m = n >> 1
sieve[m] = not sieve[m]
n += d
primes = [2, 3]
if end <= 3:
return primes[:max(0,end-2)]
for n in xrange(5 >> 1, (int(sqrt(end))+1) >> 1):
if sieve[n]:
primes.append((n << 1) + 1)
aux = (n << 1) + 1
aux *= aux
for k in xrange(aux, end, 2 * aux):
sieve[k >> 1] = False
s = int(sqrt(end)) + 1
if s % 2 == 0:
s += 1
primes.extend([i for i in xrange(s, end, 2) if sieve[i >> 1]])
return primes
def ambi_sieve_plain(n):
s = range(3, n, 2)
for m in xrange(3, int(n**0.5)+1, 2):
if s[(m-3)/2]:
for t in xrange((m*m-3)/2,(n>>1)-1,m):
s[t]=0
return [2]+[t for t in s if t>0]
def sundaram3(max_n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/2073279#2073279
numbers = range(3, max_n+1, 2)
half = (max_n)//2
initial = 4
for step in xrange(3, max_n+1, 2):
for i in xrange(initial, half, step):
numbers[i-1] = 0
initial += 2*(step+1)
if initial > half:
return [2] + filter(None, numbers)
################################################################################
# Using Numpy:
def ambi_sieve(n):
# http://tommih.blogspot.com/2009/04/fast-prime-number-generator.html
s = np.arange(3, n, 2)
for m in xrange(3, int(n ** 0.5)+1, 2):
if s[(m-3)/2]:
s[(m*m-3)/2::m]=0
return np.r_[2, s[s>0]]
def primesfrom3to(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Returns a array of primes, p < n """
assert n>=2
sieve = np.ones(n/2, dtype=np.bool)
for i in xrange(3,int(n**0.5)+1,2):
if sieve[i/2]:
sieve[i*i/2::i] = False
return np.r_[2, 2*np.nonzero(sieve)[0][1::]+1]
def primesfrom2to(n):
# https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
""" Input n>=6, Returns a array of primes, 2 <= p < n """
sieve = np.ones(n/3 + (n%6==2), dtype=np.bool)
sieve[0] = False
for i in xrange(int(n**0.5)/3+1):
if sieve[i]:
k=3*i+1|1
sieve[ ((k*k)/3) ::2*k] = False
sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False
return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]
if __name__=='__main__':
import itertools
import sys
def test(f1,f2,num):
print('Testing {f1} and {f2} return same results'.format(
f1=f1.func_name,
f2=f2.func_name))
if not all([a==b for a,b in itertools.izip_longest(f1(num),f2(num))]):
sys.exit("Error: %s(%s) != %s(%s)"%(f1.func_name,num,f2.func_name,num))
n=1000000
test(sieveOfAtkin,sieveOfEratosthenes,n)
test(sieveOfAtkin,ambi_sieve,n)
test(sieveOfAtkin,ambi_sieve_plain,n)
test(sieveOfAtkin,sundaram3,n)
test(sieveOfAtkin,sieve_wheel_30,n)
test(sieveOfAtkin,primesfrom3to,n)
test(sieveOfAtkin,primesfrom2to,n)
test(sieveOfAtkin,rwh_primes,n)
test(sieveOfAtkin,rwh_primes1,n)
test(sieveOfAtkin,rwh_primes2,n)
Exécuter les tests de script que toutes les implémentations donnent le même résultat.
Autres conseils
Le plus rapide et plus pur code Python-sage mémoire:
def primes(n):
""" Returns a list of primes < n """
sieve = [True] * n
for i in range(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 range(3,n,2) if sieve[i]]
ou à partir de la moitié tamis
def primes1(n):
""" Returns a list of primes < n """
sieve = [True] * (n//2)
for i in range(3,int(n**0.5)+1,2):
if sieve[i//2]:
sieve[i*i//2::i] = [False] * ((n-i*i-1)//(2*i)+1)
return [2] + [2*i+1 for i in range(1,n//2) if sieve[i]]
Plus rapide et plus code numpy-sage mémoire:
import numpy
def primesfrom3to(n):
""" Returns a array of primes, 3 <= p < n """
sieve = numpy.ones(n//2, dtype=numpy.bool)
for i in range(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
une variation plus rapide en commençant par un tiers d'un tamis:
import numpy
def primesfrom2to(n):
""" Input n>=6, Returns a array of primes, 2 <= p < n """
sieve = numpy.ones(n//3 + (n%6==2), dtype=numpy.bool)
for i in range(1,int(n**0.5)//3+1):
if sieve[i]:
k=3*i+1|1
sieve[ k*k//3 ::2*k] = False
sieve[k*(k-2*(i&1)+4)//3::2*k] = False
return numpy.r_[2,3,((3*numpy.nonzero(sieve)[0][1:]+1)|1)]
A (dur-à-code) la version pure python du code ci-dessus serait:
def primes2(n):
""" Input n>=6, Returns a list of primes, 2 <= p < n """
n, correction = n-n%6+6, 2-(n%6>1)
sieve = [True] * (n//3)
for i in range(1,int(n**0.5)//3+1):
if sieve[i]:
k=3*i+1|1
sieve[ k*k//3 ::2*k] = [False] * ((n//6-k*k//6-1)//k+1)
sieve[k*(k-2*(i&1)+4)//3::2*k] = [False] * ((n//6-k*(k-2*(i&1)+4)//6-1)//k+1)
return [2,3] + [3*i+1|1 for i in range(1,n//3-correction) if sieve[i]]
Malheureusement, pur-python ne pas adopter la manière plus simple et plus rapide de faire numpy cession, et d'appeler len() dans la boucle comme dans [False]*len(sieve[((k*k)//3)::2*k]) est trop lent. J'ai donc dû improviser pour corriger l'entrée (et éviter plus de maths) et faire un peu extrême (et douloureux) mathématiques magie.
Personnellement, je pense qu'il est dommage que numpy (qui est si largement utilisé) ne fait pas partie de la bibliothèque standard de Python, et que les améliorations de la syntaxe et de la vitesse semblent être complètement négligé par les développeurs Python.
Il y a un échantillon assez propre du livre de recettes Python ici - la version la plus rapide proposée sur cette URL est:
import itertools
def erat2( ):
D = { }
yield 2
for q in itertools.islice(itertools.count(3), 0, None, 2):
p = D.pop(q, None)
if p is None:
D[q*q] = q
yield q
else:
x = p + q
while x in D or not (x&1):
x += p
D[x] = p
qui donnerait
def get_primes_erat(n):
return list(itertools.takewhile(lambda p: p<n, erat2()))
La mesure à l'invite du shell (comme je préfère le faire) avec ce code dans pri.py, je constate:
$ python2.5 -mtimeit -s'import pri' 'pri.get_primes(1000000)'
10 loops, best of 3: 1.69 sec per loop
$ python2.5 -mtimeit -s'import pri' 'pri.get_primes_erat(1000000)'
10 loops, best of 3: 673 msec per loop
il semble que la solution livre de recettes est deux fois plus rapide.
Utilisation Sieve de Sundaram, je pense que je battu le record de pur-Python :
def sundaram3(max_n):
numbers = range(3, max_n+1, 2)
half = (max_n)//2
initial = 4
for step in xrange(3, max_n+1, 2):
for i in xrange(initial, half, step):
numbers[i-1] = 0
initial += 2*(step+1)
if initial > half:
return [2] + filter(None, numbers)
comparasion:
C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.get_primes_erat(1000000)"
10 loops, best of 3: 710 msec per loop
C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.daniel_sieve_2(1000000)"
10 loops, best of 3: 435 msec per loop
C:\USERS>python -m timeit -n10 -s "import get_primes" "get_primes.sundaram3(1000000)"
10 loops, best of 3: 327 msec per loop
L'algorithme est rapide, mais il a un grave défaut:
>>> sorted(get_primes(530))
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73,
79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163,
167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251,
257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443,
449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 527, 529]
>>> 17*31
527
>>> 23*23
529
Vous supposez que numbers.pop() retournerait le plus petit nombre dans l'ensemble, mais ce n'est pas garanti du tout. Les ensembles sont non-ordonnées et pop() Supprime et renvoie un arbitraire élément, donc il ne peut pas être utilisé pour sélectionner le prochain premier des numéros restants.
vraiment solution la plus rapide avec N suffisamment grand serait de télécharger un liste précalculée de nombres premiers , stocker comme un tuple et faire quelque chose comme:
for pos,i in enumerate(primes):
if i > N:
print primes[:pos]
Si N > primes[-1] uniquement puis calculer plus les nombres premiers et enregistrer la nouvelle liste dans votre code, donc la prochaine fois, il est tout aussi rapide.
Toujours penser en dehors des sentiers battus.
Si vous ne voulez pas de réinventer la roue, vous pouvez installer la bibliothèque mathématique symbolique sympy (oui il est Python 3 compatible)
pip install sympy
Et utiliser le primerange fonction
from sympy import sieve
primes = list(sieve.primerange(1, 10**6))
Si vous acceptez itertools mais pas numpy, voici une adaptation de rwh_primes2 pour Python 3 qui tourne deux fois plus vite sur ma machine. Le seul changement substantiel utilise un bytearray au lieu d'une liste pour le booléen, et en utilisant compriment au lieu d'une compréhension de la liste pour construire la liste finale. (J'ajouter cela comme un commentaire comme moarningsun si je pouvais.)
import itertools
izip = itertools.zip_longest
chain = itertools.chain.from_iterable
compress = itertools.compress
def rwh_primes2_python3(n):
""" Input n>=6, Returns a list of primes, 2 <= p < n """
zero = bytearray([False])
size = n//3 + (n % 6 == 2)
sieve = bytearray([True]) * size
sieve[0] = False
for i in range(int(n**0.5)//3+1):
if sieve[i]:
k=3*i+1|1
start = (k*k+4*k-2*k*(i&1))//3
sieve[(k*k)//3::2*k]=zero*((size - (k*k)//3 - 1) // (2 * k) + 1)
sieve[ start ::2*k]=zero*((size - start - 1) // (2 * k) + 1)
ans = [2,3]
poss = chain(izip(*[range(i, n, 6) for i in (1,5)]))
ans.extend(compress(poss, sieve))
return ans
Comparaisons:
>>> timeit.timeit('primes.rwh_primes2(10**6)', setup='import primes', number=1)
0.0652179726976101
>>> timeit.timeit('primes.rwh_primes2_python3(10**6)', setup='import primes', number=1)
0.03267321276325674
et
>>> timeit.timeit('primes.rwh_primes2(10**8)', setup='import primes', number=1)
6.394284538007014
>>> timeit.timeit('primes.rwh_primes2_python3(10**8)', setup='import primes', number=1)
3.833829450302801
Il est instructif d'écrire votre propre code de premier résultat, mais il est aussi utile d'avoir une bibliothèque rapide et fiable à portée de main. J'ai écrit un wrapper autour de la bibliothèque C ++ primesieve, a nommé primesieve-python
Essayez-pip install primesieve
import primesieve
primes = primesieve.generate_primes(10**8)
Je serais curieux de voir la vitesse par rapport.
Voici deux versions mises à jour (Python pur 3.6) de l'une des fonctions les plus rapides,
from itertools import compress
def rwh_primes1v1(n):
""" Returns a list of primes < n for n > 2 """
sieve = bytearray([True]) * (n//2)
for i in range(3,int(n**0.5)+1,2):
if sieve[i//2]:
sieve[i*i//2::i] = bytearray((n-i*i-1)//(2*i)+1)
return [2,*compress(range(3,n,2), sieve[1:])]
def rwh_primes1v2(n):
""" Returns a list of primes < n for n > 2 """
sieve = bytearray([True]) * (n//2+1)
for i in range(1,int(n**0.5)//2+1):
if sieve[i]:
sieve[2*i*(i+1)::2*i+1] = bytearray((n//2-2*i*(i+1))//(2*i+1)+1)
return [2,*compress(range(3,n,2), sieve[1:])]
Une mise en œuvre déterministe du test de primalité de Miller-Rabin sur l'hypothèse que N <9.080.191
import sys
import random
def miller_rabin_pass(a, n):
d = n - 1
s = 0
while d % 2 == 0:
d >>= 1
s += 1
a_to_power = pow(a, d, n)
if a_to_power == 1:
return True
for i in xrange(s-1):
if a_to_power == n - 1:
return True
a_to_power = (a_to_power * a_to_power) % n
return a_to_power == n - 1
def miller_rabin(n):
for a in [2, 3, 37, 73]:
if not miller_rabin_pass(a, n):
return False
return True
n = int(sys.argv[1])
primes = [2]
for p in range(3,n,2):
if miller_rabin(p):
primes.append(p)
print len(primes)
Selon l'article sur Wikipedia ( http: //en.wikipedia. org / wiki / Miller-Rabin_primality_test ) essai N <9080191 pour a = 2,3,37 et 73 est suffisant pour décider si N est composite ou non.
J'adapté le code source de l'application probabiliste du test de Miller-Rabin originale trouvé ici: http://en.literateprograms.org/Miller-Rabin_primality_test_ (python)
Si vous avez le contrôle sur N, le très moyen le plus rapide à la liste tous les nombres premiers est de les précalculer. Sérieusement. Précalcul est un moyen d'optimisation négligé.
Voici le code que je l'habitude d'utiliser pour générer des nombres premiers en Python:
$ python -mtimeit -s'import sieve' 'sieve.sieve(1000000)'
10 loops, best of 3: 445 msec per loop
$ cat sieve.py
from math import sqrt
def sieve(size):
prime=[True]*size
rng=xrange
limit=int(sqrt(size))
for i in rng(3,limit+1,+2):
if prime[i]:
prime[i*i::+i]=[False]*len(prime[i*i::+i])
return [2]+[i for i in rng(3,size,+2) if prime[i]]
if __name__=='__main__':
print sieve(100)
Il ne peut pas rivaliser avec les solutions plus rapides publiées ici, mais au moins il est python pur.
Merci de poster cette question. J'ai vraiment beaucoup appris aujourd'hui.
Pour le code le plus rapide, la solution numpy est le meilleur. Pour des raisons purement académiques, cependant, je poste ma version python pur, qui est un peu moins de 50% plus rapide que la version livre de recettes affiché ci-dessus. Depuis que je fais la liste entière en mémoire, vous avez besoin d'assez d'espace pour contenir tout, mais il semble à l'échelle assez bien.
def daniel_sieve_2(maxNumber):
"""
Given a number, returns all numbers less than or equal to
that number which are prime.
"""
allNumbers = range(3, maxNumber+1, 2)
for mIndex, number in enumerate(xrange(3, maxNumber+1, 2)):
if allNumbers[mIndex] == 0:
continue
# now set all multiples to 0
for index in xrange(mIndex+number, (maxNumber-3)/2+1, number):
allNumbers[index] = 0
return [2] + filter(lambda n: n!=0, allNumbers)
Et les résultats:
>>>mine = timeit.Timer("daniel_sieve_2(1000000)",
... "from sieves import daniel_sieve_2")
>>>prev = timeit.Timer("get_primes_erat(1000000)",
... "from sieves import get_primes_erat")
>>>print "Mine: {0:0.4f} ms".format(min(mine.repeat(3, 1))*1000)
Mine: 428.9446 ms
>>>print "Previous Best {0:0.4f} ms".format(min(prev.repeat(3, 1))*1000)
Previous Best 621.3581 ms
Une mise en œuvre légèrement différente d'un tamis en utilisant la moitié Numpy:
import math
import numpy
def prime6(upto):
primes=numpy.arange(3,upto+1,2)
isprime=numpy.ones((upto-1)/2,dtype=bool)
for factor in primes[:int(math.sqrt(upto))]:
if isprime[(factor-2)/2]: isprime[(factor*3-2)/2:(upto-1)/2:factor]=0
return numpy.insert(primes[isprime],0,2)
Quelqu'un peut-il comparer avec les autres horaires? Sur ma machine, il semble assez comparable à l'autre Numpy demi-tamis.
Il est tout écrit et testé. Il n'y a donc pas besoin de réinventer la roue.
python -m timeit -r10 -s"from sympy import sieve" "primes = list(sieve.primerange(1, 10**6))"
nous donne un record 12.2 msec
10 loops, best of 10: 12.2 msec per loop
Si cela ne suffit pas rapide, vous pouvez essayer PyPy:
pypy -m timeit -r10 -s"from sympy import sieve" "primes = list(sieve.primerange(1, 10**6))"
qui se traduit par:
10 loops, best of 10: 2.03 msec per loop
La réponse avec 247 voix mises liste 15,9 ms pour la meilleure solution. Comparer !!!
Je l'ai testé les fonctions de certains unutbu , je l'ai calculé à l'aide des millions nombre faim
Les gagnants sont les fonctions qui utilisent la bibliothèque numpy,
Remarque : il serait aussi intéressant un test d'utilisation de la mémoire:)
Exemple de code
Code complet sur mon dépôt GitHub
#!/usr/bin/env python
import lib
import timeit
import sys
import math
import datetime
import prettyplotlib as ppl
import numpy as np
import matplotlib.pyplot as plt
from prettyplotlib import brewer2mpl
primenumbers_gen = [
'sieveOfEratosthenes',
'ambi_sieve',
'ambi_sieve_plain',
'sundaram3',
'sieve_wheel_30',
'primesfrom3to',
'primesfrom2to',
'rwh_primes',
'rwh_primes1',
'rwh_primes2',
]
def human_format(num):
# https://stackoverflow.com/questions/579310/formatting-long-numbers-as-strings-in-python?answertab=active#tab-top
magnitude = 0
while abs(num) >= 1000:
magnitude += 1
num /= 1000.0
# add more suffixes if you need them
return '%.2f%s' % (num, ['', 'K', 'M', 'G', 'T', 'P'][magnitude])
if __name__=='__main__':
# Vars
n = 10000000 # number itereration generator
nbcol = 5 # For decompose prime number generator
nb_benchloop = 3 # Eliminate false positive value during the test (bench average time)
datetimeformat = '%Y-%m-%d %H:%M:%S.%f'
config = 'from __main__ import n; import lib'
primenumbers_gen = {
'sieveOfEratosthenes': {'color': 'b'},
'ambi_sieve': {'color': 'b'},
'ambi_sieve_plain': {'color': 'b'},
'sundaram3': {'color': 'b'},
'sieve_wheel_30': {'color': 'b'},
# # # 'primesfrom2to': {'color': 'b'},
'primesfrom3to': {'color': 'b'},
# 'rwh_primes': {'color': 'b'},
# 'rwh_primes1': {'color': 'b'},
'rwh_primes2': {'color': 'b'},
}
# Get n in command line
if len(sys.argv)>1:
n = int(sys.argv[1])
step = int(math.ceil(n / float(nbcol)))
nbs = np.array([i * step for i in range(1, int(nbcol) + 1)])
set2 = brewer2mpl.get_map('Paired', 'qualitative', 12).mpl_colors
print datetime.datetime.now().strftime(datetimeformat)
print("Compute prime number to %(n)s" % locals())
print("")
results = dict()
for pgen in primenumbers_gen:
results[pgen] = dict()
benchtimes = list()
for n in nbs:
t = timeit.Timer("lib.%(pgen)s(n)" % locals(), setup=config)
execute_times = t.repeat(repeat=nb_benchloop,number=1)
benchtime = np.mean(execute_times)
benchtimes.append(benchtime)
results[pgen] = {'benchtimes':np.array(benchtimes)}
fig, ax = plt.subplots(1)
plt.ylabel('Computation time (in second)')
plt.xlabel('Numbers computed')
i = 0
for pgen in primenumbers_gen:
bench = results[pgen]['benchtimes']
avgs = np.divide(bench,nbs)
avg = np.average(bench, weights=nbs)
# Compute linear regression
A = np.vstack([nbs, np.ones(len(nbs))]).T
a, b = np.linalg.lstsq(A, nbs*avgs)[0]
# Plot
i += 1
#label="%(pgen)s" % locals()
#ppl.plot(nbs, nbs*avgs, label=label, lw=1, linestyle='--', color=set2[i % 12])
label="%(pgen)s avg" % locals()
ppl.plot(nbs, a * nbs + b, label=label, lw=2, color=set2[i % 12])
print datetime.datetime.now().strftime(datetimeformat)
ppl.legend(ax, loc='upper left', ncol=4)
# Change x axis label
ax.get_xaxis().get_major_formatter().set_scientific(False)
fig.canvas.draw()
labels = [human_format(int(item.get_text())) for item in ax.get_xticklabels()]
ax.set_xticklabels(labels)
ax = plt.gca()
plt.show()
Pour Python 3
def rwh_primes2(n):
correction = (n%6>1)
n = {0:n,1:n-1,2:n+4,3:n+3,4:n+2,5:n+1}[n%6]
sieve = [True] * (n//3)
sieve[0] = False
for i in range(int(n**0.5)//3+1):
if sieve[i]:
k=3*i+1|1
sieve[ ((k*k)//3) ::2*k]=[False]*((n//6-(k*k)//6-1)//k+1)
sieve[(k*k+4*k-2*k*(i&1))//3::2*k]=[False]*((n//6-(k*k+4*k-2*k*(i&1))//6-1)//k+1)
return [2,3] + [3*i+1|1 for i in range(1,n//3-correction) if sieve[i]]
La première fois en utilisant python, donc quelques-unes des méthodes que j'utilise dans cela peut sembler un peu lourd. Je me suis converti tout droit de mon c ++ code python et c'est ce que j'ai (mais un peu peu slowww en python)
#!/usr/bin/env python
import time
def GetPrimes(n):
Sieve = [1 for x in xrange(n)]
Done = False
w = 3
while not Done:
for q in xrange (3, n, 2):
Prod = w*q
if Prod < n:
Sieve[Prod] = 0
else:
break
if w > (n/2):
Done = True
w += 2
return Sieve
start = time.clock()
d = 10000000
Primes = GetPrimes(d)
count = 1 #This is for 2
for x in xrange (3, d, 2):
if Primes[x]:
count+=1
elapsed = (time.clock() - start)
print "\nFound", count, "primes in", elapsed, "seconds!\n"
pythonw Primes.py
Trouvé 664579 nombres premiers 12.799119 secondes!
#!/usr/bin/env python
import time
def GetPrimes2(n):
Sieve = [1 for x in xrange(n)]
for q in xrange (3, n, 2):
k = q
for y in xrange(k*3, n, k*2):
Sieve[y] = 0
return Sieve
start = time.clock()
d = 10000000
Primes = GetPrimes2(d)
count = 1 #This is for 2
for x in xrange (3, d, 2):
if Primes[x]:
count+=1
elapsed = (time.clock() - start)
print "\nFound", count, "primes in", elapsed, "seconds!\n"
pythonw Primes2.py
Trouvé 664579 nombres premiers 10.230172 secondes!
#!/usr/bin/env python
import time
def GetPrimes3(n):
Sieve = [1 for x in xrange(n)]
for q in xrange (3, n, 2):
k = q
for y in xrange(k*k, n, k << 1):
Sieve[y] = 0
return Sieve
start = time.clock()
d = 10000000
Primes = GetPrimes3(d)
count = 1 #This is for 2
for x in xrange (3, d, 2):
if Primes[x]:
count+=1
elapsed = (time.clock() - start)
print "\nFound", count, "primes in", elapsed, "seconds!\n"
python Primes2.py
Trouvé 664579 nombres premiers en 7.113776 secondes!
Je sais que le concours est fermé depuis quelques années.…
Néanmoins, c'est ma suggestion pour un tamis principal en python pur, basé sur l'omission des multiples de 2, 3 et 5 en utilisant les étapes appropriées lors du traitement du tamis vers l'avant.Néanmoins, il est en fait plus lent pour N<10^9 que les solutions supérieures de @Robert William Hanks rwh_primes2 et rwh_primes1.En utilisant un tableau de tamis ctypes.c_ushort supérieur à 1,5* 10^8, il s'adapte en quelque sorte aux limites de mémoire.
10^6
$ Python -mtimeit -s "Importer PrimesievespeedCompe" "PrimesievespeedComp.primesieveseq (1000000)" 10 boucles, Best of 3:46,7 ms par boucle
Pour comparer: $ python -mtimeit -s "Import PrimesievespeedComp" "PrimesievespeedComp.rwh_primes1 (1000000)" 10 boucles, meilleur sur 3:43,2 msec par boucle pour comparer:$ python -m timeit -s "Importer PrimesievespeedCompe" "PrimesievespeedComp.rwh_primes2 (1000000)" 10 boucles, meilleur sur 3:34,5 ms par boucle
10^7
$ python -mtimeit -s "Import PrimesievespeedCompe" "PrimesievespeedComp.primesievseq (10000000)" 10 boucles, Best of 3:530 ms par boucle
Pour comparer: $ python -mtimeit -s "Importer PrimesievespeedComp" "PrimesievespeedComp.rwh_primes1 (10000000)" 10 boucles, Best of 3:494 ms par boucle pour comparer:$ python -m timeit -s "Importer PrimesievespeedComp" "PrimesievespeedComp.rwh_primes2 (10000000)" 10 boucles, Best of 3:375 ms par boucle
10^8
$ Python -mtimeit -s "Importer PrimesievespeedCompe" "PrimesievespeedComp.primesieveseq (100000000)" 10 boucles, Best of 3:5,55 secondes par boucle
comparer:$ Python -Mtimeit -s "Importer PrimesievespeedComp" "PrimesievespeedComp.rwh_primes1 (100000000)" 10 boucles, Best of 3:5,33 secondes par boucle pour comparer:$ python -m timeit -s "Importer PrimesievespeedComp" "PrimesievespeedComp.rwh_primes2 (100000000)" 10 boucles, meilleur sur 3:3,95 secondes par boucle
10^9
$ Python -mtimeit -s "Importer PrimesievespeedCompe" "PrimesievespeedComp.primesieveseq (1000000000)" 10 boucles, Best of 3: 61.2 sec par boucle
comparer:$ Python -Mtimeit -N 3 -S "Importer PrimesievespeedComp" "PrimesievespeedComp.rwh_primes1 (1000000000)" 3 boucles, Best of 3: 97.8sec par boucle
comparer:$ python -m timeit -s "Importer PrimesievespeedCompe" "PrimesievespeedComp.rwh_primes2 (1000000000)" 10 boucles, meilleur sur 3:41,9 secondes par boucle
Vous pouvez copier le code ci-dessous dans Ubuntu primeSieveSpeedComp pour revoir ces tests.
def primeSieveSeq(MAX_Int):
if MAX_Int > 5*10**8:
import ctypes
int16Array = ctypes.c_ushort * (MAX_Int >> 1)
sieve = int16Array()
#print 'uses ctypes "unsigned short int Array"'
else:
sieve = (MAX_Int >> 1) * [False]
#print 'uses python list() of long long int'
if MAX_Int < 10**8:
sieve[4::3] = [True]*((MAX_Int - 8)/6+1)
sieve[12::5] = [True]*((MAX_Int - 24)/10+1)
r = [2, 3, 5]
n = 0
for i in xrange(int(MAX_Int**0.5)/30+1):
n += 3
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
n += 2
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
n += 1
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
n += 2
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
n += 1
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
n += 2
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
n += 3
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
n += 1
if not sieve[n]:
n2 = (n << 1) + 1
r.append(n2)
n2q = (n2**2) >> 1
sieve[n2q::n2] = [True]*(((MAX_Int >> 1) - n2q - 1) / n2 + 1)
if MAX_Int < 10**8:
return [2, 3, 5]+[(p << 1) + 1 for p in [n for n in xrange(3, MAX_Int >> 1) if not sieve[n]]]
n = n >> 1
try:
for i in xrange((MAX_Int-2*n)/30 + 1):
n += 3
if not sieve[n]:
r.append((n << 1) + 1)
n += 2
if not sieve[n]:
r.append((n << 1) + 1)
n += 1
if not sieve[n]:
r.append((n << 1) + 1)
n += 2
if not sieve[n]:
r.append((n << 1) + 1)
n += 1
if not sieve[n]:
r.append((n << 1) + 1)
n += 2
if not sieve[n]:
r.append((n << 1) + 1)
n += 3
if not sieve[n]:
r.append((n << 1) + 1)
n += 1
if not sieve[n]:
r.append((n << 1) + 1)
except:
pass
return r
tamis le plus rapide de choix dans pur Python :
from itertools import compress
def half_sieve(n):
"""
Returns a list of prime numbers less than `n`.
"""
if n <= 2:
return []
sieve = bytearray([True]) * (n // 2)
for i in range(3, int(n ** 0.5) + 1, 2):
if sieve[i // 2]:
sieve[i * i // 2::i] = bytearray((n - i * i - 1) // (2 * i) + 1)
primes = list(compress(range(1, n, 2), sieve))
primes[0] = 2
return primes
I optimisé Crible d'Ératosthène pour la vitesse et la mémoire.
Indice de référence
from time import clock
import platform
def benchmark(iterations, limit):
start = clock()
for x in range(iterations):
half_sieve(limit)
end = clock() - start
print(f'{end/iterations:.4f} seconds for primes < {limit}')
if __name__ == '__main__':
print(platform.python_version())
print(platform.platform())
print(platform.processor())
it = 10
for pw in range(4, 9):
benchmark(it, 10**pw)
Sortie
>>> 3.6.7
>>> Windows-10-10.0.17763-SP0
>>> Intel64 Family 6 Model 78 Stepping 3, GenuineIntel
>>> 0.0003 seconds for primes < 10000
>>> 0.0021 seconds for primes < 100000
>>> 0.0204 seconds for primes < 1000000
>>> 0.2389 seconds for primes < 10000000
>>> 2.6702 seconds for primes < 100000000
Voici une version numpy de Sieve d'Eratosthène ayant à la fois une bonne complexité (inférieure à trier un tableau de longueur n) et vectorisation. Par rapport aux temps de @unutbu ce aussi vite que les paquets avec 46 microsecons pour trouver tous les nombres premiers inférieurs à un million.
import numpy as np
def generate_primes(n):
is_prime = np.ones(n+1,dtype=bool)
is_prime[0:2] = False
for i in range(int(n**0.5)+1):
if is_prime[i]:
is_prime[i**2::i]=False
return np.where(is_prime)[0]
synchronisations:
import time
for i in range(2,10):
timer =time.time()
generate_primes(10**i)
print('n = 10^',i,' time =', round(time.time()-timer,6))
>> n = 10^ 2 time = 5.6e-05
>> n = 10^ 3 time = 6.4e-05
>> n = 10^ 4 time = 0.000114
>> n = 10^ 5 time = 0.000593
>> n = 10^ 6 time = 0.00467
>> n = 10^ 7 time = 0.177758
>> n = 10^ 8 time = 1.701312
>> n = 10^ 9 time = 19.322478
Je suppose que le le plus rapide de toutes façons est de coder en dur les nombres premiers dans votre code.
Alors pourquoi ne pas simplement écrire un script lent qui génère un autre fichier source qui a tous les numéros cablé dans, puis importer ce fichier source lorsque vous exécutez votre programme réel.
Bien sûr, cela ne fonctionne que si vous connaissez la limite supérieure de N au moment de la compilation, mais est donc le cas pour (presque) tous les projets des problèmes d'Euler.
PS: Je peux me tromper si ssi analyse la source avec les nombres premiers câblés est plus lent que les calculer en premier lieu, mais aussi loin je sais que Python va de fichiers .pyc compilés si la lecture d'un tableau binaire avec tous les nombres premiers jusqu'à N devrait être sanglante rapide dans ce cas.
Désolé de déranger, mais erat2 () a une grave lacune dans l'algorithme.
Alors que la recherche de la prochaine composite, nous avons besoin de tester les nombres impairs seulement. q, p sont tous deux impair; alors q + p est même et n'a pas besoin d'être testé, mais q + 2 * p est toujours impair. Ceci élimine le test « si même » dans la condition de la boucle while et permet d'économiser environ 30% du temps d'exécution.
Alors que nous y sommes: au lieu de l'élégant 'D.pop (q, None)' obtenir et supprimer l'utilisation de la méthode 'si q D: p = D [q], del D [q]' qui est deux fois plus vite! Au moins sur ma machine (P3-1Ghz). Je suggère donc que cette mise en œuvre de cet algorithme intelligent:
def erat3( ):
from itertools import islice, count
# q is the running integer that's checked for primeness.
# yield 2 and no other even number thereafter
yield 2
D = {}
# no need to mark D[4] as we will test odd numbers only
for q in islice(count(3),0,None,2):
if q in D: # is composite
p = D[q]
del D[q]
# q is composite. p=D[q] is the first prime that
# divides it. Since we've reached q, we no longer
# need it in the map, but we'll mark the next
# multiple of its witnesses to prepare for larger
# numbers.
x = q + p+p # next odd(!) multiple
while x in D: # skip composites
x += p+p
D[x] = p
else: # is prime
# q is a new prime.
# Yield it and mark its first multiple that isn't
# already marked in previous iterations.
D[q*q] = q
yield q
La méthode la plus rapide que je l'ai essayé jusqu'à présent est basé sur le livre de recettes Python erat2 fonction :
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:
x = q + 2*p
while x in D:
x += 2*p
D[x] = p
Voir cette répondre pour une explication de l'accélération.
Je peux être en retard à la fête, mais devra ajouter mon propre code pour cela. Il utilise environ n / 2 dans l'espace parce que nous ne avons pas besoin de stocker les numéros pairs et je fais également l'utilisation du module python BitArray, réduisant encore de façon drastique la consommation de mémoire et permettant le calcul de tous les nombres premiers jusqu'à 1000000000
from bitarray import bitarray
def primes_to(n):
size = n//2
sieve = bitarray(size)
sieve.setall(1)
limit = int(n**0.5)
for i in range(1,limit):
if sieve[i]:
val = 2*i+1
sieve[(i+i*val)::val] = 0
return [2] + [2*i+1 for i, v in enumerate(sieve) if v and i > 0]
python -m timeit -n10 -s "import euler" "euler.primes_to(1000000000)"
10 loops, best of 3: 46.5 sec per loop
Ceci a été exécuté sur un 2.4Ghz 64bit MAC OSX 10.8.3
Je collectionnais plusieurs tamis de nombres premiers au fil du temps. Le plus rapide sur mon ordinateur est ceci:
from time import time
# 175 ms for all the primes up to the value 10**6
def primes_sieve(limit):
a = [True] * limit
a[0] = a[1] = False
#a[2] = True
for n in xrange(4, limit, 2):
a[n] = False
root_limit = int(limit**.5)+1
for i in xrange(3,root_limit):
if a[i]:
for n in xrange(i*i, limit, 2*i):
a[n] = False
return a
LIMIT = 10**6
s=time()
primes = primes_sieve(LIMIT)
print time()-s
Je suis lent à répondre à cette question, mais il me semblait comme un exercice amusant. J'utilise numpy qui pourrait être la tricherie et je doute cette méthode est la plus rapide, mais il devrait être clair. Il tamise un tableau booléen se référant à ses indices uniquement et suscite des nombres premiers à partir des indices de toutes les valeurs vraies. Pas modulo nécessaire.
import numpy as np
def ajs_primes3a(upto):
mat = np.ones((upto), dtype=bool)
mat[0] = False
mat[1] = False
mat[4::2] = False
for idx in range(3, int(upto ** 0.5)+1, 2):
mat[idx*2::idx] = False
return np.where(mat == True)[0]
Voici une technique intéressante pour générer des nombres premiers (mais pas le plus efficace) à l'aide des compréhensions liste de python:
noprimes = [j for i in range(2, 8) for j in range(i*2, 50, i)]
primes = [x for x in range(2, 50) if x not in noprimes]
Vous pouvez trouver l'exemple et quelques explications à droite href="http://www.secnetix.de/olli/Python/list_comprehensions.hawk"
J'ai trouvé un pur Python 2 premier générateur , dans un commentaire par Willy Bon, qui est plus rapide que rwh2_primes.
def primes235(limit):
yield 2; yield 3; yield 5
if limit < 7: return
modPrms = [7,11,13,17,19,23,29,31]
gaps = [4,2,4,2,4,6,2,6,4,2,4,2,4,6,2,6] # 2 loops for overflow
ndxs = [0,0,0,0,1,1,2,2,2,2,3,3,4,4,4,4,5,5,5,5,5,5,6,6,7,7,7,7,7,7]
lmtbf = (limit + 23) // 30 * 8 - 1 # integral number of wheels rounded up
lmtsqrt = (int(limit ** 0.5) - 7)
lmtsqrt = lmtsqrt // 30 * 8 + ndxs[lmtsqrt % 30] # round down on the wheel
buf = [True] * (lmtbf + 1)
for i in xrange(lmtsqrt + 1):
if buf[i]:
ci = i & 7; p = 30 * (i >> 3) + modPrms[ci]
s = p * p - 7; p8 = p << 3
for j in range(8):
c = s // 30 * 8 + ndxs[s % 30]
buf[c::p8] = [False] * ((lmtbf - c) // p8 + 1)
s += p * gaps[ci]; ci += 1
for i in xrange(lmtbf - 6 + (ndxs[(limit - 7) % 30])): # adjust for extras
if buf[i]: yield (30 * (i >> 3) + modPrms[i & 7])
Mes résultats:
$ time ./prime_rwh2.py 1e8
5761455 primes found < 1e8
real 0m3.201s
user 0m2.609s
sys 0m0.578s
$ time ./prime_wheel.py 1e8
5761455 primes found < 1e8
real 0m2.710s
user 0m2.469s
sys 0m0.219s
... sur mon ordinateur portable récent de milieu de gamme (i5 8265U cadencé à 1,6 GHz) en cours d'exécution Ubuntu sur Windows 10.
Ceci est un mod tamis de 30 roues qui sautille multiples de 2, 3 et 5. Il fonctionne très bien pour moi jusqu'à environ 2.5e9, quand mon portable commence à court de 8G RAM et échange beaucoup.
J'aime mod 30 car il a seulement 8 reliquats qui ne sont pas des multiples de 2, 3 ou 5. Cela permet d'utiliser des quarts de travail et pour la multiplication, la division « & » et mod, et devrait permettre des résultats d'emballage pour un mod 30 roue dans un octet. Je suis en train de transformer cela en un tamis segmenté pour réduire l'empreinte RAM / échange pour N vraiment grand, mais n'ont pas encore réussi. La partie la plus difficile est de trouver les multiples premiers de base corrects et indices pour commencer lorsque vous passez au segment suivant.
Il y a un Javascript Version href="https://stackoverflow.com/a/57108107/11943198"> encore plus rapide qui est segmenté et utilise un 210 mod roue (pas multiple de 2, 3 , 5 ou 7) par @GordonBGood avec une explication approfondie qui est utile pour moi. Celui-là est si rapide que je dois apprendre Javascript. Si j'ignorer tous les « nous » s il semble assez lisible.
