piergiac-1
Nuovo Utente
- Messaggi
- 50
- Reazioni
- 4
- Punteggio
- 20
Non me la prendo, perchè dovrei? La tua risposta comunque la prendo come un si.
Pubblicità
C++ - Ricerca numeri primi
- Autore discussione M1n021
- Data d'inizio
Pubblicità
Non me la prendo, perchè dovrei? La tua risposta comunque la prendo come un si.
- Messaggi
- 2,345
- Reazioni
- 1,934
- Punteggio
- 134
Non si se hai spulciato anche il sorgente per vedere nextprime: https://github.com/sympy/sympy/blob/sympy-1.13.2/sympy/ntheory/generate.py#L629
Il crivello usato è quello di Eratostene, guardando sieve() (espanso dinamicamente, come fai tu).
M1n021
Utente Attivo
- Messaggi
- 279
- Reazioni
- 97
- Punteggio
- 43
No, ma cercando su google avevo letto che utilizza un algoritmo su base probabilistica per valori maggiori di 2^64-1. Forse utilizza il crivello per numeri più piccoli e il metodo probabilistico per quelli più grandi. Per caso sei riuscito a capire anche quale algoritmo probabilistico utilizza?
- Messaggi
- 2,345
- Reazioni
- 1,934
- Punteggio
- 134
Avevo letto anche io, ma guardando il codice la cosa non mi torna molto (cioè, non mi sembra di vedere controlli di quel tipo, ma non ho verificato isprime () che check per esempio). Però l'ho guardato da smartphone, da PC provo a dare un'altra occhiata.
EDIT:
"Finally if the number is larger than 2^64, a strong
BPSW test is performed. While this is a probable prime test and we
believe counterexamples exist, there are no known counterexamples"
È proprio nel check di isprime().
Ora leggo il codice
sympy/sympy/ntheory/primetest.py at sympy-1.13.2 · sympy/sympy
A computer algebra system written in pure Python. Contribute to sympy/sympy development by creating an account on GitHub.
github.com
M1n021
Utente Attivo
- Messaggi
- 279
- Reazioni
- 97
- Punteggio
- 43
Eh sì, essendo l'algoritmo probabilistico in fin dei conti un test di primalità, è giusto che sia lì.
Non conosco il python, ma dai riferimenti riportati mi sembra di capire che il test utilizzato sia il seguente:
Baillie–PSW primality test - Wikipedia
piergiac-1
Nuovo Utente
- Messaggi
- 50
- Reazioni
- 4
- Punteggio
- 20
Test Deterministici
1. Divisione per Tentativi
- Caratteristiche: Semplice ma inefficiente per numeri grandi.
- Uso: Numeri piccoli, situazioni didattiche.
- Perché: Facile da implementare e comprendere.
2. Crivello di Eratostene
- Caratteristiche: Efficiente per trovare tutti i primi fino a un certo limite.
- Uso: Generazione di tavole di numeri primi, applicazioni matematiche.
- Perché: Ottimo per pre-calcolare set di numeri primi.
3. Test AKS (Agrawal-Kayal-Saxena)
- Caratteristiche: Primo test deterministico in tempo polinomiale.
- Uso: Principalmente di interesse teorico.
- Perché: Prova la primalità in tempo polinomiale, ma spesso troppo lento per usi pratici.
4. Test di Lucas-Lehmer
- Caratteristiche: Specifico per numeri di Mersenne.
- Uso: Ricerca di grandi numeri primi di Mersenne.
- Perché: Estremamente efficiente per questa classe specifica di numeri.
5. Test Deterministico di Miller
- Caratteristiche: Versione deterministica del test di Miller-Rabin.
- Uso: Prova rigorosa di primalità.
- Perché: Garantisce la correttezza, ma richiede la conoscenza dell'ipotesi di Riemann estesa.
Test Probabilistici
6. Test di Fermat
- Caratteristiche: Rapido ma può fallire con i numeri di Carmichael.
- Uso: Come pre-test rapido prima di test più rigorosi.
- Perché: Veloce e semplice da implementare.
7. Test di Miller-Rabin
- Caratteristiche: Molto affidabile e efficiente.
- Uso: Ampiamente utilizzato in crittografia e software matematico.
- Perché: Offre un buon equilibrio tra velocità e affidabilità.
8. Test di Solovay-Strassen
- Caratteristiche: Basato sul simbolo di Jacobi.
- Uso: Meno comune del Miller-Rabin, ma ancora utilizzato in alcuni contesti.
- Perché: Storicamente importante, ancora utile in alcune applicazioni specifiche.
9. Test di Baillie-PSW
- Caratteristiche: Combina Miller-Rabin e un test di Lucas probabilistico.
- Uso: Verifica di primalità in software matematici.
- Perché: Molto affidabile, nessun falso positivo noto.
10. Test di Frobenius
- Caratteristiche: Basato sulle proprietà dell'anello degli interi modulo n.
- Uso: Alternative al Miller-Rabin in alcuni contesti.
- Perché: Può essere più veloce per alcuni numeri specifici.
Test Ibridi e Avanzati
11. ECPP (Elliptic Curve Primality Proving)
- Caratteristiche: Utilizza curve ellittiche per provare la primalità.
- Uso: Prova rigorosa per numeri molto grandi.
- Perché: Può fornire certificati di primalità per numeri estremamente grandi.
12. APR-CL (Adleman-Pomerance-Rumely, Cohen-Lenstra)
- Caratteristiche: Test ciclotomico.
- Uso: Prova rigorosa per numeri di dimensioni moderate.
- Perché: Più veloce dell'ECPP per alcuni intervalli di numeri.
13. Test di Proth
- Caratteristiche: Specifico per numeri della forma k·2^n + 1.
- Uso: Ricerca di numeri primi di Proth.
- Perché: Efficiente per questa forma specifica di numeri.
Ultima modifica da un moderatore:
- Messaggi
- 2,345
- Reazioni
- 1,934
- Punteggio
- 134
Davo per scontato fosse lì, non poteva che essere così considerando che il risultato effettivo lo restituisce quella funzione e ogni blocco del codice ne fa uso.Eh sì, essendo l'algoritmo probabilistico in fin dei conti un test di primalità, è giusto che sia lì.
Non conosco il python, ma dai riferimenti riportati mi sembra di capire che il test utilizzato sia il seguente:
Baillie–PSW primality test - Wikipedia
en.wikipedia.org
Si, comunque è quello sicuramente usato.
Grazie per averci incollato la risposta di GPT...
M1n021
Utente Attivo
- Messaggi
- 279
- Reazioni
- 97
- Punteggio
- 43
Per curiosità sono andato a vedere l'implementazione del crivello, anche per fare un paragone con la mia, ma sinceramente ad un lettura veloce, considerando anche il fatto che non conosco il python, non ci ho capito più di tanto!
Più che altro mi interessava capire come affrontava questa specifica questione:
- Messaggi
- 25,070
- Reazioni
- 12,679
- Punteggio
- 254
potresti migliorare quanto hai già fatto usando una frazione della memoria, ci stavo pensando su: invece di memorizzare interi, potresti usare le sequenze di bit (in C++ è il tipo di dato bitset); in pratica un 1/true dice che un certo intero è primo, 0/false che non lo è;
per semplificare cosa voglio dire, faccio una corrispondenza 1:1 tra vettore di bit e numeri interi
l'indice i del bitset corrispone all'intero i
il vettore va inizializzato con tutti 1 (presumi falsamente che tutti gli interi siano primi), poi si mettono a zero i primi 2 bit perché 0 e 1 non sono primi; supponendo che il bitseti si chiami p
0011111111111111111111....
dice che 0 e 1 non sono primi (perché p[0]=0, p[1]=0) mentre tutti gli altri sono primi
l'indice successivo è il 2, p[2] dice che il numero 2 è primo --> crivelli il seguito del bitset eliminando multipli di 2 (eliminare = metti il bit a 0)
il bit successivo ti dice che 3 è primo --> crivelli i multipli di 3
il bit successivo (corrisponde al numero 4) è stato crivellato dal 2 --> non primo perché p[4]=0
e così via...
alla fine hai un bitset simile a questo (i primi 2 bit sono lo 0 e l'1 che non sono primi):
001101010001010001010001000001010000010001010001000...
che corrsisponde ai numeri
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 (primi <=50)
in questo modo, per interi a 32 bit (o più: 64...), con il bitset usi un solo bit per numero, invece se usi un array di int ne usi almeno 32
questo significa che riduci l'uso della memoria di 32 volte
io l'esempio l'ho fatto usando anche i pari e includendo 0 e 1, potresti fare la stessa cosa usando il tuo approccio che usa solo i dispari ed il calcolo degli indici;
ovviamente vanno prese accorgimenti: con il primo bitset costruisci un vector "di base" (da usare in seguito) dove gli interi vanno memorizzati per esteso (basta fare pushback degli indici: nell'esempio precedente l'indice 29 è il primo 29); per le espansioni usi il bitset, e alla fine lo scorri per ricreare gli interi.
Questo metodo perde un po' in efficienza perché alla fine di ogni espansione va comunque aggiornato il vector iniziale con i nuovi primi trovati), per le espansioni successive risparmi memoria.
- Messaggi
- 2,345
- Reazioni
- 1,934
- Punteggio
- 134
M1n021
Utente Attivo
- Messaggi
- 279
- Reazioni
- 97
- Punteggio
- 43
@BAT In realtà il consumo di memoria in sé non è un gran problema nel mio programma, visto che non ogni espansione del crivello corrisponde ad un'espansione del vettore di interi, in quanto la densità dei numeri primi tende man mano a diminuire. Questo significa che il massimo spreco di memoria rispetto allo spazio necessario per immagazzinare i numeri primi sarà (ipotizzando di utilizzare un vettore di
In ogni caso sarebbe interessante vedere se cambia qualcosa dal punto di vista delle prestazioni, quindi appena ho tempo lo implemento e ti faccio sapere!
Ecco il programma aggiornato con i bitset:
- ho definito il bitset come variabile statica perché le variabili automatiche (o dinamiche che dir si voglia) sono molto più restrittive in termini di memoria a disposizione. Ciò però implica anche che il bitset vada resettato ogni volta nelle chiamate a
- nella prima chiamata a
- tenendo conto del teorema dei numeri primi ho aggiunto la seguente riga di codice
- a parità di memoria utilizzata per immagazzinare i nuovi valori da sottoporre a crivellatura nella funzione
L'ho scritto abbastanza velocemente, se riscontrate qualche errore o vi viene in mente qualche miglioria, fatemi sapere.
uint64_t ) di N*64bit, che per N=1.000.000 (così come l'ho impostato nel programma) significa "solo" 8MB.In ogni caso sarebbe interessante vedere se cambia qualcosa dal punto di vista delle prestazioni, quindi appena ho tempo lo implemento e ti faccio sapere!
Grazie! Più che altro mi interessa capire come affronta la specifica questione che ho riportato nel mio precedente post.
Ecco il programma aggiornato con i bitset:
C++:
#include <iostream>
#include <bitset>
#include <cstdint>
#include <cmath>
#include <vector>
using namespace std;
void espandi_e_applica_crivello(vector<uint32_t> &v, uint32_t &last)
{
static const uint32_t N = 50'000'000;
static bitset<N> b;
uint32_t first = last + 2;
last += 2 * N;
if(v.size() == 1)
{
for(uint32_t i = 0, n, q; ; ++i)
{
if(!b[i])
{
n = 2 * i + first;
if((q = n * n) > last)
{
break;
}
for(uint32_t j = (q - first) / 2; j < N; b.set(j), j += n);
}
}
}
else
{
b.reset();
for(uint32_t i = 1, j, n, q; (q = (n = v[i]) * v[i]) <= last ; ++i)
{
if(q < first)
{
uint32_t r = first % n;
j = r ? (r % 2 ? (n - r) / 2 : n - r / 2) : 0;
}
else
{
j = (q - first) / 2;
}
for(; j < N; b.set(j), j += n);
}
}
for(uint32_t i = 0; i < N; ++i)
{
if(!b[i])
{
v.push_back(2 * i + first);
}
}
}
uint32_t next_prime_number(const uint32_t n)
{
static vector<uint32_t> v = {2};
//v.reserve(n / log(n));
static uint32_t last = 1;
while(v.back() <= n)
{
espandi_e_applica_crivello(v, last);
}
cout << "\n" << v.size() << " NUMERI PRIMI <= " << last << " (" << (double)v.size() / last * 100 << "%)\n";
uint32_t i_dx = v.size() - 1;
for(uint32_t i_sx = 0, i; i_dx != i_sx + 1; (n < v[i = (i_sx + i_dx) / 2] ? i_dx : i_sx) = i);
return v[i_dx];
}
int main()
{
uint32_t n = 347'306'346;
uint32_t p = next_prime_number(n);
cout << "\nNEXT PRIME(" << n << ") = " << p << endl;
}- ho definito il bitset come variabile statica perché le variabili automatiche (o dinamiche che dir si voglia) sono molto più restrittive in termini di memoria a disposizione. Ciò però implica anche che il bitset vada resettato ogni volta nelle chiamate a
espandi_e_applica_crivello() successive alla prima;- nella prima chiamata a
espandi_e_applica_crivello() i valori di cui eliminare i multipli sono prelevati dal bitset stesso, mentre nelle chiamate successive vengono prelevati dal vettore di interi contenente i numeri primi, visto che sicuramente il massimo valore di cui eliminare i multipli si troverà in esso e non nel bitset;- tenendo conto del teorema dei numeri primi ho aggiunto la seguente riga di codice
v.reserve(n/log(n)) per evitare troppe reallocazioni del vettore v , ma alla fine l'ho disattivata in quanto non ho riscontrato miglioramenti rilevanti;- a parità di memoria utilizzata per immagazzinare i nuovi valori da sottoporre a crivellatura nella funzione
espandi_e_applica_crivello() , i bitset consentono di testare ad ogni chiamata 32 e 64 volte i valori immagazzinabili in variabili intere rispettivamente a 32 e 64 bit. In queste condizioni ciò ovviamente significa un minor numero di chiamate a espandi_e_applica_crivello() e quindi un bel passo in avanti in termini di velocità.L'ho scritto abbastanza velocemente, se riscontrate qualche errore o vi viene in mente qualche miglioria, fatemi sapere.
- Messaggi
- 25,070
- Reazioni
- 12,679
- Punteggio
- 254
tuttavia per lo stesso intervallo con un bitset bastano circa 122,5 KiB e tenendo conto che usi solo i dispari e che quindi l'intervallo è di 2 milioni di numeri, con il bitset bastano 250 KiB.
La cosa seccante del bitset è che la sua dimensione deve essere nota a tempo di compilazione, quindi serve definire una costante per la sua dimensione, purtroppo l'ho letto 5 minuti fa
quindi per l'espansione bisogna usare un vettore di booleani, sto vedendo che c'è una classe specializzata vector<bool> che si comporta in modo simile ad un bitset. Altra limitazione è che la dimensione non può essere eccessiva, quindi va scelta con criterio (almeno per la primissima inizializzazione).Qualcosa l'ho fatto
Includo un codice di "assaggio" del bitset, con miglioramenti rispetto a quanto dicevo nel mio post precedente:
- rimane la corrispondenza biunivoca che p(i)=1 se i è primo (utile perché evita di fare calcoli sugli indici), tuttavia gli indici pari non vengono analizzati (il numero 2, unico primo pari, è inserito direttamente nel risultato finale)
- quando i è dispari e p(i)=1 -> i è primo, viene inserito nel risultato e si crivellano gli indici dispari multipli di i (per es, per il 3 si ha p[3]=1, 3 si inserisce nel risultato e si marcano con 0 le componenti p[9], p[15], p[21]); quindi in generale se p=1 si marcano p[3*i],=0, p[5*i]=0, p[7*i]=0 ecc.
- Ci saranno inevitabilmente indici dispari marcati più volte con 0 (significa che il numero dispari corrispondente è divisibile per almeno 2 numeri primi).
- Ogni intero non primo ha almeno un divisore primo <= alla sua radice quadrata -> si crivella solo per i primi <= della radice quadrata del massimo intero dell'intervallo considerato.
- Il risultato finale è un vector con tutti i numeri primi nell'intervallo [0; n-1] (dove n è la dimensione del bitset)
Nel codice esempio ho deifinito una dimensione costante grande (2 097 152 = 2^21) ma l'idea è quella di usare un bitset molto più piccolo e fare le espansioni a partire da esso. Per la cronaca sul mio PC-catorcio viene crivellato l'intervallo [0 - 2 097 151], trovati i suoi 155611 primi in un tempo variabile tra i 47 ed i 70 ms (dipende da quanto rimane in cache).
C++:
#include <iostream>
#include <iomanip> // serve per setw
#include <bitset>
#include <cmath> // per radice quadrata
#include <vector>
using namespace std;
// dimensione intervallo e bitset
const unsigned dim = 2097152; // interi <= 2097151
vector<int> inizializza(){
vector<int> primi; // il vettore da ritornare
primi.push_back(2); // il 2 è l'unico primo pari
bitset<dim> p = p.set(); // imposta a 1 tutti i bit
int limiteCrivello = (int)floorl(sqrtl(dim-1)); // max indice per crivello
for(int i = 3; i<dim; i += 2) // solo indici dispari
if(p[i]){ // se i è primo
primi.push_back(i); // aggiunge il primo al vettore risultato
if(i<=limiteCrivello)
for(int j = 3*i; j<dim; j += 2*i) // crivello
p[j] = 0;
}
return primi;
}
// riga orizzontale da 80 caratteri
const string strHR =
"--------------------------------------------------------------------------------\n";
int main(void) {
cout << strHR;
cout << "Crivello di Eratostene per l'intervallo [0; " << dim-1 << "]" << endl;
cout << strHR;
vector<int> primi = inizializza();
cout <<"Primi trovati: " << primi.size() << endl;
}--------------------------------------------------------------------------------
Crivello di Eratostene per l'intervallo [0; 2097151]
--------------------------------------------------------------------------------
Primi trovati: 155611
Process returned 0 (0x0) execution time : 0.047 s
Press any key to continue.
Crivello di Eratostene per l'intervallo [0; 2097151]
--------------------------------------------------------------------------------
Primi trovati: 155611
Process returned 0 (0x0) execution time : 0.047 s
Press any key to continue.
concordo, è così che va fatto, infatti l'intenzione è proprio di usare il bitset come prima inizializzazione, poi il vector risultato, combinato con un'espansione fatta con vector<bool> o con bitset (se solo ne capissi di C++)
non lo avrei fatto per un altro motivo: il teorema dei numeri primi dà una stima del numero di primi <=N ma non il numero esatto, quindi qualche riallocazione ulteriore in genere sarebbe necessaria, a questo punto tanto vale non farla affatto.
P.S.
casomai non si fosse capito, il C++ non lo conosco proprio, mi sono dovuto arrangiare
M1n021
Utente Attivo
- Messaggi
- 279
- Reazioni
- 97
- Punteggio
- 43
In realtà il bitset va bene. L'espansione riguarda i nuoviLa cosa seccante del bitset è che la sua dimensione deve essere nota a tempo di compilazione, quindi serve definire una costante per la sua dimensione, purtroppo l'ho letto 5 minuti fa
N valori da sottoporre al crivello, dove N è una costante.Come detto nel precedente post, se si definisce la variabile come statica si ha a disposizione molta più memoria.
Però non tieni conto del fatto che, considerando per esempio
i=13, è inutile partire 13*3, perché 13*3, 13*5, 13*7, 13*9 e 13*11 sono già stati cancellati in precedenza. Quindi tanto vale partire da i^2 . Inoltre a questo punto quello che tu chiami
limiteCrivello diventa inutile, in quanto la condizione per terminare direttamente il ciclo più esterno diventa i*i>dim-1 (i nomi delle variabili si riferiscono al codice da te postato).Non ho mai approfondito per bene come funziona l'incremento automatico di capacità dei vector, ma ipotizzo che quella riga di codice dovrebbe in teoria scongiurare la maggior parte delle reallocazioni che si verificherebbero in sua assenza, non saprei...
L'importante è l'algoritmo di base, i linguaggi sono un di più!casomai non si fosse capito, il C++ non lo conosco proprio, mi sono dovuto arrangiare
- Messaggi
- 25,070
- Reazioni
- 12,679
- Punteggio
- 254
eh sì, le ottimizzazioni sul linguaggio le lascio a te, il C++ non lo conosco per cui il codice è arrangiato in base alle scarne info che ho, preferisco concentrarmi sull'algoritmo. Abbiamo postato quasi in contemporanea, credo a distanza di un paio di secondi, non avevo letto il tuo ultimo post e mi riferivo al precedente
Vero! e all'inizio è così che avevo fatto ed è andato tutto bene fino ad una certa dimensione, poi perde una manciata di primi
: ho fatto un po' di controlli con Mathematica online (https://www.wolframalpha.com/input/?i=calculator) e quando iniziavo con i=p^2 su intervalli grandi (a memoria con estremo superiore >16 milioni) perde per strada 4 numeri primi, segno che è "scappato" qualcosa, probabilmente ho fatto un errore di codifica che non vedo solo io, magari ci riprovo perché con i=p^2 tagliavo un bel po' di marcature inutili.Per la cronaca, con Mahtematica il controllo sulla quantità di numeri primi la fai mettendo come input PrimePi[n]
deformazione mentale:, quando sono in gioco numeri grandi non lo faccio mai:
in questo caso porta inefficienza perché la moltiplicazione i*i del ciclo potrebbe essere eseguita (decine di) migliaia/milioni di volte --> conviene usare la funzione di libreria radice quadrata che esegue al più qualche centinaio di istruzioni una sola volta (bisogna usare un funzione che abbia una precisione in cifre maggiore degli interi in gioco);
in generale i*i>n può portare errori runtime/overflow quando i si avvicina al massimo intero rappresentabile, per es. se usi un (u)int a 32 bit hai un limite massimo di (4)2 milardi e spiccioli: eleva al quadrato un numero "piccolo" come 100 000 e vai fuori range!
In questo caso particolare non succederebbe per la dimensione relativamente piccola del bitset che volevo usare, però in generale può succedere.
Ho un sacco di problemi di connessione, non so se potrò ricollegarmi più tardi
M1n021
Utente Attivo
- Messaggi
- 279
- Reazioni
- 97
- Punteggio
- 43
Magari se posti il codice in questione possiamo darci un'occhiata.
Considera
i*i>dim-1 come condizione di uscita del ciclo più esterno:- se non si verifica, allora vai ad eliminare i multipli di
i a partire da i^2 , e quindi il calcolo di i*i ti serve comunque;- se invece si verifica allora abbiamo finito con la crivellatura.
Quindi il calcolo di
i^2 risulta "inutile" una sola volta e l'alternativa sarebbe il calcolo di una radice quadrata. Questo dimostra che il metodo che ho descritto (ossia l'utilizzo della condizione i*i>dim-1 ) risulta più efficiente.Riguardo al rischio di "overflow" nel calcolo di
i*i hai ragione; se il crivello lavora su uint32_t basta inserire il suddetto prodotto in un uint64_t.
Ultima modifica:
Pubblicità