Questa pagina illustra parte della matematica e degli algoritmi utilizzati per una efficiente ricerca dei numeri primi di Mersenne. Poiché George Woltman è più un programmatore che un matematico, si preferirà offrire una serie di puntatori al riguardo senza approfondire eccessivamente gli argomenti.

1. Creazione di una lista

Nel trattare numeri di Mersenne, è semplice dimostrare che se 2^P-1 è primo, allora P deve essere primo. Così, il primo passo per la ricerca consiste nel creare una lista di esponenti primi da testare. Fortunatamente questo lavoro è già stato svolto e sono noti almeno tutti i numeri primi compresi nell'intervallo [1- 79'300'000]

2. Fattorizzazione

Il passo successivo tenderà ad eliminare almeno alcuni candidati cercando dei fattori. Esistono algoritmi molto efficienti per determinare se un numero divida o meno il numero (2^P-1). Verifichiamo,  ad esempio, se 47 divide 2²³-1. Occorre convertire l'esponente 23 in formato binario, ovvero 10111. Iniziando da 1, occorre elevare ripetutamente al quadrato, rimuovere il bit maggiore dell'esponente e, se questo è pari ad 1, moltiplicare il valore ottenuto per 2. Infine calcolare il resto della divisione per 47.

Come appare evidente, 2²³ = 1 mod 47. Sottraiamo 1 ad entrambi i membri. 2²³-1 = 0 mod 47. Poiché abbiamo mostrato che 47 è un fattore, 2²³-1 non è primo.

Una proprietà molto elegante dei numeri di Mersenne è che qualsiasi fattore q di 2^P-1 dev'essere della forma 2kp+1. Inoltre, q deve essere 1 o 7 mod 8. Una prova può trovarsi all'indirizzo indicato. Infine, un programma efficiente può avvantaggiarsi del fatto che ogni fattore q potenziale dev'essere primo.

Il programma di fattorizzazione implementato in Prime95 crea un crivello di Eratostene modificato, in cui ciascun bit rappresenta un potenziale 2kp+1. Il crivello elimina quindi tutti quei potenziali fattori divisibili per un numero primo inferiore a 40.000 circa. Inoltre vengono azzerati i bit che rappresentano potenziali fattori di 3 o 5 modulo 8. Questo procedimento elimina approssimativamente il 95% dei potenziali fattori; i rimanenti verranno elaborati attraverso l'algoritmo di elevamento a potenza mostrato sopra.

Ora, l'unico problema da esaminare è quale livello di fattorizzazione raggiungere. La risposta dipende da tre variabili: il costo di fattorizzazione (in termini di tempo), la probabilità di trovare un fattore ed il costo (sempre in termini di tempo) di un test di primalità di Lucas-Lehmer. La formula utilizzata è dunque la seguente:

	costo_fattorizzazione < Probabilità_di_trovare_un_fattore * 2 * costo_test_primalità
 

Ciò significa che il tempo utilizzato nella fattorizzazione dev'essere inferiore del tempo risparmiato previsto. Se viene trovato un fattore, potremo evitare di eseguire tanto il test di primalità quanto il test di "double-check".
Esaminando i dati di fattorizzazione già in nostro possesso vediamo che la probabilità di trovare un fattore tra 2^X e 2^X+1 è circa 1/X. I costi dei test di fattorizzazione e di primalità vengono calcolati cronometrando il programma. Attualmente il programma ricerca fattori tra questi limiti:

3. Fattorizzazione P-1

Esiste un altro metodo di fattorizzazione utilizzato dal programma alla base del progetto G.I.M.P.S. per cercare fattori ed evitare così dispendiosi test di primalità. Questo metodo è chiamato metodo di Pollard (P-1). Se q è fattore di un numero, allora il metodo P-1 troverà il fattore q se q-1 è altamente composito - vale a dire, non ha che fattori piccoli.

Questo metodo risulta ancora più efficace se applicato a numeri di Mersenne. Come si ricorderà, il fattore q è della forma 2kp+1. Risulta quindi semplice modificare il metodo P-1 in modo tale da trovare il fattore q ogniqualvolta k sia altamente composito.

Il metodo P-1 è decisamente semplice: per prima cosa si sceglie un limite B1. P-1 cercherà il fattore q sin quando tutti i fattori di k siano inferiori a B1 (k viene chiamato B1-smooth). Si calcola quindi E, il prodotto di tutti i numeri primi inferiori a B1. Infine si calcola x = 3^E*2*P. Come ultimo passaggio si controlla il M.C.D. tra x-1 e 2^P-1 per controllare se ci sia un fattore.

Esiste un miglioramento dell'algoritmo di Pollard, che introduce l'utilizzo un secondo limite, B2, e che troverà il fattore q se k ha un solo fattore tra B1 e B2 e tutti i rimanenti fattori sono inferiori a B1.

Il progetto G.I.M.P.S. ha usato questo metodo per trovare alcuni fattori decisamente interessanti, come ad esempio:

 22944999-1 è divisibile per 314584703073057080643101377.
 
 314584703073057080643101377 = 2 * 53409984701702289312 * 2944999 + 1.
 
 K = 53409984701702289312 è altamente composito:
 
 53409984701702289312 = 25 * 3 * 19 * 947 * 7187 * 62297 * 69061
 

Come è possibile scegliere B1 e B2 in modo ottimale? E' sufficiente utilizzare una modifica della formula impiegata nella fattorizzazione. Occorre infatti massimizzare la formula:

	probabilità_di_trovare_un_fattore * 2 * costo_test_primalità - costo_fattorizzazione
 

La probabilità di trovare un fattore ed il costo di fattorizzazione variano entrambi al variare di B1 e B2. La funzione di Dickman (vedere Knuth, Art of Computer Programming, vol. 2) viene utilizzata per determinare la probabilità di trovare un fattore k sia B1-smooth o B1-smooth con un solo fattore tra B1 e B2. Il programma ricerca diversi valori di B1 e, se la memoria disponibile è sufficiente, diversi valori di B2, selezionando quelle coppie di B1 e B2 che massimizzino la formula riportata sopra.

4. Test di Lucas-Lehmer

Il test di primalità di Lucas-Lehmer è decisamente semplice: dichiara che per P > 2, 2^P-1 è primo se e solo se S_p-2 è zero in questa sequenza: S₀ = 4, S_N = S_N-1² - 2 mod 2^P-1. Ad esempio, per provare che 2⁷ - 1 sia primo:

	S0 = 4

	S1 = 4 * 4 - 2 mod 127 = 14

	S2 = 14 * 14 - 2 mod 127 = 67

	S3 = 67 * 67 - 2 mod 127 = 42

	S4 = 42 * 42 - 2 mod 127 = 111

	S5 = 111 * 111 - 2 mod 127 = 0
 

Per implementare in modo efficiente il test di Lucas-Lehmer occorre trovare il modo più rapido per elevare al quadrato grandi numeri modulo 2^P-1. Dalla fine degli anni Sessanta, l'algoritmo più rapido per elevare al quadrato grandi numeri consiste nel dividere il numero grande in pezzi più piccoli che andranno a popolare un grande vettore, quindi effettuare una Fast Fourier Transform (FFT, trasformate di Fourier veloci), un'elevazione al quadrato ed una Inverse Fast Fourier Transform (IFFT). E' possibile avere altre informazioni nella sezione "How Fast Can We Multiply?" ("Quanto rapidamente possiamo moltiplicare?") in Knuth, Art of Computer Programming vol. 2. Nell'articolo di  Richard Crandall e Barry Fagin del gennaio 1994, apparso sulla rivista Mathematics of Computation ed intitolato "Discrete Weighted Transforms and Large-Integer Arithmetic", viene introdotto il nuovo concetto "irrational base FFT". Tale miglioria ha più che raddoppiato la velocità dell'elevamento a potenza consentendo l'utilizzo di una FFT più piccola eseguendo inoltre il passaggio mod 2^P-1 in modo trasparente e senza ulteriore carico di CPU. Sebbene il progetto G.I.M.P.S. utilizzi una FFT tradizionale, ovvero in virgola mobile, per ragioni intrinseche all'architettura del processore Intel Pentium®, Peter Montgomery ha dimostrato che una trasformata all-integer weighted offre gli stessi risultati.

Come accennato nel paragrafo precedente, il progetto G.I.M.P.S. utilizza FFTs in virgola mobile scritte in linguaggio assembly che sfrutta intensivamente la cache interna e l'architettura pipelined dei moderni processori. Dal momento che i calcoli in virgola mobile tendono ad accumulare un margine di inesattezza, dopo ciascuna iterazione i valori in virgola mobile vengono arrotondati nuovamente ad interi. La discrepanza tra il risultato intero esatto e quello ottenuto tramite i calcoli in virgola mobile viene chiamato errore di convoluzione. Se l'errore di convoluzione dovesse eccedere 0,5 il passaggio di arrotondamento produrrà risultati non corretti, a dimostrazione del fatto che si doveva utilizzare una FFT più ampia. Uno degli algoritmi di correzione degli errori utilizzato dal progetto G.I.M.P.S. consiste nel garantire che il massimo errore di convoluzione non superi 0,4. Sfortunatamente questo sistema di correzione è particolarmente dispendioso in termini di tempo e di risorse e non viene eseguito ad ogni elevazione al quadrato ma solo dopo un certo numero di iterazioni. Esiste invece un altro tipo di controllo dell'errore che risulta piuttosto semplice da implementare. Una delle proprietà della moltiplicazione con FFT dichiara che:

	(somma dei valori FFT in input)2 = (Somma dei valori IFFT in output)
 

Dal momento che stiamo utilizzando numeri in virgola mobile dobbiamo cambiare il simbolo di "uguale" con quello di "approssimativamente uguale", per via dell'errore di convoluzione, appunto. Se i due valori differiscono di una quantità sostanziale, il programma produrrà un errore di tipo SUMINPUT!= SUMOUT (come descritto nel file readme.txt.). Se la somma dei valori in input alla FFT corrisponde ad un valore in virgola mobile illegale, come ad esempio un infinito, abbiamo invece un errore di tipo ILLEGAL SUMOUT. Sfortunatamente, questi controlli di errore non riescono a bloccare tutti gli errori, il che ci porta alla sezione seguente.

Quali sono le probabilità che un test di Lucas-Lehmer riconosca un nuovo numero primo di Mersenne? Un semplice approccio consiste nell'applicare ripetutamente l'osservazione secondo la quale la probabilità di trovare un fattore tra 2^X e 2^X+1 è circa 1/X. Ipotizziamo ad esempio di stare testando 2^10000139-1 che il lavoro di fattorizzazione ha dimostrato non possedere fattori inferiori a 2⁶⁴. La probabilità che esso sia primo è la (probabilità di nessun fattore a 65 bit) * (probabilità di nessun fattore a 66 bit) * ... * (probabilità di nessun fattore a 5000070 bit). Vale a dire:

	64   65         5000069

	-- * -- * ... * -------

	65   66         5000070
 

Che, semplificando, dà 64/5000070 o 1 su 78126. Questo semplice approccio non è del tutto corretto: occorrerebbe la formula di Quanto_Profondamente_Fattorizzato diviso per (esponente diviso 2). Comunque un lavoro più rigoroso ha mostrato che la formula corrisponde a (Profondità _Di_Fattorizzazione-1) / (esponente * costante di Eulero (0.577...)). In questo caso, 1 su 91623. Tuttavia anche queste formule più rigorose non sono state definitivamente provate.

5. Test di verifica (Double-Checking)

Per verificare che il test di primalità di Lucas-Lehmer sia stato effettuato senza errori, il progetto G.I.M.P.S. prevede l'esecuzione del controllo di primalità una seconda volta. Durante ciascun test vengono stampati i 64 bit di ordine inferiore del valore finale S_P-2, chiamato residuo. Se questi valori coincidono, il progetto G.I.M.P.S. può dichiarare che l'esponente è stato correttamente verificato. Se invece tali residui non coincidono, allora l'intero test di primalità viene eseguito di nuovo sin quando venga raggiunta l'equivalenza dei due risultati.

G.I.M.P.S. va in realtà oltre per garantirsi da errori di programmazione. Prima di iniziare il test di Lucas-Lehmer, sul valore S₀ viene eseguita una operazione di shift a sinistra di un numero casuale di posizioni. Ciascun elevamento al quadrato raddoppia sostanzialmente il valore delle posizioni di spiazzamento di S. Notare che il passo relativo al mod 2^P-1 ruota semplicemente il p-esimo bit ed i superiori verso i bit meno significativi, in modo da non perdere informazioni. Ma per quale motivo è necessario intraprendere questo lavoro? Se ci fosse un errore nel codice relativo alla FFT, allora lo slittamento dei valori di S assicura che le FFT nel primo test di primalità siano relative a dati completamente diversi rispetto a quelle relative al test di verifica. Sarebbe praticamente impossibile, per un errore di programma, riprodurre la medesima sequenza come residuo a 64 bit.

Storicamente, il tasso di errore per un test di Lucas-Lehmer nel quale non si siano verificati seri errori durante il calcolo è di circa 1.5%. Per i testi di Lucas-Lehmer nei quali si fosse riscontrato un errore, il tasso di errore si aggira intorno al 50%. Ai fini statistici l'errore di tipo "ILLEGAL SUMOUT" non viene contato come errore grave in quanto isolabile e corregibile dal codice stesso durante l'elaborazione.