C / C++ / C++.NET : Le crible d'Eratosthène classique (2024)

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 !

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 "petit*" 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}

Malheureusem*nt, 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.cppEratoA: (classique)m=10, n=4, t=0 msm=100, n=25, t=0 msm=1000, n=168, t=0 msm=10000, n=1229, t=0 msm=100000, n=9592, t=0 msm=1000000, n=78498, t=0 msm=10000000, n=664579, t=94 msm=100000000, n=5761455, t=1326 msm=1000000000, n=50847534, t=16411 msEratoB: (classique +)m=10, n=4, t=0 msm=100, n=25, t=0 msm=1000, n=168, t=0 msm=10000, n=1229, t=0 msm=100000, n=9592, t=0 msm=1000000, n=78498, t=0 msm=10000000, n=664579, t=94 msm=100000000, n=5761455, t=1248 msm=1000000000, n=50847534, t=15085 msEratoC: (classique avec multiplication)m=10, n=4, t=0 msm=100, n=25, t=0 msm=1000, n=168, t=0 msm=10000, n=1229, t=0 msm=100000, n=9592, t=0 msm=1000000, n=78498, t=0 msm=10000000, n=664579, t=78 msm=100000000, n=5761455, t=1232 msm=1000000000, n=50847534, t=15071 msEratoD: (personnel avec multiplication)m=10, n=4, t=0 msm=100, n=25, t=0 msm=1000, n=168, t=0 msm=10000, n=1229, t=0 msm=100000, n=9592, t=0 msm=1000000, n=78498, t=0 msm=10000000, n=664579, t=31 msm=100000000, n=5761455, t=421 msm=1000000000, n=50847534, t=4680 msEratoE: (personnel sans multiplication)m=10, n=4, t=0 msm=100, n=25, t=0 msm=1000, n=168, t=0 msm=10000, n=1229, t=0 msm=100000, n=9592, t=0 msm=1000000, n=78498, t=0 msm=10000000, n=664579, t=31 msm=100000000, n=5761455, t=437 msm=1000000000, n=50847534, t=4627 msOutput de EratoOdd.cppEratoOdd: (Sans nombres pairs)m=10, n=4, t=0 msm=100, n=25, t=0 msm=1000, n=168, t=0 msm=10000, n=1229, t=0 msm=100000, n=9592, t=0 msm=1000000, n=78498, t=0 msm=10000000, n=664579, t=15 msm=100000000, n=5761455, t=265 msm=1000000000, n=50847534, t=2808 ms