Tamiz de Eratóstenes – Encontrando Primes Python

Solo para aclarar, esto no es un problema de tarea 🙂

Quería encontrar números primos para una aplicación de matemáticas que estoy construyendo y me encontré con el enfoque de Tamiz de Eratóstenes .

He escrito una implementación de ella en Python. Pero es terriblemente lento. Por ejemplo, si quiero encontrar todos los números primos de menos de 2 millones. Se tarda> 20 minutos. (Lo paré en este punto). ¿Cómo puedo acelerar esto?

def primes_sieve(limit): limitn = limit+1 primes = range(2, limitn) for i in primes: factors = range(i, limitn, i) for f in factors[1:]: if f in primes: primes.remove(f) return primes print primes_sieve(2000) 

ACTUALIZACIÓN: Terminé haciendo perfiles en este código y encontré que se gastó bastante tiempo en eliminar un elemento de la lista. Bastante comprensible, teniendo en cuenta que tiene que atravesar la lista completa (en el peor de los casos) para encontrar el elemento y luego eliminarlo y luego reajustar la lista (¿tal vez alguna copia continúa?). De todos modos, eliminé la lista para el diccionario. Mi nueva implementación –

 def primes_sieve1(limit): limitn = limit+1 primes = dict() for i in range(2, limitn): primes[i] = True for i in primes: factors = range(i,limitn, i) for f in factors[1:]: primes[f] = False return [i for i in primes if primes[i]==True] print primes_sieve1(2000000) 

No estás implementando el algoritmo correcto:

En su primer ejemplo, primes_sieve no mantiene una lista de indicadores de primalidad para marcar / desarmar (como en el algoritmo), sino que cambia el tamaño de una lista de enteros continuamente, lo que es muy costoso: eliminar un elemento de una lista requiere cambiar todos los siguientes artículos por uno.

En el segundo ejemplo, primes_sieve1 mantiene un diccionario de indicadores de primalidad, que es un paso en la dirección correcta, pero itera sobre el diccionario en un orden indefinido, y elimina factores de factores (en lugar de solo factores de primos, como en el algoritmo). Puede solucionar esto clasificando las claves y omitiendo los números primos (lo que ya hace que sea un orden de magnitud más rápido), pero aún es mucho más eficiente usar una lista directamente.

El algoritmo correcto (con una lista en lugar de un diccionario) se ve algo así como:

 def primes_sieve2(limit): a = [True] * limit # Initialize the primality list a[0] = a[1] = False for (i, isprime) in enumerate(a): if isprime: yield i for n in xrange(i*i, limit, i): # Mark factors non-prime a[n] = False 

(Tenga en cuenta que esto también incluye la optimización algorítmica de iniciar el marcado no primo en el cuadrado del primo ( i*i ) en lugar de su doble).

 def eratosthenes(n): multiples = [] for i in range(2, n+1): if i not in multiples: print (i) for j in range(i*i, n+1, i): multiples.append(j) eratosthenes(100) 

Eliminar desde el principio de una matriz (lista) requiere mover todos los elementos después de que esté abajo. Eso significa que eliminar todos los elementos de una lista de esta manera comenzando desde el frente es una operación O (n ^ 2).

Puedes hacer esto mucho más eficientemente con conjuntos:

 def primes_sieve(limit): limitn = limit+1 not_prime = set() primes = [] for i in range(2, limitn): if i in not_prime: continue for f in range(i*2, limitn, i): not_prime.add(f) primes.append(i) return primes print primes_sieve(1000000) 

… o alternativamente, evita tener que reorganizar la lista:

 def primes_sieve(limit): limitn = limit+1 not_prime = [False] * limitn primes = [] for i in range(2, limitn): if not_prime[i]: continue for f in xrange(i*2, limitn, i): not_prime[f] = True primes.append(i) return primes 

Me doy cuenta de que esto no es realmente responder a la pregunta de cómo generar números primos rápidamente, pero tal vez a algunos les resulte interesante esta alternativa: como Python proporciona una evaluación perezosa a través de generadores, el tamiz de eratóstenes se puede implementar exactamente como se indica:

 def intsfrom(n): while True: yield n n += 1 def sieve(ilist): p = next(ilist) yield p for q in sieve(n for n in ilist if n%p != 0): yield q try: for p in sieve(intsfrom(2)): print p, print '' except RuntimeError as e: print e 

El bloque de prueba está allí porque el algoritmo se ejecuta hasta que explota la stack y, sin el bloque de prueba, se muestra el retroceso empujando la salida real que desea ver fuera de la pantalla.

Al combinar las contribuciones de muchos entusiastas (incluidos Glenn Maynard y MrHIDEn de los comentarios anteriores), se me ocurrió el siguiente fragmento de código en Python 2:

 def simpleSieve(sieveSize): #creating Sieve. sieve = [True] * (sieveSize+1) # 0 and 1 are not considered prime. sieve[0] = False sieve[1] = False for i in xrange(2,int(math.sqrt(sieveSize))+1): if sieve[i] == False: continue for pointer in xrange(i**2, sieveSize+1, i): sieve[pointer] = False # Sieve is left with prime numbers == True primes = [] for i in xrange(sieveSize+1): if sieve[i] == True: primes.append(i) return primes sieveSize = input() primes = simpleSieve(sieveSize) 

El tiempo necesario para el cálculo en mi máquina para diferentes entradas con una potencia de 10 es:

  • 3: 0.3 ms
  • 4: 2.4 ms
  • 5: 23 ms
  • 6: 0.26 s
  • 7: 3.1 s
  • 8: 33 s

Mucho mas rápido:

 import time def get_primes(n): m = n+1 #numbers = [True for i in range(m)] numbers = [True] * m #EDIT: faster for i in range(2, int(n**0.5 + 1)): if numbers[i]: for j in range(i*i, m, i): numbers[j] = False primes = [] for i in range(2, m): if numbers[i]: primes.append(i) return primes start = time.time() primes = get_primes(10000) print(time.time() - start) print(get_primes(100)) 

Un truco de velocidad simple: cuando define la variable “números primos”, establezca el paso en 2 para omitir todos los números pares automáticamente y establezca el punto de inicio en 1.

Luego puede optimizar aún más por en lugar de for i en primes, use para i en primes [: round (len (primes) ** 0.5)]. Eso boostá dramáticamente el rendimiento. Además, puedes eliminar los números que terminan con 5 para boost aún más la velocidad.

Mi implementación:

 import math n = 100 marked = {} for i in range(2, int(math.sqrt(n))): if not marked.get(i): for x in range(i * i, n, i): marked[x] = True for i in range(2, n): if not marked.get(i): print i 

Aquí hay una versión que es un poco más eficiente en memoria (y: un tamiz adecuado, no divisiones de prueba). Básicamente, en lugar de mantener una matriz de todos los números y tachar aquellos que no son primos, esto mantiene una gran variedad de contadores, uno por cada primo que se descubra, y saltándolos por delante del presunto primo. De esa manera, utiliza el almacenamiento proporcional al número de primos, no hasta el primo más alto.

 import itertools def primes(): class counter: def __init__ (this, n): this.n, this.current, this.isVirgin = n, n*n, True # isVirgin means it's never been incremented def advancePast (this, n): # return true if the counter advanced if this.current > n: if this.isVirgin: raise StopIteration # if this is virgin, then so will be all the subsequent counters. Don't need to iterate further. return False this.current += this.n # pre: this.current == n; post: this.current > n. this.isVirgin = False # when it's gone, it's gone return True yield 1 multiples = [] for n in itertools.count(2): isPrime = True for p in (m.advancePast(n) for m in multiples): if p: isPrime = False if isPrime: yield n multiples.append (counter (n)) 

Notará que primes() es un generador, por lo que puede mantener los resultados en una lista o puede usarlos directamente. Aquí están los primeros n primos:

 import itertools for k in itertools.islice (primes(), n): print (k) 

Y, para completar, aquí hay un temporizador para medir el rendimiento:

 import time def timer (): t, k = time.process_time(), 10 for p in primes(): if p>k: print (time.process_time()-t, " to ", p, "\n") k *= 10 if k>100000: return 

En caso de que se lo pregunte, también escribí primes() como un simple iterador (usando __iter__ y __next__ ), y se ejecutó casi a la misma velocidad. ¡También me sorprendió!

Prefiero NumPy debido a la velocidad.

 import numpy as np # Find all prime numbers using Sieve of Eratosthenes def get_primes1(n): m = int(np.sqrt(n)) is_prime = np.ones(n, dtype=bool) is_prime[:2] = False # 0 and 1 are not primes for i in range(2, m): if is_prime[i] == False: continue is_prime[i*i::i] = False return np.nonzero(is_prime)[0] # Find all prime numbers using brute-force. def isprime(n): ''' Check if integer n is a prime ''' n = abs(int(n)) # n is a positive integer if n < 2: # 0 and 1 are not primes return False if n == 2: # 2 is the only even prime number return True if not n & 1: # all other even numbers are not primes return False # Range starts with 3 and only needs to go up the square root # of n for all odd numbers for x in range(3, int(n**0.5)+1, 2): if n % x == 0: return False return True # To apply a function to a numpy array, one have to vectorize the function def get_primes2(n): vectorized_isprime = np.vectorize(isprime) a = np.arange(n) return a[vectorized_isprime(a)] 

Compruebe la salida:

 n = 100 print(get_primes1(n)) print(get_primes2(n)) [ 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] [ 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] 

Compara la velocidad del Tamiz de Eratóstenes y la fuerza bruta en el Cuaderno Jupyter. Tamiz de eratóstenes en 539 veces más rápido que la fuerza bruta para millones de elementos.

 %timeit get_primes1(1000000) %timeit get_primes2(1000000) 4.79 ms ± 90.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each) 2.58 s ± 31.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) 

Pensé que debía ser posible usar simplemente la lista vacía como condición de terminación para el bucle y se me ocurrió esto:

 limit = 100 ints = list(range(2, limit)) # Will end up empty while len(ints) > 0: prime = ints[0] print prime ints.remove(prime) i = 2 multiple = prime * i while multiple <= limit: if multiple in ints: ints.remove(multiple) i += 1 multiple = prime * i 
 import math def sieve(n): primes = [True]*n primes[0] = False primes[1] = False for i in range(2,int(math.sqrt(n))+1): j = i*i while j < n: primes[j] = False j = j+i return [x for x in range(n) if primes[x] == True]