rango - programa que muestre los numeros primos en python
La forma más rápida de enumerar todos los números primos debajo de N (30)
Aquí está el código que normalmente uso para generar números primos 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)
No puede competir con las soluciones más rápidas publicadas aquí, pero al menos es pura python.
Gracias por publicar esta pregunta. Realmente aprendí mucho hoy.
Este es el mejor algoritmo que pude encontrar.
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
¿Se puede hacer aún más rápido?
Este código tiene un defecto: dado que los numbers son un conjunto desordenado, no hay garantía de que numbers.pop() eliminará el número más bajo del conjunto. Sin embargo, funciona (al menos para mí) para algunos números de entrada:
>>> 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
Aquí hay dos versiones actualizadas (de Python 3.6 puro) de una de las funciones más rápidas,
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:])]
Aquí hay una muestra bastante clara del libro de cocina de Python: la versión más rápida propuesta en esa URL es:
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
por lo que daría
def get_primes_erat(n):
return list(itertools.takewhile(lambda p: p<n, erat2()))
Al medir en el indicador de shell (como prefiero hacer) con este código en pri.py, observo:
$ 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
Así que parece que la solución de Cookbook es más del doble de rápida.
El algoritmo es rápido, pero tiene un defecto grave:
>>> 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
Usted asume que numbers.pop() devolvería el número más pequeño en el conjunto, pero esto no está garantizado en absoluto. Los conjuntos no están ordenados y pop() elimina y devuelve un elemento arbitrary , por lo que no se puede usar para seleccionar el siguiente primo de los números restantes.
Es instructivo escribir su propio código de búsqueda principal, pero también es útil tener una biblioteca rápida y confiable a mano. Escribí una envoltura alrededor de la biblioteca de C ++ primesieve , la llamé primesieve-python
Pruébalo pip install primesieve
import primesieve
primes = primesieve.generate_primes(10**8)
Tendría curiosidad por ver la velocidad comparada.
Para el código más rápido, la solución numpy es la mejor. Sin embargo, por razones puramente académicas, estoy publicando mi versión de python pura, que es un poco menos del 50% más rápida que la versión de libro de recetas publicada anteriormente. Ya que hago la lista completa en la memoria, necesitas suficiente espacio para guardar todo, pero parece que se escala bastante 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)
Y los resultados:
>>>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
Pregunta relacionada (que trata con los generadores de primos y que incluye puntos de referencia):
¿Acelerar las operaciones de cadena de bits / bit en Python?
Para las versiones de python 3 que utilizan bytearray y compress, consulte: https://.com/a/46635266/350331
Más rápido y más código de Python puro en cuanto a memoria:
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]]
o comenzando con medio tamiz
def primes1(n):
""" 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]]
Un código numpy más rápido y con más memoria:
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
una variación más rápida a partir de un tercio de un tamiz:
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 xrange(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)]
Una versión de código puro (difícil de codificar) del código anterior sería:
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 xrange(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 xrange(1,n/3-correction) if sieve[i]]
Desafortunadamente, pure-python no adopta la forma más simple y rápida de realizar Asignaciones y llamar a len () dentro del bucle como en [False]*len(sieve[((k*k)/3)::2*k]) es demasiado lento. Así que tuve que improvisar para corregir la entrada (y evitar más matemáticas) y hacer algo de magia mágica extrema (y dolorosa).
Personalmente, creo que es una pena que numpy (que se usa tanto) no sea parte de la biblioteca estándar de Python (2 años después de la versión de Python 3 y no tenga compatibilidad con numpy), y que las mejoras en la sintaxis y la velocidad parecen ser pasadas por alto por completo. Desarrolladores de Python.
Si lo aceptas, pero no por cantidad, aquí hay una adaptación de rwh_primes2 para Python 3 que se ejecuta casi el doble de rápido en mi máquina. El único cambio sustancial es usar un bytearray en lugar de una lista para el booleano, y usar comprimir en lugar de una lista de comprensión para construir la lista final. (Agregaría esto como un comentario como moarningsun si pudiera).
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
Comparaciones
>>> 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
y
>>> 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
Si no desea reinventar la rueda, puede instalar la biblioteca simbólica de matemáticas sympy (sí, es compatible con Python 3)
pip install sympy
Y usa la función de primer primerange .
from sympy import sieve
primes = list(sieve.primerange(1, 10**6))
Si tiene control sobre N, la forma más rápida de enumerar todos los números primos es precalcularlos. Seriamente. La precomputación es una optimización pasada por alto.
Una implementación determinista de la prueba de primalidad de Miller-Rabin en el supuesto de 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)
Según el artículo en Wikipedia ( http://en.wikipedia.org/wiki/Miller–Rabin_primality_test ), la prueba N <9,080,191 para a = 2,3,37, y 73 es suficiente para decidir si N es compuesto o no.
Y adapté el código fuente de la implementación probabilística de la prueba original de Miller-Rabin que se encuentra aquí: http://en.literateprograms.org/Miller-Rabin_primality_test_(Python)
Una solución verdaderamente más rápida con una N suficientemente grande sería descargar una lista de números primos calculada previamente , almacenarla como una tupla y hacer algo como:
for pos,i in enumerate(primes):
if i > N:
print primes[:pos]
Si N > primes[-1] solo luego calcule más primos y guarde la nueva lista en su código, de modo que la próxima vez sea igual de rápido.
Siempre piensa fuera de la caja.
Usando Sundaram''s Sieve , creo que rompí el récord de puro 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
Advertencia: los resultados de timeit pueden variar debido a diferencias en el hardware o la versión de Python.
A continuación se muestra un script que compara una serie de implementaciones:
- ambi_sieve_plain,
- rwh_primes ,
- rwh_primes ,
- rwh_primes ,
- sieveOfAtkin ,
- sieveOfEratosthenes ,
- sundaram3
- sieve_wheel_30 ,
- ambi_sieve (requiere numpy)
- rwh_primes (requiere numpy)
- rwh_primes (requiere numpy)
Muchas gracias a stephan por traerme sieve_wheel_30 a mi atención. El crédito va a rwh_primes para primesfrom2to, primesfrom3to, rwh_primes, rwh_primes1, y rwh_primes2.
De los métodos de Python comprobados, con psyco , para n = 1000000, rwh_primes1 fue el más rápido probado.
+---------------------+-------+
| 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 |
+---------------------+-------+
De los métodos de Python probados, sin psyco , para n = 1000000, rwh_primes2 fue el más rápido.
+---------------------+-------+
| 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 todos los métodos probados, que permiten numpy , para n = 1000000, primesfrom2to fue el más rápido probado.
+---------------------+-------+
| Method | ms |
+---------------------+-------+
| primesfrom2to | 15.9 |
| primesfrom3to | 18.4 |
| ambi_sieve | 29.3 |
+---------------------+-------+
Los tiempos fueron medidos usando el comando:
python -mtimeit -s"import primes" "primes.{method}(1000000)"
con {method} reemplazado por cada uno de los nombres de los métodos.
primes.py:
#!/usr/bin/env python
import psyco; psyco.full()
from math import sqrt, ceil
import numpy as np
def rwh_primes(n):
# https://.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://.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://.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: <[email protected]>, 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, <[email protected]>, 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://.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://.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://.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)
Ejecutando el script prueba que todas las implementaciones dan el mismo resultado.
Aquí hay una técnica interesante para generar números primos (pero no el más eficiente) usando la lista de comprensión 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]
Puedes encontrar el ejemplo y algunas explicaciones here
Esta es la forma en que se puede comparar con los demás.
# You have to list primes upto n
nums = xrange(2, n)
for i in range(2, 10):
nums = filter(lambda s: s==i or s%i, nums)
print nums
Tan sencillo...
Lo siento molestar pero erat2 () tiene una falla seria en el algoritmo.
Mientras buscamos el siguiente compuesto, solo necesitamos probar los números impares. q, p ambos son impares; entonces q + p es par y no necesita ser probado, pero q + 2 * p siempre es impar. Esto elimina la prueba "if even" en la condición de bucle while y ahorra aproximadamente el 30% del tiempo de ejecución.
Mientras estamos en ello: en lugar del elegante método ''D.pop (q, None)'' obtener y eliminar use ''si q en D: p = D [q], del D [q]'', que es el doble de rápido ! Al menos en mi máquina (P3-1GHz). Así que sugiero esta implementación de este algoritmo inteligente:
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
Para 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]]
Sé que la competencia está cerrada por algunos años. ...
No obstante, esta es mi sugerencia para un tamiz primario puro de python, basado en omitir los múltiplos de 2, 3 y 5 utilizando los pasos apropiados al procesar el tamiz hacia adelante. No obstante, en realidad es más lento para N <10 ^ 9 que las soluciones superiores de @Robert William Hanks rwh_primes2 y rwh_primes1. Al usar una matriz de tamices ctypes.c_ushort por encima de 1.5 * 10 ^ 8, de alguna manera se adapta a los límites de memoria.
10 ^ 6
$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (1000000)" 10 bucles, lo mejor de 3: 46.7 mseg por bucle
para comparar: $ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (1000000)" 10 bucles, lo mejor de 3: 43.2 seg por bucle para comparar: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeed (1000000) "10 bucles, lo mejor de 3: 34.5 mseg por bucle
10 ^ 7
$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (10000000)" 10 bucles, lo mejor de 3: 530 mseg por bucle
para comparar: $ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (10000000)" 10 bucles, lo mejor de 3: 494 ms por bucle para comparar: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeed (10000000) "10 bucles, lo mejor de 3: 375 mseg por bucle
10 ^ 8
$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (100000000)" 10 bucles, el mejor de 3: 5,55 seg por bucle
para comparar: $ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (100000000)" 10 bucles, lo mejor de 3: 5.33 seg por bucle para comparar: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeed (100000000) "10 bucles, el mejor de 3: 3,95 seg por bucle
10 ^ 9
$ python -mtimeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.primeSieveSeq (1000000000)" 10 bucles, el mejor de 3: 61.2 segundos por bucle
para comparar: $ python -mtimeit -n 3 -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes1 (1000000000)" 3 bucles, el mejor de 3: 97.8 seg por bucle
para comparar: $ python -m timeit -s "import primeSieveSpeedComp" "primeSieveSpeedComp.rwh_primes2 (1000000000)" 10 bucles, el mejor de 3: 41.9 seg por bucle
Puede copiar el siguiente código en ubuntus primeSieveSpeedComp para revisar estas pruebas.
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
Una implementación ligeramente diferente de un medio tamiz usando 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)
¿Alguien puede comparar esto con los otros tiempos? En mi máquina parece bastante comparable con la otra media criba Numpy.
active algunas active con el número de millones hambrientos
Los ganadores son las funciones que utilizan la biblioteca numpy,
Nota : También sería interesante hacer una prueba de utilización de memoria :)
Código de muestra
Código completo en mi repositorio 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://.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()
El método más rápido que he probado hasta ahora se basa en la función de here :
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
Vea this respuesta para una explicación de la aceleración.
En general, si necesita un cálculo rápido de números, Python no es la mejor opción. Hoy en día hay muchos algoritmos más rápidos (y complejos). Por ejemplo, en mi computadora obtuve 2.2 segundos para su código, con Mathematica obtuve 0.088005.
En primer lugar: ¿necesitas set?
Esta es una solución elegante y más simple para encontrar números primos utilizando una lista almacenada. Comienza con 4 variables, solo tiene que probar los números primos impares para los divisores, y solo tiene que probar hasta la mitad de qué número está probando como primo (no tiene sentido probar si 9, 11, 13 se dividen en 17) . Prueba primos previamente almacenados como divisores.
# Program to calculate Primes
primes = [1,3,5,7]
for n in range(9,100000,2):
for x in range(1,(len(primes)/2)):
if n % primes[x] == 0:
break
else:
primes.append(n)
print primes
Estoy respondiendo lentamente a esta pregunta, pero me pareció un ejercicio divertido. Estoy usando Numpy, que podría estar haciendo trampa y dudo que este método sea el más rápido, pero debería ser claro. Tamiza una matriz booleana que se refiere solo a sus índices y obtiene números primos de los índices de todos los valores verdaderos. No se necesita módulo.
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]
La primera vez que utilizo python, algunos de los métodos que uso en esto pueden parecer un poco incómodos. Acabo de convertir directamente mi código c ++ a python y esto es lo que tengo (aunque un poco lento 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
¡Encontrado 664579 números primos en 12.799119 segundos!
#!/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
¡Encontrado 664579 números primos en 10.230172 segundos!
#!/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
¡Encontrado 664579 números primos en 7.113776 segundos!
Mi conjetura es que la más rápida de todas las formas es codificar los números primos en su código.
Entonces, ¿por qué no escribir un script lento que genere otro archivo de origen con todos los números conectados y luego importar ese archivo de origen cuando ejecuta su programa real?
Por supuesto, esto funciona solo si conoce el límite superior de N en el momento de la compilación, pero así es el caso de (casi) todos los problemas de Euler del proyecto.
PD: Podría estar equivocado, aunque si el análisis de la fuente con primos fijos es más lento que calcularlos en primer lugar, pero por lo que sé, Python se ejecuta desde .pycarchivos compilados , por lo que leer una matriz binaria con todos los primos hasta N debería ser sangriento Rápido en ese caso.
Puede que llegue tarde a la fiesta, pero tendré que agregar mi propio código para esto. Utiliza aproximadamente n / 2 en el espacio porque no necesitamos almacenar números pares y también utilizo el módulo de bitarray python, reduciendo drásticamente el consumo de memoria y permitiendo el cálculo de todos los primos hasta 1,000,000,000
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
Esto se ejecutó en un 64bit 2.4GHZ MAC OSX 10.8.3
Recolecté varios tamices de números primos con el tiempo. El más rápido en mi computadora es este:
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
Todo está escrito y probado. Así que no hay necesidad de reinventar la rueda.
python -m timeit -r10 -s"from sympy import sieve" "primes = list(sieve.primerange(1, 10**6))"
¡Nos da un récord de 12.2 mseg !
10 loops, best of 10: 12.2 msec per loop
Si esto no es lo suficientemente rápido, puedes probar PyPy:
pypy -m timeit -r10 -s"from sympy import sieve" "primes = list(sieve.primerange(1, 10**6))"
lo que resulta en:
10 loops, best of 10: 2.03 msec per loop
La respuesta con 247 votos ascendentes enumera 15.9 ms para la mejor solución. Compara esto !!!