[Python] numeri primi

[Enrico Franchi]

Marco Beri wrote:
> Ma ce le fai vedere o no queste 5 e poi 15 righe? O sono scienza
> segreta?  :-)

def erat_basic(buffer, end):
     limit = int(math.ceil(math.sqrt(end)))
     buffer[0] = buffer[1] = 0
     for i in xrange(2, limit):
         if buffer[i]:
             buffer[np.arange(2*i, end, step=i)] = 0

Buffer glielo passo da fuori ed e' una robba tipo:

np.ones(size+1, dtype=np.bool)


La versione da 15 fa un briciolo di unrolling a mano:

def erat(buffer, end):
     limit = int(math.ceil(math.sqrt(end)))
     buffer[0] = buffer[1] = 0
     buffer[np.arange(4, end, step=2)] = 0
     buffer[np.arange(9, end, step=3)] = 0
     buffer[np.arange(25, end, step=5)] = 0
     buffer[np.arange(49, end, step=7)] = 0
     buffer[np.arange(121, end, step=11)] = 0
     buffer[np.arange(169, end, step=13)] = 0
     buffer[np.arange(289, end, step=17)] = 0
     buffer[np.arange(361, end, step=19)] = 0
     for i in xrange(23, limit):
         if buffer[i]:
             buffer[np.arange(2*i, end, step=i)] = 0



Siccome pero' a me le funzioni di piu' di 10 righe mi infastidiscono,
sto lavorando su una cosa come:

def erat(buffer, end, small_primes=[2, 3, 5, 7]):
     limit = int(math.ceil(math.sqrt(end)))
     buffer[0] = buffer[1] = 0

     for p in small_primes:
         buffer[np.arange(p*p, end, step=p)] = 0

     for i in xrange(p+1, limit):
         if buffer[i]:
             buffer[np.arange(2*i, end, step=i)] = 0


Che da risultati comparabili a patto di scegliere bene l'array di primi 
iniziale. Qualcosa come:

def small_primes(end):
     limit = int(math.ceil(math.sqrt(end)))
     buffer = np.ones(end, dtype=np.bool)
     buffer[0] = buffer[1] = 0
     for i in xrange(2, limit):
         if buffer[i]:
             buffer[np.arange(2*i, end, step=i)] = 0
     return np.nonzero(buffer)[0]

funziona.

% time python erat2.py 100000000 10000
(array([       2,        3,        5, ..., 99999959, 99999971, 99999989]),)
python erat2.py 100000000 10000  5.99s user 0.88s system 97% cpu 7.022 total

Qui per dire gli do in pasto i primi 10000 primi. Bisogna un po' capire 
quale e' un buon valore per l'array iniziale.


-- 
.
..: -enrico-