Mejorar puro Python tamiz prime por la fórmula de recurrencia

Question

Estoy tratando de optimizar aún más la solución campeón en hilo número primo mediante la adopción de la fórmula compleja para la longitud de sub-lista. len() de la misma subsecuencia es demasiado lenta ya que len es costosa y generar la subsecuencia es costoso. Esto parece acelerar un poco la función, pero aún no pude quitar la división, aunque hago la división solo dentro de la declaración de condición. Por supuesto que podría tratar de simplificar el cálculo de la longitud mediante la suscripción de la optimización de partida marcado para n en lugar de n * n ...

que sustituyen división/división entera por // para que sea compatible con Python 3 o

from __future__ import division 

También me gustaría ser interesante si esta fórmula de recurrencia podría ayudar a acelerar la solución numpy, pero no tengo experiencia en usar numpy mucho.

Si habilita psyco para el código, la historia se vuelve completamente diferente, sin embargo, y el código tamiz de Atkins se vuelve más rápido que esta técnica especial de corte.

import cProfile 

def rwh_primes1(n): 
    # http://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 primes(n): 
    # http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188 
    # recurrence formula for length by amount1 and amount2 Tony Veijalainen 2010 
    """ Returns a list of primes < n """ 
    sieve = [True] * (n//2) 
    amount1 = n-10 
    amount2 = 6 

    for i in xrange(3,int(n**0.5)+1,2): 
     if sieve[i//2]: 
      ## can you make recurrence formula for whole reciprocal? 
      sieve[i*i//2::i] = [False] * (amount1//amount2+1) 
     amount1-=4*i+4 
     amount2+=4 

    return [2] + [2*i+1 for i in xrange(1,n//2) if sieve[i]] 

numprimes=1000000 
print('Profiling') 
cProfile.Profile.bias = 4e-6 
for test in (rwh_primes1, primes): 
    cProfile.run("test(numprimes)") 

Profiling (tanta diferencia no entre las versiones)

  3 function calls in 0.191 CPU seconds 

    Ordered by: standard name 

    ncalls tottime percall cumtime percall filename:lineno(function) 
     1 0.006 0.006 0.191 0.191 <string>:1(<module>) 
     1 0.185 0.185 0.185 0.185 myprimes.py:3(rwh_primes1) 
     1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects} 


     3 function calls in 0.192 CPU seconds 

    Ordered by: standard name 

    ncalls tottime percall cumtime percall filename:lineno(function) 
     1 0.006 0.006 0.192 0.192 <string>:1(<module>) 
     1 0.186 0.186 0.186 0.186 myprimes.py:12(primes) 
     1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects} 

Curiosamente aumentando el límite a 10 ** 8 y poniendo decorador de temporización a las funciones de eliminación de la elaboración de perfiles:

rwh_primes1 took 23.670 s 
primes took 22.792 s 
primesieve took 10.850 s 

Curiosamente, si no se produce una lista de primos, sino que se devuelve el tamiz, el tiempo es aproximadamente la mitad de la versión de la lista de números.

Answer 1

que podría hacer una optimización de la rueda. Los múltiplos de 2 y 3 no son números primos, así que no los almacene en absoluto. A continuación, puede empezar a partir de 5 y salta múltiplos de 2 y 3 mediante el incremento en los pasos de 2,4,2,4,2,4 etc ..

A continuación se muestra un código C++ para ello. Espero que esto ayude.

void sieve23() 
{ 
    int lim=sqrt(MAX); 
    for(int i=5,bit1=0;i<=lim;i+=(bit1?4:2),bit1^=1) 
    { 
     if(!isComp[i/3]) 
     { 
      for(int j=i,bit2=1;;) 
      { 
       j+=(bit2?4*i:2*i); 
       bit2=!bit2; 
       if(j>=MAX)break; 
       isComp[j/3]=1; 
      } 
     } 
    } 
} 

Answer 2

Si usted puede decidir que se va a C++ para mejorar la velocidad, me portado el tamiz de Python a C++. La discusión completa se puede encontrar aquí: Porting optimized Sieve of Eratosthenes from Python to C++.

En Intel Q6600, Ubuntu 10.10, compilado con g++ -O3 y con N = 100000000 esto lleva 415 ms.

#include <vector> 
#include <boost/dynamic_bitset.hpp> 

// http://vault.embedded.com/98/9802fe2.htm - integer square root 
unsigned short isqrt(unsigned long a) { 
    unsigned long rem = 0; 
    unsigned long root = 0; 

    for (short i = 0; i < 16; i++) { 
     root <<= 1; 
     rem = ((rem << 2) + (a >> 30)); 
     a <<= 2; 
     root++; 

     if (root <= rem) { 
      rem -= root; 
      root++; 
     } else root--; 

    } 

    return static_cast<unsigned short> (root >> 1); 
} 

// https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188 
// https://stackoverflow.com/questions/5293238/porting-optimized-sieve-of-eratosthenes-from-python-to-c/5293492 
template <class T> 
void primesbelow(T N, std::vector<T> &primes) { 
    T i, j, k, sievemax, sievemaxroot; 

    sievemax = N/3; 
    if ((N % 6) == 2) sievemax++; 

    sievemaxroot = isqrt(N)/3; 

    boost::dynamic_bitset<> sieve(sievemax); 
    sieve.set(); 
    sieve[0] = 0; 

    for (i = 0; i <= sievemaxroot; i++) { 
     if (sieve[i]) { 
      k = (3*i + 1) | 1; 
      for (j = k*k/3; j < sievemax; j += 2*k) sieve[j] = 0; 
      for (j = (k*k+4*k-2*k*(i&1))/3; j < sievemax; j += 2*k) sieve[j] = 0; 
     } 
    } 

    primes.push_back(2); 
    primes.push_back(3); 

    for (i = 0; i < sievemax; i++) { 
     if (sieve[i]) primes.push_back((3*i+1)|1); 
    } 

}