Le crible d'Eratosthène classique
Description
La probabilité de découvrir aujourd'hui une amélioration significative de l'algorithme classique du Crible d'Eratosthène est très mince car rares sont les programmeurs ou développeurs qui ne se sont jamais intéressés de près ou de loin à ce procédé.
Malgré tout, je me décide de présenter ici un développement propre (inédit ?) dont les derniers codes divisent les temps de calcul carrément par 3 ou plus !
Par classique, on entend (souvent) un code ressemblant à EratoA ci-dessous, et qui comporte les particularités suivantes:
a) Le crible est complet, sans suppression préliminaire (p.ex. nombres pairs).
b) L'élimination des candidats s'arrête lorsque p*p<=mx.
c) L'élimination des candidats se limite au multiples de nombres premiers if (isPrime[p]) ...
Une amélioration, qui devrait être "classique" aussi, mais souvent oubliée, est de commencer l'élimination avec le carré du nombre premier considéré: for (uint q=p*p; ....
Il est vrai que le gain obtenu n'est que de quelques pourcents:
Ce code très sobre est par exemple proposé dans primesieve pour déterminer les "petits" nombres premiers, qui sont ensuite utilisés pour déterminer le crible segmenté.
Malgré la simplicité ces quelques lignes, on rencontre souvent des codes qui remplissent des pages !
Selon mon "intuition", le prochain code EratoC, qui, sous une forme "équivalente" à EratoB, utilise la multiplication p*q, devrait être moins performant.
Mais les mesures montrent qu'ils sont tout à fait comparables.
En travaillant dans d'autres articles sur les nombres premiers (cycles et multiplications), il m'est venu l'idée de visualiser les multiples de p pour lesquels l'élimination isPrime[p*q]=false est vraiment nécessaire.
Par exemple, pour p=5, j'ai obtenu p*q=25,35,55,65,85,etc. En regardant de plus près, on peut constater que p*q=5*5,5*7,5*11,5*13,5*17,... , donc les produits de p=5 par les nombres premiers >= 5.
Cela ne nous avance pas à grand chose, car ces nombres premiers, on est justement entrain de les chercher !
Mais pour éviter d'éliminer tous les multiples de p à partir de p*p, on peut se limiter aux produits p*q pour lesquels q n'est pas un multiple d'un nombre premier inférieur à p, c'est-à-dire pour ceux qui restent candidats (isPrime[q]=true).
De cette manière, on s'approche de l'idée du Crible d'Euler qui n'élimine un candidat qu'une seule fois.
Malheureusement, ce code donne des résultats faux !!!. En effet, par exemple pour p=2, après avoir éliminé 4, on n'éliminera plus 8=2*4 !
Après avoir laissé "dormir" ce projet pendant quelques temps, je me suis demandé si l'erreur ci-dessus ne provenait pas simplement de l'ordre dans lequel les candidats sont éliminés.
Et, en effet, en inversant la boucle intérieure, la méthode devient correcte ("Ouf !!!"):
Par rapport à la version classique EratoB, je m'attendais pour EratoD à une diminution du temps d'exécution de quelques %, ou peut-être même de quelques dizaines de %.
Je suis donc extrêmement surpris que ce gain est nettement plus élevé (voir Mesures à la fin de ce texte).
En résumé, on obtient un algorithme 3 à 4 fois plus rapide en remplaçant la boucle intérieure:
for (uint q=p*p; q<=mx; q+=p) isPrime[q]=false;
par:
for (uint q=mx/p; q>=p; --q) if (isPrime[q]) isPrime[p*q]=false;
La tentative d'éviter la multiplication dans la boucle intérieure n'améliore pas (ou peu) le temps de calcul.
La soustraction pq-=p doit en effet être faite pour tous les multiples de p, alors que le produit p*q, incontestablement plus lourd, n'est effectué que pour les candidats q (multiples de p) non éliminés.
Pour essayer de rester "simple", je me suis limité jusqu'ici au crible d'Ératosthène classique (complet).
Pour diviser par deux la mémoire utilisée, ne considérons maintenant que les nombres pairs à l'aide de la relation suivante entre l'index i du crible et le nombre correspondant p:
Le temps de calcul par rapport à EratoD est divisé par environ 1.65 !
Maintenant que j'ai le nouveau code sous les yeux (EratoD), l'amélioration apportée me semble tellement "simple" et "évidente", que j'ai de la peine de croire qu'elle n'ait pas été "découverte" et publiée jusqu'à ce jour.
Merci d'avance de me faire parvenir les coordonnées de développements semblables ou complémentaires.
Une prochaine étape sera d'optimiser la mémoire en n'utilisant qu'un bit par élément du crible, ainsi que d'éliminer d'avance non seulement les nombres pairs, mais aussi les multiples de 3, 5, 7, ... .
Il me reste aussi à adapter mes articles actuels et futurs concernant les nombres premiers, dont une liste, complétée de mesures des temps d'exécution, peut être obtenue ici: Nombres premiers.
Les fichiers sieve.cpp et EratoOdd.cpp vous permettent d'essayer EratoA à Erato 4 et EratoOdd,
dont voici ce que me donne mon ordi (i5 CPU 650 3.2 GHz):
Malgré tout, je me décide de présenter ici un développement propre (inédit ?) dont les derniers codes divisent les temps de calcul carrément par 3 ou plus !
Le crible d'Ératosthène classique
Par classique, on entend (souvent) un code ressemblant à EratoA ci-dessous, et qui comporte les particularités suivantes:
a) Le crible est complet, sans suppression préliminaire (p.ex. nombres pairs).
b) L'élimination des candidats s'arrête lorsque p*p<=mx.
c) L'élimination des candidats se limite au multiples de nombres premiers if (isPrime[p]) ...
bool *EratoA(uint mx) {
bool *isPrime=(bool*)malloc(mx+1);
memset(isPrime,true,mx+1); isPrime[0]=false; isPrime[1]=false;
for (uint p=2; p*p<=mx; ++p)
if (isPrime[p])
for (uint q=p+p; q<=mx; q+=p) isPrime[q]=false;
return isPrime;
}
Une amélioration, qui devrait être "classique" aussi, mais souvent oubliée, est de commencer l'élimination avec le carré du nombre premier considéré: for (uint q=p*p; ....
Il est vrai que le gain obtenu n'est que de quelques pourcents:
bool *EratoB(uint mx) {
bool *isPrime=(bool*)malloc(mx+1);
memset(isPrime,true,mx+1); isPrime[0]=false; isPrime[1]=false;
for (uint p=2; p*p<=mx; ++p)
if (isPrime[p])
for (uint q=p*p; q<=mx; q+=p) isPrime[q]=false;
return isPrime;
}
Ce code très sobre est par exemple proposé dans primesieve pour déterminer les "petits" nombres premiers, qui sont ensuite utilisés pour déterminer le crible segmenté.
Malgré la simplicité ces quelques lignes, on rencontre souvent des codes qui remplissent des pages !
Essais personnels
Selon mon "intuition", le prochain code EratoC, qui, sous une forme "équivalente" à EratoB, utilise la multiplication p*q, devrait être moins performant.
Mais les mesures montrent qu'ils sont tout à fait comparables.
bool *EratoC(uint mx) {
bool *isPrime=(bool*)malloc(mx+1);
memset(isPrime,true,mx+1); isPrime[0]=false; isPrime[1]=false;
for (uint p=2; p*p<=mx; ++p)
if (isPrime[p])
for (uint q=p; q<=mx/p; ++q) isPrime[p*q]=false;
return isPrime;
}
En travaillant dans d'autres articles sur les nombres premiers (cycles et multiplications), il m'est venu l'idée de visualiser les multiples de p pour lesquels l'élimination isPrime[p*q]=false est vraiment nécessaire.
Par exemple, pour p=5, j'ai obtenu p*q=25,35,55,65,85,etc. En regardant de plus près, on peut constater que p*q=5*5,5*7,5*11,5*13,5*17,... , donc les produits de p=5 par les nombres premiers >= 5.
Cela ne nous avance pas à grand chose, car ces nombres premiers, on est justement entrain de les chercher !
Mais pour éviter d'éliminer tous les multiples de p à partir de p*p, on peut se limiter aux produits p*q pour lesquels q n'est pas un multiple d'un nombre premier inférieur à p, c'est-à-dire pour ceux qui restent candidats (isPrime[q]=true).
De cette manière, on s'approche de l'idée du Crible d'Euler qui n'élimine un candidat qu'une seule fois.
bool *EratoFaux(uint mx) { // Faux
bool *isPrime=(bool*)malloc(mx+1);
memset(isPrime,true,mx+1); isPrime[0]=false; isPrime[1]=false;
for (uint p=2; p*p<=mx; ++p)
if (isPrime[p])
for (uint q=p; q<=mx/p; ++q) if (isPrime[q]) isPrime[p*q]=false; // Faux
return isPrime; // Faux
}
Malheureusement, ce code donne des résultats faux !!!. En effet, par exemple pour p=2, après avoir éliminé 4, on n'éliminera plus 8=2*4 !
Amélioration personnelle (inédite ?)
Après avoir laissé "dormir" ce projet pendant quelques temps, je me suis demandé si l'erreur ci-dessus ne provenait pas simplement de l'ordre dans lequel les candidats sont éliminés.
Et, en effet, en inversant la boucle intérieure, la méthode devient correcte ("Ouf !!!"):
bool *EratoD(uint mx) {
bool *isPrime=(bool*)malloc(mx+1);
memset(isPrime,true,mx+1); isPrime[0]=false; isPrime[1]=false;
for (uint p=2; p*p<=mx; ++p)
if (isPrime[p])
for (uint q=mx/p; q>=p; --q) if (isPrime[q]) isPrime[p*q]=false;
return isPrime;
}
Par rapport à la version classique EratoB, je m'attendais pour EratoD à une diminution du temps d'exécution de quelques %, ou peut-être même de quelques dizaines de %.
Je suis donc extrêmement surpris que ce gain est nettement plus élevé (voir Mesures à la fin de ce texte).
En résumé, on obtient un algorithme 3 à 4 fois plus rapide en remplaçant la boucle intérieure:
for (uint q=p*p; q<=mx; q+=p) isPrime[q]=false;
par:
for (uint q=mx/p; q>=p; --q) if (isPrime[q]) isPrime[p*q]=false;
La tentative d'éviter la multiplication dans la boucle intérieure n'améliore pas (ou peu) le temps de calcul.
La soustraction pq-=p doit en effet être faite pour tous les multiples de p, alors que le produit p*q, incontestablement plus lourd, n'est effectué que pour les candidats q (multiples de p) non éliminés.
bool *EratoE(uint mx) {
bool *isPrime=(bool*)malloc(mx+1);
memset(isPrime,true,mx+1); isPrime[0]=false; isPrime[1]=false;
for (uint p=2; p*p<=mx; ++p)
if (isPrime[p])
for (uint q=mx/p,pq=p*q; q>=p; --q,pq-=p) if (isPrime[q]) isPrime[pq]=false;
return isPrime;
}
Suppression préliminaire des nombres pairs
Pour essayer de rester "simple", je me suis limité jusqu'ici au crible d'Ératosthène classique (complet).
Pour diviser par deux la mémoire utilisée, ne considérons maintenant que les nombres pairs à l'aide de la relation suivante entre l'index i du crible et le nombre correspondant p:
// correspondance: i=p/2; p=2 si i=0 et p=2*i+1 si i>0;
// i: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 ...
// p: 2 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 ...
bool *EratoOdd(uint mx) {
bool *isPrime=(bool*)malloc((mx+1)/2);
memset(isPrime,true,(mx+1)/2);
for (uint i=1,p=3; p*p<=mx; ++i,p+=2)
if (isPrime[i]) {
uint a=mx/p; // a must be odd !
for (uint j=((a&1)?a:a-1)/2; j>=i; --j)
if (isPrime[j]) isPrime[i+p*j]=false;
}
return isPrime;
}
Le temps de calcul par rapport à EratoD est divisé par environ 1.65 !
Conclusion
Maintenant que j'ai le nouveau code sous les yeux (EratoD), l'amélioration apportée me semble tellement "simple" et "évidente", que j'ai de la peine de croire qu'elle n'ait pas été "découverte" et publiée jusqu'à ce jour.
Merci d'avance de me faire parvenir les coordonnées de développements semblables ou complémentaires.
Une prochaine étape sera d'optimiser la mémoire en n'utilisant qu'un bit par élément du crible, ainsi que d'éliminer d'avance non seulement les nombres pairs, mais aussi les multiples de 3, 5, 7, ... .
Il me reste aussi à adapter mes articles actuels et futurs concernant les nombres premiers, dont une liste, complétée de mesures des temps d'exécution, peut être obtenue ici: Nombres premiers.
Mesures
Les fichiers sieve.cpp et EratoOdd.cpp vous permettent d'essayer EratoA à Erato 4 et EratoOdd,
dont voici ce que me donne mon ordi (i5 CPU 650 3.2 GHz):
Output de Erato.cpp EratoA: (classique) m=10, n=4, t=0 ms m=100, n=25, t=0 ms m=1000, n=168, t=0 ms m=10000, n=1229, t=0 ms m=100000, n=9592, t=0 ms m=1000000, n=78498, t=0 ms m=10000000, n=664579, t=94 ms m=100000000, n=5761455, t=1326 ms m=1000000000, n=50847534, t=16411 ms EratoB: (classique +) m=10, n=4, t=0 ms m=100, n=25, t=0 ms m=1000, n=168, t=0 ms m=10000, n=1229, t=0 ms m=100000, n=9592, t=0 ms m=1000000, n=78498, t=0 ms m=10000000, n=664579, t=94 ms m=100000000, n=5761455, t=1248 ms m=1000000000, n=50847534, t=15085 ms EratoC: (classique avec multiplication) m=10, n=4, t=0 ms m=100, n=25, t=0 ms m=1000, n=168, t=0 ms m=10000, n=1229, t=0 ms m=100000, n=9592, t=0 ms m=1000000, n=78498, t=0 ms m=10000000, n=664579, t=78 ms m=100000000, n=5761455, t=1232 ms m=1000000000, n=50847534, t=15071 ms EratoD: (personnel avec multiplication) m=10, n=4, t=0 ms m=100, n=25, t=0 ms m=1000, n=168, t=0 ms m=10000, n=1229, t=0 ms m=100000, n=9592, t=0 ms m=1000000, n=78498, t=0 ms m=10000000, n=664579, t=31 ms m=100000000, n=5761455, t=421 ms m=1000000000, n=50847534, t=4680 ms EratoE: (personnel sans multiplication) m=10, n=4, t=0 ms m=100, n=25, t=0 ms m=1000, n=168, t=0 ms m=10000, n=1229, t=0 ms m=100000, n=9592, t=0 ms m=1000000, n=78498, t=0 ms m=10000000, n=664579, t=31 ms m=100000000, n=5761455, t=437 ms m=1000000000, n=50847534, t=4627 ms Output de EratoOdd.cpp EratoOdd: (Sans nombres pairs) m=10, n=4, t=0 ms m=100, n=25, t=0 ms m=1000, n=168, t=0 ms m=10000, n=1229, t=0 ms m=100000, n=9592, t=0 ms m=1000000, n=78498, t=0 ms m=10000000, n=664579, t=15 ms m=100000000, n=5761455, t=265 ms m=1000000000, n=50847534, t=2808 ms
Vous n'êtes pas encore membre ?
inscrivez-vous, c'est gratuit et ça prend moins d'une minute !
Les membres obtiennent plus de réponses que les utilisateurs anonymes.
Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.
Le fait d'être membre vous permet d'avoir des options supplémentaires.