[C++] Crivello di Eratostene: che problema ha?

[giuscri]

Stamattina per caso ho trovato sulla pagina di Wiki un algoritmo chiamato "crivello di Eratostone": è un metodo "non molto efficace" (come dice l'Enciclopedia) per trovare tutti i numeri primi presenti nell'intervallo $[a,b] \in NN$.

Descrizione:

Il procedimento è il seguente: si scrivono tutti i naturali a partire da 2 fino n in un elenco detto setaccio (in programmazione spesso l'elenco è implementato da un array). Poi si cancellano (setacciano) tutti i multipli del primo numero del setaccio (escluso lui stesso). Si prosegue così fino ad arrivare in fondo. I numeri che restano sono i numeri primi minori od uguali a n.

È come se si utilizzassero dei setacci a maglie via via più larghe: il primo lascia passare solo i numeri non multipli di 2, il secondo solo i non multipli di 3, e così via.

Mi sembra abbastanza semplice da implementare in C++. Il programma che segue compila, ma non funziona:

#include <iostream>

using namespace std;

int main(){

  const int N = 29; //Voglio cercarmi i numeri primi nell'intervallo [2, 30];

  int v[N];

  // Inizializzo il vettore:

  for(int j = 0, number = 2; j < N; j++, number++)
    v[j] = number;

  // Inizializzazione completata.

  //Crivello:

  // Durante l'esecuzione del crivello
  // i numeri multipli di altri valori già letti
  // vengono "azzerati".

  for(int j = 0; j < N; j++)

// Se v[j] == 0 allora v[j] era la posizione di un numero non primo;
// non mi interessa:

    if(v[j] != 0){


// Fissato il numero presente al posto j, vado a scorrere l'array
// cercando tutti i suoi multipli:

      for(int k = j; k < N; k++)
        if(v[k] != 0){

// Appena trovato un multiplo (che non sia zero) annullo il valore
// in quella posizione:

          if(v[k] % v[j] == 0)
            v[k] = 0;

        }

    }

  // Stampa primi:

  // stampa tutti gli elemtenti di v, a patto che siano diversi da zero:

  // ...

  return 0;

}

A parte l'annidamento di così tanti $"if"$, mi sembra che logicamente non ci dovrebbe essere alcun problema. Dov'è l'errore che faccio?

[claudio862]

Il ciclo più interno può impostare v[j] = 0 (quando j == k), e quando tenti di calcolare il resto di una divisione per zero ottieni un errore.
Ti basta far partire il ciclo più interno da k = j + 1 (anche perché se partisse da zero alla fine avresti tutti gli elementi di v nulli.

[giuscri]

"claudio86":Il ciclo più interno può impostare v[j] = 0 (quando j == k), e quando tenti di calcolare il resto di una divisione per zero ottieni un errore.
Ti basta far partire il ciclo più interno da k = j + 1 (anche perché se partisse da zero alla fine avresti tutti gli elementi di v nulli.

Come dici tu funziona.
Ma non capisco perché. Scrivendo il codice ho evitato che nella divisione $(v[k]) / (v[j])$ entrambi fossero nulli.

Il codice relativo alla ricerca dei multipli di $v[j]$ parte solo se questo non è nullo; d'altro canto anche $v[k]$ viene pescato soltanto se non è nullo. Quando arriva l'istruzione di dividerli dovrei essere sicuro che non ci sia lo zero come divisore..

Ci ripenso.

[claudio862]

if(v[j] != 0) {
    // Ok, qua v[j] != 0

    for(int k = j; k < N; k++) {
        // Alla prima iterazione k = j

        if(v[k] != 0) {
            // Alla prima iterazione v[k] == v[j] != 0

            if(v[k] % v[j] == 0) {
                // Alla prima iterazione il resto tra v[k] e v[j] è zero.
                // v[k] è multiplo di v[j]
                // (ovvio, dato che sono lo stesso elemento)

                // Alla seconda iterazione v[j] == 0, perché così è stato
                // impostato alla prima iterazione. Ma allora stai facendo
                // v[k] % 0, che è indefinito.

                v[k] = 0;
                // Alla prima iterazione v[j] viene impostato a 0
            }
        }
    }
}

Tu controlli che v[j] sia diverso da zero prima di entrare nel ciclo su k. Dentro il ciclo v[j] diventa zero, e poi viene usato come divisore.

[apatriarca]

Il crivello di Eratostene può essere abbastanza efficiente. Viene in effetti usato, con qualche ottimizzazione e variante, per calcolare i primi compresi in intervalli di numeri anche molto grossi. Non è però così come l'hai scritto che si implementa il crivello di Eratostene. Il principale vantaggio dell'utilizzo del crivello di Eratostene rispetto ad altri metodi è che il numero di multipli di un primo all'interno di un qualche intervallo non è poi molto alto. Molto inferiore al numero di valori in tale intervallo. Nel tuo codice però, invece di prendere in considerazione solo i multipli, stai iterando su tutti i valori. Il tuo algoritmo diventa quindi del tutto equivalente ad una ricerca esaustiva dei divisori di ogni numero nell'intervallo (in alcuni casi potrebbe anzi essere più veloce questo ultimo algoritmo).

Come si implementa quindi un crivello di Eratostene? La prima cosa da osservare è che il tuo array è prima di tutto inizializzato con interi che sono praticamente uguali all'indice corrispondente. Sono infatti uguali all'indice meno due. In seguito cambi alcuni di questi valori a zero. Ogni elemento di questo array contiene quindi solo due valori possibili; uno dei quali dipende dall'indice. Lavorare con gli indici è spesso preferibile. Possiamo infatti a questo punto considerare semplicemente un array di valori booleani che contengono true se l'indice corrispondente più due è un primo e false in caso contrario (in realtà si possono usare anche strutture dati più complesse ma visto che sei alle prime armi lasciamo perdere). Chiamiamo questo array is_prime. Per sapere quindi se il numero \(m\) è primo sarà sufficiente verificare che is_prime[ m - 2 ] sia true.

A questo punto vediamo altre ottimizzazioni classiche. Osserviamo prima di tutto che ogni numero \(m\) non primo è il prodotto di almeno due primi, uno dei quali inferiore a \(\sqrt m.\) Tutti i multipli dei primi maggiori di \( \sqrt N \) (\(N\) intero massimo sotto il quale vogliamo trovare i primi) saranno già stati quindi eliminati nei cicli di eliminazione dei primi più piccoli di tale valore. Un discorso simile vale quando devi eliminare i multipli di un primo \(p\) fissato. Tutti i multipli di \(p\) inferiori a \(p^2\) saranno già stati eliminati da qualche primo più piccolo. Possiamo quindi considerare solo i multipli di \(p\) maggiori di \(p^2.\)

Altre ottimizzazioni classiche comprendono l'eliminazione di tutti i numeri pari dall'array. In questo caso all'indice \(i\) corrisponde il numero \( 3 + 2\,i. \) In modo simile si possono eliminare anche le potenze di \(3\) considerando solo i numeri interi congruenti a \( 1 \bmod 6 \) oppure \( 5 \bmod 6 \). Ovviamente la relazione tra indici e valori diventa in questo caso molto più complicata.

P.S. Per quanto riguarda l'errore nel tuo codice, ti ha già risposto claudio86. Se non ti è chiaro la motivazione ti consiglio di usare un debugger o di provare ad eseguire il codice una riga per volta a mano.

[giuscri]

"claudio86":

if(v[j] != 0) {
    // Ok, qua v[j] != 0

    for(int k = j; k < N; k++) {
        // Alla prima iterazione k = j

        if(v[k] != 0) {
            // Alla prima iterazione v[k] == v[j] != 0

            if(v[k] % v[j] == 0) {
                // Alla prima iterazione il resto tra v[k] e v[j] è zero.
                // v[k] è multiplo di v[j]
                // (ovvio, dato che sono lo stesso elemento)

                // Alla seconda iterazione v[j] == 0, perché così è stato
                // impostato alla prima iterazione. Ma allora stai facendo
                // v[k] % 0, che è indefinito.

                v[k] = 0;
                // Alla prima iterazione v[j] viene impostato a 0
            }
        }
    }
}

Tu controlli che v[j] sia diverso da zero prima di entrare nel ciclo su k. Dentro il ciclo v[j] diventa zero, e poi viene usato come divisore.

Ah! Cavoli. Ora capisco qual è il senso di partire da $k = j+1$. Grazie mille!

[giuscri]

"apatriarca":Il principale vantaggio dell'utilizzo del crivello di Eratostene rispetto ad altri metodi è che il numero di multipli di un primo all'interno di un qualche intervallo non è poi molto alto. Molto inferiore al numero di valori in tale intervallo. [...]
Come si implementa quindi un crivello di Eratostene? La prima cosa da osservare è che il tuo array è prima di tutto inizializzato con interi che sono praticamente uguali all'indice corrispondente. Sono infatti uguali all'indice meno due. In seguito cambi alcuni di questi valori a zero. Ogni elemento di questo array contiene quindi solo due valori possibili; uno dei quali dipende dall'indice. Lavorare con gli indici è spesso preferibile. Possiamo infatti a questo punto considerare semplicemente un array di valori booleani che contengono true se l'indice corrispondente più due è un primo e false in caso contrario (in realtà si possono usare anche strutture dati più complesse ma visto che sei alle prime armi lasciamo perdere). Chiamiamo questo array is_prime. Per sapere quindi se il numero \(m\) è primo sarà sufficiente verificare che is_prime[ m - 2 ] sia true.

Bellissima questa idea! Più tardi provo a buttare giù qualcosa.

EDIT: ci ho ripensato sù questa sera e le implementazioni che mi sono venute in mente erano identiche a quelle di oggi pomeriggio -a parte usare variabili booleane al posto di interi.

Ho pensato che avessi potuto farmi un discorso di memoria occupata (le variabili booleane occupano un solo bit, vero?), ma vista l'ultima discussione che abbiamo avuto non credo fosse quello il senso del tuo consiglio.

C'è da credere che le migliorie siano nelle informazioni che mi hai dato nel resto del messaggio.

Comunque domani provo a postare qualche riga di codice perché ad essere sincero ci ho pensato approssimativamente e con leggerezza.

Post scriptum: non è la prima volta che mi dici che è preferibile usare gli indici di un array rispetto a caricarci valori sopra -quando è possibile. A me sembra equivalente, a volermi chiudere gli occhi e a non capire che usare gli indici è molto più "compatto"; ma c'è una qualche ragione più profonda, oltre alla praticità?