HPC - Crivello di Eratostene

Moreno Marzolla moreno.marzolla@unibo.it

Ultimo aggiornamento: 2021-10-13

Il crivello di Eratostene è un algoritmo che calcola i numeri primi appartenenti ad un dato intervallo che normalmente è l'insieme \(\{2, \ldots, n\}\) . Un intero \(p \geq 2\) è primo se e solo se gli unici suoi divisori sono 1 e \(p\) (2 è un numero primo).

Per illustrare il funzionamento del crivello di Eratostene consideriamo il caso \(n=20\). Iniziamo tabulando gli interi \(2, \ldots n\):

Il primo valore della tabella (2) è un numero primo; marchiamo tutti i suoi multipli:

Il successivo valore non marcato (3) è anch'esso un numero primo. Marchiamo tutti i suoi multipli, partendo da $3 3£ (infatti \(3 \times 2\) è già stato considerato tra i multipli di 2), ottenendo:

Il successivo valore non marcato (5) è primo. Il più piccolo multiplo di 5 non ancora marcato è \(5 \times 5\), dato che \(5 \times 2\) e \(5 \times 4\) sono già stati considerati tra i multipli di 2, e \(5 \times 3\) tra i multipli di 3. Poiché \(5 \times 5 > 20\), il procedimento termina e tutti i valori non marcati sono primi:

Il file omp-sieve.c contiene un programma seriale che, dato un intero \(n \geq 1\), calcola il numero \(\pi(n)\) di primi appartenenti all'insieme \(\{2, \ldots n\}\) usando il crivello di Eratostene1. Il programma potrebbe essere reso più efficiente, ma ho deciso di favorire la comprensibilità a scapito dell'efficienza. La tabella è rappresentata dall'array isprime[] di lunghezza \(n+1\); durante l'esecuzione, isprime[k] vale 0 se e solo se k è stato marcato, cioè \(k\) è un numero composto (\(2 \leq k \leq n\)); isprime[0] e isprime[1] non sono utilizzati.

Viene fornita una funzione int mark(char *isprime, int from, int to, int p) che marca tutti i multipli di \(p\) appartenenti all'insieme \(\{\texttt{from}, \ldots \texttt{to}-1\}\). La funziona restituisce il numero di valori che sono stati marcati per la prima volta. Si richiede di parallelizzare il corpo della funzione mark() sfruttando le indicazioni che seguono.

Il corpo principale della funzione è costituito dalle istruzioni seguenti:

count = n - 1;
for (i=2; i*i <= n; i++) {
    if (isprime[i]) {
        count -= mark(isprime, i*i, n+1, i);
    }
}

Per calcolare \(\pi(n)\) si inizia ponendo count al numero di elementi dell'insieme \(\{2, \ldots n\}\); ogni volta che si marca un valore per la prima volta, si decrementa count in modo che alla fine si abbia \(\pi(n) = \texttt{count}\).

Si noti che non è possibile parallelizzare il ciclo for del frammento di codice sopra, perché il contenuto dell'array isprime[] viene modificato dalla funzione mark(), e questo causa una loop-carried dependency. Si può invece parallelizzare il corpo della funzione mark(), ad esempio dividendo l'intervallo \([i \times i, n]\) tra i thread in modo che ciascuno marchi solo i multipli di \(i\) nel proprio sottointervallo. A tale scopo si inizi usando un costrutto omp parallel determinando manualmente gli estremi dei cicli eseguiti dai singoli thread (attenzione: è complicato farlo nel modo corretto).

Compilare con:

    gcc -std=c99 -Wall -Wpedantic -fopenmp omp-sieve.c -o omp-sieve

Eseguire con:

    ./omp-sieve [n]

Esempio:

    OMP_NUM_THREADS=2 ./omp-sieve 1000

Sono riportati nella tabella seguente i valori di \(\pi(n)\) per alcuni valori di \(n\)

\(n\) \(\pi(n)\)
1 0
10 4
100 25
1000 168
10000 1229
100000 9592
1000000 78498
10000000 664579
100000000 5761455
1000000000 50847534
10000000000 Non provare sul server alloca ~10GB di RAM!!

File


  1. La funzione \(\pi(n)\) è chiamata prime-counting function