2009-06-21 14 views
21

Al momento sto facendo un progetto e ho bisogno di un metodo efficiente per calcolare i numeri primi. Ho usato il sieve of Eratosthenes ma, ho cercato in giro e ho trovato che il sieve of Atkin è un metodo più efficiente. Ho trovato difficile trovare una spiegazione (che io sia stato in grado di capire!) Di questo metodo. Come funziona? Il codice di esempio (preferibilmente in C o python) sarebbe brillante.Spiegazione del setaccio Atkin

Modifica: Grazie per il vostro aiuto, l'unica cosa che ancora non capisco è a cosa si riferiscono le variabili xey nello pseudo codice. Qualcuno potrebbe per favore far luce su questo per me?

+0

Questo assomiglia molto a una domanda università ... – Spence

+17

o una domanda Project Euler –

+0

correlati: http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n -in-python – jfs

risposta

13

Il wiki page è sempre un buon punto di partenza, poiché spiega l'algoritmo per intero e fornisce uno pseudocodice commentato. (NB: C'è un sacco di dettagli, e dal momento che il sito wiki è affidabile, io non cito qui.)

Per i riferimenti nelle lingue specifiche che hai menzionato:

Speranza che aiuta.

+5

-1. il codice copiato è per setaccio di Eratostene. E fornire collegamenti è * non * una spiegazione. –

+0

'def sieveOfErat (fine):' non era un'indicazione che questo fosse il setaccio di Eratostene e non di Atkin ?! –

2

Ecco un C++ implementazione di Crivello di Atkin che potrebbero aiutarvi ...

#include <iostream> 
#include <cmath> 
#include <fstream> 
#include<stdio.h> 
#include<conio.h> 
using namespace std; 

#define limit 1000000 

int root = (int)ceil(sqrt(limit)); 
bool sieve[limit]; 
int primes[(limit/2)+1]; 

int main (int argc, char* argv[]) 
{ 
    //Create the various different variables required 
    FILE *fp=fopen("primes.txt","w"); 
    int insert = 2; 
    primes[0] = 2; 
    primes[1] = 3; 
    for (int z = 0; z < limit; z++) sieve[z] = false; //Not all compilers have false as the  default boolean value 
    for (int x = 1; x <= root; x++) 
    { 
     for (int y = 1; y <= root; y++) 
     { 
      //Main part of Sieve of Atkin 
      int n = (4*x*x)+(y*y); 
      if (n <= limit && (n % 12 == 1 || n % 12 == 5)) sieve[n] ^= true; 
      n = (3*x*x)+(y*y); 
      if (n <= limit && n % 12 == 7) sieve[n] ^= true; 
      n = (3*x*x)-(y*y); 
      if (x > y && n <= limit && n % 12 == 11) sieve[n] ^= true; 
     } 
    } 
     //Mark all multiples of squares as non-prime 
    for (int r = 5; r <= root; r++) if (sieve[r]) for (int i = r*r; i < limit; i += r*r) sieve[i] = false; 
    //Add into prime array 
    for (int a = 5; a < limit; a++) 
    { 
      if (sieve[a]) 
      { 
        primes[insert] = a; 
        insert++; 
      } 
    } 
    //The following code just writes the array to a file 
    for(int i=0;i<1000;i++){ 
      fprintf(fp,"%d",primes[i]); 
      fprintf(fp,", "); 
    } 

    return 0; 
} 
+0

Sembra una traduzione dello pseudocodice di WP, ma guardate la pagina di discussione in cui è in discussione che lo pseudocodice non va bene. Vedi questa sottosezione http://en.wikipedia.org/wiki/Talk:Sieve_of_Atkin#How_is_this_faster_than_the_Sieve_of_Eratosthenes.3F esp. verso la sua fine –

+3

Vale la pena citare una frase di "EJ" = Emil Jeřábek come "la velocità relativa del setaccio di Atkin contro il setaccio di Eratostene è estremamente sensibile a vari trucchi di ottimizzazione utilizzati nell'implementazione" dove spiega che mentre il SoA è un po 'più teoricamente più veloce del SoE, questo è molto difficile da realizzare nella pratica. Un segmento di codice breve come pubblicato qui non avrà quelle ottimizzazioni massime, che sono molto più complesse delle ottimizzazioni di base del SoE, quindi in genere SoA è più lento del SoE per implementazioni rapide e facili. – GordonBGood

3

Edit: l'unica cosa che io continuo a non capire è ciò che il variabili x e y si riferiscono al lo pseudo codice. Qualcuno potrebbe per favore far luce su questo per me?

Per una spiegazione di base dell'uso comune dei 'x' e variabili 'Y' nella pagina di pseudo-codice Wikipedia (o meglio implementazioni del crivello di Atkin), si potrebbe trovare my answer to another related question utile.

1
// Title : Seive of Atkin (Prime number Generator) 

#include <iostream> 
#include <cmath> 
#include <vector> 

using namespace std; 

int main() 
{ 
    ios_base::sync_with_stdio(false); 
    long long int n; 
    cout<<"Enter the value of n : "; 
    cin>>n; 
    vector<bool> is_prime(n+1); 
    for(long long int i = 5; i <= n; i++) 
    { 
     is_prime[i] = false; 
    } 
    long long int lim = ceil(sqrt(n)); 

    for(long long int x = 1; x <= lim; x++) 
    { 
     for(long long int y = 1; y <= lim; y++) 
     { 
      long long int num = (4*x*x+y*y); 
      if(num <= n && (num % 12 == 1 || num%12 == 5)) 
      { 
       is_prime[num] = true; 
      } 

      num = (3*x*x + y*y); 
      if(num <= n && (num % 12 == 7)) 
      { 
       is_prime[num] = true; 
      } 

      if(x > y) 
      { 
       num = (3*x*x - y*y); 
       if(num <= n && (num % 12 == 11)) 
       { 
        is_prime[num] = true; 
       } 
      } 
     } 
    } 
    // Eliminating the composite by seiveing 
    for(long long int i = 5; i <= lim; i++) 
    { 
     if(is_prime[i]) 
      for(long long int j = i*i; j <= n; j += i) 
      { 
       is_prime[j] = false; 
      } 
    } 
    // Here we will start printing of prime number 
    if(n > 2) 
    { 
     cout<<"2\t"<<"3\t"; 
    } 
    for(long long int i = 5; i <= n; i++) 
    { 
     if(is_prime[i]) 
     { 
      cout<<i<<"\t"; 
     } 
    } 
    cout<<"\n"; 
    return 0; 
} 
9

Wikipedia article ha una spiegazione:

  • "L'algoritmo ignora completamente qualsiasi numeri divisibili per due, tre o cinque Tutti i numeri con un modulo-sessanta resto anche sono divisibile per due e non. primo numero Tutti i numeri con resto del modulo-sixty divisibile per tre sono anche divisibili per tre e non primi Tutti i numeri con resto del modulo-sixty divisibile per cinque sono divisibili per cinque e non primi Tutti questi restanti vengono ignorati. "

Cominciamo con il famoso

primes = sieve [2..] = sieve (2:[3..]) 
    -- primes are sieve of list of 2,3,4... , i.e. 2 prepended to 3,4,5... 
sieve (x:xs) = x : sieve [y | y <- xs, rem y x /= 0] -- set notation 
    -- sieve of list of (x prepended to xs) is x prepended to the sieve of 
    --     list of `y`s where y is drawn from xs and y % x /= 0 

Vediamo come si procede per un paio di primi passi:

primes = sieve [2..] = sieve (2:[3..]) 
        = 2 : sieve p2  -- list starting w/ 2, the rest is (sieve p2) 
p2 = [y | y <- [3..], rem y 2 /= 0]  -- for y from 3 step 1: if y%2 /= 0: yield y 

p2 è quello di contenere non multipli di .IOW conterrà solo i numeri coprimi con — nessun numero in p2 ha come il suo fattore. Per trovare p2 non abbiamo effettivamente bisogno di testare dividere per 2 ogni numero in [3..] (che è e fino, 3,4,5,6,7, ...), perché possiamo enumerare tutti i multipli di di anticipo:

rem y 2 /= 0 === not (ordElem y [2,4..])  -- "y is not one of 2,4,6,8,10,..." 

ordElem è come elem (cioè è elemento test), semplicemente assume che la sua tesi elenco è un'ordinata crescente lista di numeri, quindi in grado di rilevare la non presenza s afely così come la presenza:

ordElem y xs = take 1 (dropWhile (< y) xs) == [y] -- = elem y (takeWhile (<= y) xs) 

L'Ordinario elem non si assume nulla, quindi deve ispezionare ogni elemento della sua tesi lista, quindi non in grado di gestire infinite liste. ordElem can. Così, dunque,

p2 = [y | y <- [3..], not (ordElem y [2,4..])] -- abstract this as a function `diff`: 
    = diff [3..] [2,4..]    -- diff ys xs = [y | y <- ys, not (ordElem y xs)] 
    -- 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 
    -- . 4 . 6 . 8 . 10 . 12 . 14 . 16 . 18 . 20 . 22 . 
    = diff [3..] (map (2*) [2..])   -- y > 2, so [4,6..] is enough 
    = diff [2*k+x | k <- [0..], x <- [3,4]] -- "for k from 0 step 1: for x in [3,4]: 
      [2*k+x | k <- [0..], x <- [ 4]] --       yield (2*k+x)" 
    = [  2*k+x | k <- [0..], x <- [3 ]]     -- 2 = 1*2 = 2*1 
    = [3,5..]            -- that's 3,5,7,9,11,... 

p2 è solo un elenco di tutte le avversità sopra . Bene, sapevamo che era. E il prossimo passo?

sieve p2 = sieve [3,5..] = sieve (3:[5,7..]) 
         = 3 : sieve p3 
p3 = [y | y <- [5,7..], rem y 3 /= 0] 
    = [y | y <- [5,7..], not (ordElem y [3,6..])]   -- 3,6,9,12,... 
    = diff [5,7..] [6,9..]   -- but, we've already removed the multiples of 2, (!) 
    -- 5 . 7 . 9 . 11 . 13 . 15 . 17 . 19 . 21 . 23 . 25 . 27 . 
    -- . 6 . . 9 . . 12 . . 15 . . 18 . . 21 . . 24 . . 27 . 
    = diff [5,7..] (map (3*) [3,5..])      -- so, [9,15..] is enough 
    = diff [2*k+x | k <- [0..], x <- [5]] (map (3*) 
      [2*k+x | k <- [0..], x <- [3]]) 
    = diff [6*k+x | k <- [0..], x <- [5,7,9]]    -- 6 = 2*3 = 3*2 
      [6*k+x | k <- [0..], x <- [ 9]] 
    = [  6*k+x | k <- [0..], x <- [5,7 ]]    -- 5,7,11,13,17,19,... 

E il prossimo,

sieve p3 = sieve (5 : [6*k+x | k <- [0..], x <- [7,11]]) 
     = 5 : sieve p5 
p5 = [y | y <-  [6*k+x | k <- [0..], x <- [7,11]], rem y 5 /= 0] 
    = diff [ 6*k+x | k <- [0..], x <- [7,11]] (map (5*) 
      [ 6*k+x | k <- [0..], x <- [5, 7]] )     -- no mults of 2 or 3! 
    = diff [30*k+x | k <- [0..], x <- [7,11,13,17,19,23,25,29,31,35]] -- 30 = 6*5 = 5*6 
      [30*k+x | k <- [0..], x <- [     25,  35]] 
    = [  30*k+x | k <- [0..], x <- [7,11,13,17,19,23, 29,31 ]] 

Questa è la sequenza con la quale il setaccio di Atkin sta lavorando. Non ci sono multipli di 2, 3 o . Atkin e Bernstein lavorano modulo , cioè raddoppiano la gamma:

p60 = [ 60*k+x | k <- [0..], x <- [1, 7,11,13,17,19,23,29,31, 37,41,43,47,49,53,59]] 

Avanti, usano alcuni teoremi sofisticati per conoscere alcune proprietà di ciascuno di quei numeri e affrontare ogni conseguenza. Il tutto viene ripetuto (a-la "ruota") con un periodo di .

  • "All numeri n con modulo-sessanta resto 1, 13, 17, 29, 37, 41, 49, o 53 (...) sono primo se e solo se il numero di soluzioni da 4x^2 + y^2 = n è dispari e il numero è libero da quadrati. "

Che cosa significa? Se sappiamo che il numero di soluzioni a tale equazione è dispari per un numero di questo tipo, allora è primo se non è quadrato. Ciò significa che non ha fattori ripetuti (come 49, 121, ecc.).

Atkin e Bernstein utilizzano per ridurre il numero di operazioni generale: per ogni numero primo (da e superiore), annoveriamo (e contrassegno per la rimozione) multipli di suo quadrato, quindi sono molto più a parte che nel setaccio di Eratostene, quindi ce ne sono meno in un dato intervallo.

Il resto delle regole sono:

  • "Tutti i numeri n con modulo-sessanta restante 7, 19, 31, o 43 (...) sono primi se e solo se il numero di le soluzioni a 3x^2 + y^2 = n sono dispari e il numero è privo di quadrati ".

  • "Tutti n numeri con modulo-sessanta restante 11, 23, 47, o 59 (...) sono primo se e solo se il numero di soluzioni da 3x^2 − y^2 = n è dispari e il numero è di quadrati."

Questo si prende cura di tutti i numeri di base 16 nella definizione di p60.

leggi anche: Which is the fastest algorithm to find prime numbers?


La complessità temporale del crivello di Eratostene in numeri primi produrre fino a N è O (N log log N), mentre quella del crivello di Atkin è O (N) (mettendo da parte le ottimizzazioni aggiuntive log log N che non si adattano bene). La saggezza accettata è però che i fattori costanti nel setaccio di Atkin sono molto più alti e quindi potrebbe pagare solo fuori in alto sopra i numeri a 32 bit (N> 2), if at all.