Algoritmi per numeri primi: risultati, benchmark, considerazioni

Giulio_M

Ho deciso di "sperimentare" qualche nuovo algoritmo per il calcolo dei numeri primi, ovvero partendo da un algoritmo classico, cercare varie strade di ottimizzazione. Lo scopo è la massima efficienza (parliamo di numeri molto grandi, quindi in totale un numero grande di operazioni da valutare...), quindi:

linguaggio C, programmazione vicina al "linguaggio nativo", variabili forzatamente memorizzate nei registri della CPU anziché nella RAM, in modo da avere maggiore velocità di esecuzione
condizioni di valutazione, "filtrare i dati" a priori, in modo da avere minori operazioni da processare dopo (quando si parla di grandi numeri, è cercare conveniente questa strada)

In particolare, ho voluto risolvere questo tipo di problema: trovare il numero primo che segue 10¹³ (spoiler: 1000000000039), che tra l'altro è definito un "numero primo curioso" (fonte) in quanto è il più piccolo numero primo di 13 cifre la cui somma delle cifre sia uguale a 13.

Premessa: ho anche cercato online suggerimenti sui vari algoritmi, alcuni come il Crivello di Erastotene apparentemente ingegnosi e "carini" come idea (più che altro si applicano al calcolo di tutti i numeri primi, dall'origine fino ad N), peccato che dal punto di vista pratico, computazionale, siano assolutamente improponibili oltre un certo valore. Ad esempio con 10¹³, una variabile memorizzata come long int (sizeof(long int) = 8 byte), richiederebbe 8 * 10¹³ byte ovvero per capirci 80 HDD da 1 TB (1000 GB ognuno); per non parlare poi della ridotta velocità di accesso a tale memoria. Quindi semplicemente improponibile.

Riporto quindi i miei algoritmi, comunque partendo da una base di ottimizzazione già buona, sia di algoritmo che computazionale, gestione variabili e memoria (codice meno ottimizzato avrebbe poco senso considerarlo).

algoritmo 1: abbastanza "classico", ovvero per ogni numero "j" potenziale numero primo, lo divido per un numero "i" che va da 2 fino a (int)(N/2). Ovviamente la parte da N/2 a N, la possiamo già escludere. Ecco il codice:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main(){
register short cont=0;
register long int j=1000000000000;
j+=1;
register long int i=2;
clock_t start = clock();
while(cont!=1){
j++;
for(i=2;i<(int)(j/2);i++){
if(j%i==0){break;}
}
if(j%i!=0){cont++;}
}
clock_t end = clock();
printf("%f s\n",(float)(end-start)/CLOCKS_PER_SEC);
printf("%ld\n",j);
return 0;
}

algoritmo 2: rispetto a prima, ho pensato: perché cavolo analizzare ben N/2 operazioni (N>10¹³) per tutti i numeri pari?? Quindi si aggiunge una condizione che se "j" è pari, salta direttamente al prossimo step, in questo modo anche l'indice "i" parte da 3 e può saltare di 2 anziché di 1! Si dimezzano le operazioni da valutare. Ecco come cambia il codice:
```
if(j%2==0){continue;}
for(i=3;i<(int)(j/2);i+=2){
```
algoritmo 3: la parte for(i=3;i<(int)(j/2);i+=2){ la lasciamo così, prima invece facciamo questa considerazione: non solo "i" ma ovviamente anche "j" può saltare avanti di 2, ovvero valutare solamente i numeri dispari, quindi, sebbene poco rilevante in questo caso (si tratta della metà di 39 iterazioni... da cui esce subito verificando j%2==0)
possiamo aggiungere questa modifica e togliere la valutazione di j%2==0 (ovvero se il numero è pari, passa subito al numero successivo, il simbolo % indica appunto il resto della divisione); inoltre, per ogni numero j, come accennavo all'inizio (condizioni di valutazione, "filtrare i dati"), per questo caso ovvero non creare una lista di tutti i numeri primi fino ad N, bensì trovare il numero primo successivo ad un dato numero, può avere senso una valutazione a priori, quindi rispetto a prima ciò che cambia è questo:
```
j+=2;
if(j%3==0){continue;}else{
if(j%5==0){continue;}else{
if(j%7==0){continue;}else{
if(j%11==0){continue;}else{
if(j%13==0){continue;}else{
if(j%17==0){continue;}else{
if(j%19==0){continue;}else{
if(j%23==0){continue;}else{
if(j%29==0){continue;}
}}}}}}}}
```
Per ovvie ragioni mi sono fermato, si poteva continuare con j%31==0 e così via; diciamo che già questo "alleggerisce" il carico da valutare dopo.

Vediamo ora i risultati:

tempo di esecuzione algoritmo 1: 17,5 secondi
tempo di esecuzione algoritmo 2: 9 secondi
tempo di esecuzione algoritmo 3: 8,7 secondi

Cosa ne pensate? Idee di un algoritmo per rendere ancora più efficiente (tendenzialmente significa ridurre il numero di operazioni complessive per la CPU)? 🙂

✨Cassiopea ❣️ Singing A Song✨

Giulio_M Cosa ne pensate? I

Tutto ciò è affascinante. Che robe complesse. Wow! 😍

Matematico_risponde

Giulio_M All'esame di Informatica Applicata (C++) presi anche 28 con gli automi a stati finiti (interi algoritmi da scrivere a mano e su un foglio) dopo essermi esercitato per 4 mesi circa e aver letto un libro intero di programmazione. Peccato che ora non ricordi nulla. Ma la tentazione di imparare di nuovo non manca.

Antonio_Scalzo

È parecchio che non programmo più, tuttavia mi ricordo (quando i PC erano lentissimi) il metodo che avevo escogitato per snellire i calcoli: aggiungevo una terza variabile k=sqrt(j) tale che i <= K

Cioè, per provare se un numero intero J è primo, è sufficiente che non abbia divisori tra 2 e la sua radice quadrata, giacché tutti i divisori di J maggiori di sqrt(J) ... devono restituire numeri minori di sqrt(J) e che a loro volta sono dei divisori di J.

In particolare, adottando questo metodo, per provare se il numero mille miliardi (1000 000 000 000) è primo ... basta fare massimo un milione di iterazioni e non mille miliardi di iterazioni.

Giulio_M

Antonio_Scalzo beh, fantastico! Infatti l'ho testato e devo dire che (empiricamente) funziona alla grande, dato che se confrontiamo anche la mia ottimizzazione, N/2 >> sqrt(N) quando andiamo a prendere N sufficientemente grande.

In particolare, dal mio codice migliore (algoritmo 3) di 8,7 secondi, con questa ottimizzazione ho ottenuto 0,029 secondi! Un miglioramento decisamente ottimo 🙂

Non ho modo però di dimostrarne la validità generale (ok, possiamo dimostrarlo empiricamente per una lista più o meno grande di numeri primi...), così come affinare ulteriormente il metodo ovvero, c'è una ragione precisa (che magari mi sfugge) per cui debba arrivare alla radice quadrata del numero come valore limite? Oppure si può dimostrare la validità anche andando spingendosi oltre (quindi anche fermandosi ad un valore più piccolo di sqrt(j))?

Attualmente riuscirei solo a fare una verifica su base empirica. Intendo dire, la frase:

Antonio_Scalzo giacché tutti i divisori di J maggiori di sqrt(J) ... devono restituire numeri minori di sqrt(J) e che a loro volta sono dei divisori di J.

è sempre un'osservazione ricavata a posteriori, su base empirica, oppure una dimostrazione della validità generale? E come dicevo, in questo caso, perché proprio esattamente sqrt(j) come valore limite?

Edit: ho provato anche a spingermi leggermente oltre, ovvero anziché j^0.5 ho imposto j^0.4 ma su larga scala, andando a ricostruire la lista di numeri primi, in alcuni casi però non viene rispettata la condizione. Con sqrt(j) invece ho testato qualche migliaio di numeri primi ed effettivamente risulta verificato, coincidente con la lista calcolata in modo "tradizionale"! Al massimo posso testarlo su tutta la lista fino ad N, chiaro che risolverlo in modo empirico non è sufficiente come dimostrazione di validità generale.

Antonio_Scalzo

Giulio_M è sempre un'osservazione ricavata a posteriori, su base empirica, oppure una dimostrazione della validità generale?

Certo che è dimostrabile matematicamente!

Posta l'eguaglianza xy = qq

E ipotizzando che x > q e y > q ... avremo:

x - 1 >= q e y - 1 >= q ... e quindi: (x-1) * (y-1) >= q*q

da cui ricaviamo ... xy - x - y + 1 >= qq

Ossia xy + 1 >= qq + x + y ... Il che contraddice l'ipotesi di partenza.

zeunig356

Giulio_M E' un teorema. Se p é un divisore di N, anche q = N/p lo é.
e naturalmente se p > sqrt(N) allora q < sqrt(N). Per cui se non hai trovato divisori primi di N entro la sua radice quadrata non ne troverai più. Sarebbe interessante sapere se questo risultato può essere migliorato.

Giulio_M

Antonio_Scalzo pensandoci, mi è venuto in mente dopo, si può anche semplicemente dimostrare in questo modo:
Se x è divisore di N, il massimo valore possibile di x si può ottenere con x*x=N, da cui esclusa la soluzione negativa, abbiamo appunto x=sqrt(N).

Edit: più precisamente, la ragione di quanto ho scritto sopra, è questa:
Dobbiamo avere j(N-j) = max. Quindi d[j(N-j)]/dj=0, che fornisce come risultato N-2j=0 ovvero j=N/2.

Giulio_M

Aggiornamento - efficienza dell'algoritmo notevolmente migliorata:

Ho introdotto qualche modifica al codice, prendendo spunto da altri lavori. Dalla migliore prestazione precedente ovvero 0,029 secondi, siamo ora arrivati a 0,007175 secondi, vale a dire 4 volte più veloce del precedente e (incredibilmente) 1212 volte più veloce dell' "algoritmo 3". L'algoritmo è talmente veloce ed efficiente che non occorre nemmeno fare riferimento al calcolo parallelo (es. #pragma omp parallel { /* codice */ }), sarebbe controproducente in termini di prestazioni (impiega più tempo nella suddivisione e gestione dei vari thread rispetto all'esecuzione in sé, rapidissima come abbiamo visto!

Ecco tutto il codice con le nuove ottimizzazioni, anche più compatto e pulito, avendo fatto riferimento ad alcune librerie da includere. Il fattore chiave ovviamente è la nuova funzione definita is_prime(), che lavora con dati booleani (true-false) quindi la massima leggerezza possibile in termini di dimensioni variabili.

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <stdbool.h>

bool is_prime(unsigned long int numero) {
double radice = sqrt(numero);
if(numero%2==0){return false;}
for (register unsigned int i = 3; i <=radice; i+=2){
if (numero % i == 0) {
return false;
}
}
return true;
}

int main(){
register unsigned long int j=1000000000000;
clock_t start = clock();
while(is_prime(j)==false){
j++;
}
clock_t end = clock();
printf("%f s\n",(float)(end-start)/CLOCKS_PER_SEC);
printf("%ld\n",j);
return 0;
}

Giulio_M

Giusto un confronto, benchmark prestazioni fra il precedente codice in C (Giulio_M) e un linguaggio interpretato come Python, comodo per molti aspetti da con le evidenti criticità di efficienza (avendo comunque inserito le possibili ottimizzazioni di algoritmo come fatto per il programma scritto in C):

codice C: come scritto prima, il tempo di esecuzione è di 0,007175 secondi
codice Python: tempo di esecuzione di 0,15 secondi ovvero 21 volte più lento di C

Da qui vediamo che per problemi computazionali dove viene richiesta maggiore efficienza (principalmente per questioni di tempo di esecuzione del codice), è importante oltre all'algoritmo anche la scelta di un linguaggio compilato come il C, che consente le migliori prestazioni possibili, anche ordini di grandezza rispetto ad altri linguaggi di programmazione. Qui abbiamo tempi molto ridotti, inferiori al secondo, quando invece abbiamo a che fare con diversi minuti oppure ore... La differenza è decisamente considerevole!

Per completezza riporto il codice Python3:

import time
j=1000000000000;
i=3
def is_prime(j):
    if(j%2==0):
        return 0
    else:
        i=3
        while(i<j**0.5):
            if(j%i==0):
                return 0
            i+=2
    return 1
start=time.time()
if(j%2==0):
    j+=1
while(is_prime(j)==0):
    j+=2
end=time.time()
print(j)
print("tempo di esecuzione = ",round(end-start,2),"s")

Giulio_M

zeunig356 certo, in termini di efficienza, l'ultimo aggiornamento che ho messo (il codice iniziale infatti aveva molti margini di miglioramento, ho "scoperto" dopo questo fatto della radice quadrata e con N grande, srqt(N) vs N/2, il numero di operazioni richieste aumenta notevolmente). Infatti come avevo poi ragionato su questo:

Se x è divisore di N, il massimo valore possibile di x si può ottenere con x*x=N, da cui esclusa la soluzione negativa, abbiamo appunto x=sqrt(N).

Per diminuire il numero di operazioni totali, oltre alle varie ottimizzazioni del codice sono arrivato anche a questo:

if(numero%2==0){return false;}
for (i = 3; i <=radice; i+=2){
if (numero % i == 0) ...

Ovvero anziché i++ che scorre tutti i numeri, verifico prima se è peri, poi parto da 3 e salto di 2. In pratica dimezza le operazioni totali e quindi il tempo di esecuzione! Come migliorarlo ulteriormente non saprei (oltre al fatto che per questo caso, verifica del numero primo che segue un numero, per grande che sia come 10¹³, ho raggiunto la soluzione in 0,007175 secondi.

Anche per quanto riguarda proprio la lista di numeri primi fino ad un numero N (ad esempio con un algoritmo analogo all'ultimo che ho messo, la lista del primo milione di numeri primi l'ho trovata in 2,77 secondi tramite calcolo parallelo, con una macchina normale Dual Core 2.6 GHz 4 thread, anche qui dubito di ottimizzarlo più di così) un altro algoritmo efficiente è il Crivello di Eratostene, tuttavia nella pratica di implementazione è conveniente solo fintanto che N non è troppo grande, altrimenti diventa impraticabile (overflow di memoria): prima salva tutta la lista di numeri e poi li filtra. Con 10¹³ dovendo usare "long int" abbiamo 8byte ogni elemento, ovvero 80000 GB di dati solo per memorizzare le variabili? Improponibile ahah, almeno in questa forma "standard".

Giulio_M

Numero primo più grande di sempre, scoperto da un ex dipendente NVIDIA

Luke Durant, ragazzo di trentasei anni ex dipendente dell'azienda NVIDIA, ha scoperto (anche con un po' di fortuna!) il numero primo più grande di sempre, dopo che il precedente record (oltre 20 milioni di cifre) è durato sei anni.
Questo nuovo numero primo, scoperto ad ottobre 2024 ha oltre 41 milioni di cifre ed è pari a 2^136279841-1 (anche indicato con la sigla M136279841 rientrando nella categoria dei numeri di Mersenne). Luke Durant si è unito al progetto cooperativo GIMPS (Great Internet Mersenne Prime Search), progetto di calcolo distribuito; il numero è stato salvato in un file di testo dal peso di circa 40 MB.

Pubblicazione ufficiale su mersenne.org

La ricompensa non è certamente enorme (<<$3000 award offered to anyone lucky enough to find a new prime>>), ma di sicuro entrare nella storia in questo modo, specie per un non-professionista del settore, è una bella soddisfazione! Secondo quanto chiaratato, Luke Durant non vuole nemmeno tenere per sé questi soldi (dato che l'importo non cambia di certo la vita), ma li dona in ricerca: <<He plans to donate the award to the math department at the Alabama School of Math and Science>>.