Voilà plusieurs fois que je tombe sur des programmes C cherchant à savoir si un nombre est premier ou pas. La plupart du temps, on peut aller beaucoup plus loin que ce qui est proposé.
Petit tour d'horizon des optimisations à appliquer ; pas toujours évidentes, mais qui ont l'avantage d'être bien plus rapides que l'implémentation naïve.
Définition du cahier des charges
Que doit faire le programme ?
Dans cet article, j'examinerai différents algorithmes qui auront pour but de trouver tous les nombres premiers inférieurs à une constante, MAX_NUM et de les renvoyer à l'intérieur d'un tableau d'entiers.
Les contraintes suivantes devront être respectées :
- Code le plus performant possible.
- Pas d'algorithmes compliqués : les codes ne dépasseront pas quarante lignes.
- MAX_NUM restera inférieur à la taille maximale codable sur une variable de type unsigned long.
- On ne se soucie pas de l'affichage du tableau à la fin
- Pour rester général, les codes seront en C.
Télécharger les codes de tests.
Précisons les choses
Nombre premier
Avant d'aller plus loin, définissons le terme de nombre
premier. Il s'agit d'un nombre qui n'admet aucun autre
diviseur que 1 et lui-même ; par convention on ne
considére pas 1 comme premier.
Par exemple, 8 n'est pas premier car il se décompose en
4 × 2, tandis que 2, 3,
5, 7 le sont.
Protocole de test
Chacun des codes ne sera constitué que d'un seul fichier.
Les codes seront testés en C compilé avec gcc sans
aucune option de compilation. La durée d'exécution sera
calculée à l'aide de la fonction UNIX
time.
Ligne de commande de compilation :
gcc fichier.c -lm && time ./a.out
Détails pratiques
Pour les calculs, on règlera MAX_NUM à
1 000 000 ; cela suffit pour faire travailler
le processeur un minimum sans avoir à attendre la fin de la
journée (au total, on devra trouver
78 498 nombres premiers
dans cet intervalle).
Les temps que j'indique proviennent d'un PC de test qui ne
fait tourner que l'application.
Dédicace
Je dédicace cet article à un tyran fort pédagogue qui se reconnaîtra sûrement en lisant cet article dont il est l'inspirateur… hém hém.
Approches intuitives
Je suis jeune et naïf, je code vite et pas optimisé
Première approche : on veut récupérer tous les nombres
premiers inférieurs à MAX_NUM, je vais donc faire une
boucle sur ces nombres. Pour chacun des nombres, j'appellerais
une fonction qui me renverra true si le nombre est premier,
false sinon.
Cette fonction testera tous les nombres
i inférieurs à son paramètre
p et regardera si
i divise
p ; si c'est le cas on
renvoie directement true ; sinon on sort de la boucle et
on renvoie false.
Implémentation : code no1
#include <stdio.h> #define MAX_NUM 1000000 /** * Cette fonction prend en paramètre un nombre et essaie de calculer * s'il est premier ou non. * @param p Le nombre dont on veut tester la primalité * @return 0 si p n'est pas premier, 1 sinon. */ int premier(unsigned long p) { unsigned long i; //Pour chacun des nombres inférieurs à p for(i=2;<p;i++) { //i est-il un diviseur de p ? //NOTE : % dénote l'opération modulo, //qui renvoie le reste de la division euclidienne de p par i. //Exemple : 5 % 2 = 1 (car 5 = 2*2 + 1) if(p % i == 0) return 0; } //Aucun des nombres ne le divise, p est premier return 1; } int main() { unsigned long Premiers[MAX_NUM]; unsigned long premiers_trouves = 0; unsigned long nombre; //On commence la boucle à deux car 1 n'est pas premier. for(nombre=2; nombre<MAX_NUM;nombre++) { if(premier(nombre) == 1) { //Le nombre i est premier, l'ajouter. Premiers[premiers_trouves++] = nombre; } } printf("Total : %lu\n",premiers_trouves); return 0; }
Effectivement, ça marche. Mais c'est long, très long…
Temps du code : 8m54,778s
Je suis jeune et un peu moins naïf, j'ai déjà fait des maths
L'idée : s'arrêter à la racine carrée pendant les tests.
Le code précédent effectue un nombre incroyable de tests
inutiles.
Pour comprendre pourquoi, il va falloir faire un minimum de
maths…
Admettons que notre nombre x s'écrive comme le produit de y et
de z (x = y × z). Que
peut-on dire de y et de
z par rapport à
x ? L'un des deux est
forcément inférieur ou égal à la racine carrée de
x. Il suffit donc de modifier
notre fonction premier pour arrêter la boucle for à
la racine carrée : après, on a forcément déjà testé.
Comme on travaille avec des entiers et que la racine est
rarement entière, on rajoute 1 pour être sûr de ne pas
avoir un arrondi malheureux.
Accessoirement, j'ai aussi modifié la condition :
spécifier == 1 est inutile (mais cela ne change
normalement rien, le compilateur faisant efficacement son
boulot).
Implémentation : code no2
#include <stdio.h> #include <math.h> #define MAX_NUM 1000000 int premier(unsigned long p) { unsigned long i; //Pour chacun des nombres inférieurs à racine(p) unsigned long sqrt_p = sqrt(p) + 1; for(i=2;i<sqrt_p;i++) { if(p % i == 0) return 0; } //Aucun des nombres ne le divise, p est premier return 1; } int main() { unsigned long Premiers[MAX_NUM]; unsigned long premiers_trouves = 0; unsigned long nombre; for(nombre=2; nombre<MAX_NUM;nombre++) { if(premier(nombre)) { Premiers[premiers_trouves++] = nombre; } } printf("Total : %lu\n",premiers_trouves); return 0; }
Mieux, beaucoup mieux. Cette fois, il ne faut qu'un peu plus d'une seconde pour faire les mêmes calculs ; c'est une énorme amélioration ! Mais c'est loin d'être fini.
Temps du code : 0m1,013 s
Je suis jeune et encore moins naïf, j'ai déjà fait des maths et de l'algorithmie
L'idée : les nombres pairs ne sont jamais premiers.
Avez-vous remarqué que l'on continue de faire des calculs
inutiles ? Si si, regardez : on teste pour savoir si
les nombres 4, 6, 8, etc. sont premiers. Mais c'est
stupide ! Le seul nombre pair premier, c'est deux ;
calculer les autres est un gâchis de précieux temps.
De même à l'intérieur de la fonction premier :
pourquoi vérifier si le nombre est divisible par 4 alors
qu'il n'est pas divisible par deux ? Autant sauter le
problème ! On va donc modifier l'itérateur de la boucle
en remplaçant le i++ par i+=2. Et plutôt que
de commencer à 2, on va directement à 3 : on examinera
donc 3, 5, 7…
Cette nouvelle approche fait cependant disparaître le
2 de la liste des nombres ; il va donc falloir
l'ajouter manuellement au tout début.
Implémentation : code no3
#include <stdio.h> #include <math.h> #define MAX_NUM 1000000 int premier(unsigned long p) { unsigned long i; //Pour chacun des nombres impairs inférieurs à racine(p) unsigned long sqrt_p = sqrt(p) + 1; for(i=3;i<sqrt_p;i += 2) { if(p % i == 0) return 0; } return 1; } int main() { unsigned long Premiers[MAX_NUM/2]; unsigned long premiers_trouves = 0; unsigned long nombre; //Placer le 2 dans la liste car la boucle le saute directement. Premiers[premiers_trouves++]=2; //On commence la boucle à trois et on va de nombre impair en nombre impair for(nombre=3; nombre<MAX_NUM;nombre+=2) { if(premier(nombre)) { Premiers[premiers_trouves++] = nombre; } } printf("Total : %lu\n",premiers_trouves); return 0; }
Notez aussi que la taille du tableau contenant les nombres
premiers a été réduite par deux elle-aussi. En théorie on
aurait pu faire une meilleure approximation (en prenant une
borne supérieure de x/ln(x) par exemple), en pratique on se
contentera de cela, n'allons pas tout compliquer.
De façon très logique, on trouve un temps divisé par
deux… qui en est surpris, considérant les changements
de ce nouveau code ?
Temps du code : 0m0,504 s
Je suis jeune, pas trop naïf et j'ai de la suite dans les idées
La dernière fois, nous avons modifié notre code pour sauter
automatiquement les multiples de deux. Mais pourquoi s'arrêter
en si bon chemin ? Sautons aussi les multiples de
3 !
Comme je suppose que la formule ne vous vient pas
intuitivement, examinons la séquence formées des nombres qui
ne sont multiples ni de deux ni de trois : 5, 7, 11, 13,
17, 19, 23, 25…
Attention ! Notez bien qu'il ne s'agit pas de nombres
premiers (regardez 25 par exemple) ; mais de nombres
qui ne sont pas divisibles par deux ou trois.
Vous voyez une relation ? Non ? Mais si, regardez la
différence entre deux nombres consécutifs ! Une fois
deux, une fois quatre, une fois deux, et ainsi de suite. Notre
i+=2 se transforme donc en i+=pas. Mais
comment faire alterner ce pas ? Avec un test
if ? On pourrait, mais on va faire beaucoup plus
rapide – ce qui reste le but principal.
Je vous laisse chercher un peu la formule, c'est toujours la
même quand on veut faire alterner un nombre entre deux
valeurs : (somme des choix) – valeur actuelle.
Ainsi, on initalisera au début notre pas à 21, puis on lui appliquera la formule
(4+2) - pas pour obtenir le
nouveau pas, et ainsi de suite. Essayez, ça marche !
Comme la dernière fois, cette stratégie nous fait gagner du
temps… mais nous empêche de mémoriser 2 et 3,
qu'il faut donc ajouter manuellement dans le tableau.
Implémentation : code no4
#include <stdio.h> #include <math.h> #define MAX_NUM 1000000 int premier(unsigned long p) { unsigned long i; //Pour chacun des nombres impairs inférieurs à racine(p) unsigned long sqrt_p = sqrt(p) + 1; int pas = 4; for(i=5;i<sqrt_p;i += pas) { if(p % i == 0) return 0; //Mettre à jour le pas pas = 6 - pas; } return 1; } int main() { unsigned long Premiers[MAX_NUM/2]; unsigned long premiers_trouves = 0; unsigned long nombre; //Placer les 2 et 3 dans la liste car la boucle les saute directement. Premiers[premiers_trouves++]=2; Premiers[premiers_trouves++]=3; //On commence la boucle à sept et on va de nombre impair en nombre impair unsigned int pas = 4; for(nombre=5; nombre<MAX_NUM; nombre+=pas) { if(premier(nombre)) { Premiers[premiers_trouves++] = nombre; } //Mettre à jour le pas pas = 6 - pas; } printf("Total : %lu\n",premiers_trouves); return 0; }
Vous devriez vous dire qu'on peut continuer ainsi longtemps, en enlevant les diviseurs. Après 2 et 3, on vire le 5, puis le 7… sauf que le calcul du pas augmente alors en complexité, et on y perd tout avantage. En fait, l'approche naïve est maintenant bien optimisée : on peut difficilement faire mieux2.
Temps du code : 0m0,361 s
Passons aux choses sérieuses
La première approche est maintenant rapide, on ne l'améliorera
pas significativement. Pour progresser encore, il faut revoir
notre stratégie et se poser la question :
« qu'est-ce qui fait ramer ? ».
Dans un projet normal, la réponse est rarement aisée à
trouver. Ici, le peu de lignes de codes nous aide : on va
analyser les opérations potentiellement longues.
- Le calcul de la racine carrée : ce n'est pas un calcul primordial (il n'est pas dans la boucle principale), mais on pourrait l'optimiser – par exemple en calculant le développement limité de la fonction racine autour du point considéré, ou plus simplement autour du centre (500 000) et en corrigeant l'erreur pour les points lointains. Mais nous ne sommes pas ici pour faire des maths ! Conclusion : on n'y touche pas.
- L'appel de fonction : c'est sûrement l'un des problèmes principaux : chaque appel de fonction oblige le processeur à empiler et dépiler des valeurs inutilement. Le surcoût est ici difficilement tolérable, car tout le reste du code est constitué d'instructions simples (alors que dans la normale, il est courant d'avoir des instructions dont la durée est significativement plus grande qu'un appel de fonction).
- Les tests en folie : finalement, la seule chose qui reste en dehors des boucles for, c'est ces tests if. Hé bien, nous allons justement les condenser pour en extraire la substantifique moelle plus rapidement…
Retour à la base avec la nouvelle approche
Pour ne pas aller trop vite, je vais reprendre le tout premier code en lui appliquant la nouvelle architecture. Il sera donc plus lent que le code no4, mais plus rapide que le test no1 avec lequel il doit être comparé.
Implémentation : code no2-1
#include <stdio.h> #define MAX_NUM 1000000 int main() { unsigned long Premiers[MAX_NUM/2]; unsigned long premiers_trouves = 0; unsigned long nombre; unsigned long i; for(nombre=2; nombre<MAX_NUM; nombre++) { //Test de primalité direct : for(i=2; nombre % i;i++); //----Notez l'absence de corps : //tout s'effectue dans les conditions du for. if(i == nombre) { Premiers[premiers_trouves++]=i; } } printf("Total : %lu\n",premiers_trouves); return 0; }
Que dire ? D'abord, le code est beaucoup plus concis que
son équivalent « naïf » : le nombre
d'instructions est beaucoup plus limité, et la « zone
critique » (la boucle for imbriquée) n'est pas encombrée
de tests : en fait, la boucle for n'a même pas de corps
(il y a un point virgule à la fin, pas un bloc entouré par des
accolades comme on en a l'habitude).
Note : à chaque passage dans la boucle for, le processeur
teste la condition de sortie et sort de la boucle si elle vaut
0. On profite de ce comportement pour ne même pas effectuer la
comparaison implicite i % nombre ! = 0.
Une petite précision s'impose pour expliquer l'absence de
condition de sortie i < nombre : c'est le
modulo qui joue le même rôle (puisque
Nombre mod Nombre = 0), et ce en
plus de tester la divisibilité. En sortie de boucle, il nous
reste donc à déterminer à qui est dû la fin de la
boucle : la présence d'un diviseur, ou son absence (c'est
le test if).
Temps du code : 8m6,025s. Ça n'en a pas l'air, mais c'est une amélioration de presque 10% !
Toujours plus loin !
L'idée : replacer toutes les améliorations précédentes
dans la nouvelle approche.
Maintenant que le concept de base est maîtrisé (se servir des
conditions des boucles for), on peut ré-implémenter nos idées
précédentes : s'arrêter à la racine carrée, sauter les
multiples de deux et de trois…
Mais pas question de rajouter trop de code ! On peut
faire ça rapidement :
Implémentation : code no2-2
#include <stdio.h> #include <math.h> #define MAX_NUM 1000000 int main() { unsigned long Premiers[MAX_NUM/2]; unsigned long premiers_trouves = 0; unsigned long nombre; unsigned long i; Premiers[premiers_trouves++]=2; Premiers[premiers_trouves++]=3; int pasNombre = 4; for(nombre=5; nombre<MAX_NUM; nombre+=(pasNombre=6-pasNombre)) { int pasI = 4; unsigned long sqrt_nombre=sqrt(nombre) + 1; //Test de primalité direct : for(i=5; (i<sqrt_nombre) && (nombre % i);i+=(pasI = 6 - pasI)); if(i >= sqrt_nombre) { Premiers[premiers_trouves++]=nombre; } } printf("Total : %lu\n",premiers_trouves); return 0; }
Notez la concision du code, au détriment de la lisibilité… en particulier, la deuxième boucle devient peu digeste.
Temps du code : 0m0,348s. C'est quelques millisecondes de plus par rapport à notre précédent record !
Tant et plus ?
Vous pensez qu'on ne fera pas mieux ?
Détrompez-vous !
Comme je le disais plus haut, la stratégie « enlever les
nombres premiers déjà trouvés dans le pas » est payante
mais complexe : la formule de base (+2) s'est transformée
en (+2, +4), et devient de pire en pire. Plutôt que de
continuer ainsi, on va plutôt profiter de notre liste de
nombres premiers en partant du principe que si aucun nombre
premier ne divise x, alors x est premier. Pas besoin de tester
les autres nombres, qui sont tous des multiples de nombre
premiers (par définition même des nombres premiers).
Encore une fois, il suffit de tester les nombres premiers
inférieurs à x1/2.
Implémentation : code no2-3
#include <stdio.h> #include <math.h> #define MAX_NUM 1000000 int main() { unsigned long Premiers[MAX_NUM/2]; unsigned long premiers_trouves = 0; unsigned long nombre; unsigned long i; Premiers[premiers_trouves++]=2; Premiers[premiers_trouves++]=3; int pasNombre = 4; for(nombre=5; nombre<MAX_NUM; nombre+=(pasNombre=6-pasNombre)) { unsigned long pointeurCourant=0; unsigned long sqrt_nombre=sqrt(nombre) + 1; for(i=Premiers[pointeurCourant]; (i<sqrt_nombre) && (nombre % i);i=Premiers[pointeurCourant++]); if(i >= sqrt_nombre) { Premiers[premiers_trouves++]=nombre; } } printf("Total : %lu\n",premiers_trouves); return 0; }
Temps du code : 0m0,211s.
Et pour finir en beauté…
Pour finir en beauté, on va encore pousser l'approche
précédente.
En fait, on va devoir modifier des grosses sections pour faire
fonctionner le code, mais l'idée reste une simple amélioration
du concept précédent.
Soyons simple pour l'instant et dé-considérons la suppression
des multiples de 2 et de 3 ; on revient sur une
boucle for « standard ».
À chaque fois qu'on rencontre un nombre premier, on va stocker
dans un tableau son prochain multiple. Ainsi, si on voit deux,
on stockera dans un tableau 4. Si on voit trois, on mettra 6.
Pour tester si un nombre donné est premier, on ne va pas faire
des tests de divisibilité, juste vérifier qu'il n'est pas dans
la liste des multiples.
Si le nombre qu'on teste est supérieur au multiple
actuellement considéré (par exemple, on teste 5 contre le
multiple de 2 qu'est 4), on ajoute au multiple la valeur du
nombre premier associé (4 se transforme alors en 6 et
redevient utile. Le 4 disparaît mais on n'en a plus
besoin) : ainsi, la table des multiples se remplit au fur
et à mesure des nombres qui ne sont pas premiers ; tous
les autres sont ajoutés dans les deux listes (listes de
nombres premiers et liste des multiples).
Vous avez le principe ? Complexifions : on reprend notre pas qui évite directement les multiples de 2 et 3. En conséquence, on ne trouvera jamais 2 × NombrePremier qui est de façon évidente un multiple de 2. De même, 3 × NombrePremier sera sauté (multiple de 3). 5 × NombrePremier sera trouvé quand on le comparera au multiple de 5 actuellement en mémoire, 7 × NombrePremier aussi… finalement, le prochain multiple qu'on peut rencontrer pour un nombre n est n2 (testez pour n=5 : 10 est éliminé par le pas, 15 aussi, 20 de même, seul 25 arrive). À la découverte d'un nombre premier, on ajoute donc dans la table des multiples ledit nombre premier au carré. Et à chaque fois qu'on passe au dessus (25 est éliminé par comparaison à la table des multiples, valeur suivante : 29. 29>25), on met à jour le multiple. Dans la version simple, on ajoutait simplement la valeur de base (5), mais ici c'est inutile (le multiple est impair [les multiples de deux ne sont pas considérés] et la valeur ajoutée est impaire : on obtiendrait un multiple pair, ce qui ne sert à rien) : on ajoutera donc au multiple deux fois la valeur de base (25 + 2 × 5 = 35).
Toujours pas convaincu par cette « simplification » des multiples ? Prenons 7 : 14 est multiple de 2 et n'est pas considéré, 21 multiple de 3, 28 de deux, 35 n'est multiple ni de deux ni de trois, mais est dans la table des multiples [valeur associée à 5 une fois passé le cap des 25], la prochaine valeur intéressante à considérer est donc bien 72=49).
Allez, c'est pas facile, je vais donner un exemple concret.
Pour initaliser, on place dans la liste des nombres premiers
5 et son premier multiple à considérer : 253.
On lance ensuite la boucle for (qui commence à 7, 5 étant
manuellement ajouté).
nombre | Premiers | Multiples | Commentaire |
---|---|---|---|
nombre | Premiers | Multiples | Commentaire |
7 | {5} | {25} | 7 est-il dans la liste des multiples (qui ne contient pour l'instant que 25) ? Non. Le nombre premier associé à 25 (5) est-il supérieur à la racine carrée de 7 ? Oui ; on sort donc et on ajoute 7 et 49 dans les listes. |
11 | {5,7} | {25,49} | Pareil que pour 7 : 11 n'est pas dans la liste des multiples, on l'ajoute. |
13 | {5,7, 11} | {25,49,121} | |
17 | {5,7, 11,13} | {25,49,121,169} | On continue : pour l'instant tout le monde est premier. |
19 | {5,7, 11,13,17} | {25,49,121,169,…} | J'arrête de noter la table des multiples, vous avez compris le principe. |
23 | {5,7, 11,13,17,19} | {25,49,121,169,…} | |
25 | {5,7, 11,13,17,19,…} | {25,49,121,169,…} | Ah ! Cette fois, on trouve bien 25 dans la liste des multiples. On peut donc affirmer qu'il n'est pas premier ! |
29 | {5,7, 11,13,17,19,…} | {35,49,121,169,…} | 29>25 : on met donc à jour la table des multiples. Quant à 29, il n'est pas dans nos multiples, il est donc premier. |
31 | {5,7, 11,13,17,19,…} | {35,49,121,169,…} | |
35 | {5,7, 11,13,17,19,…} | {35,49,121,169,…} | Dans la table des multiples, donc non premier. |
37 | {5,7, 11,13,17,19,…} | {45,49,121,169,…} | On met à jour le prochain multiple, et on continue. |
Implémentation : code no2-4
#include <stdio.h> #define MAX_NUM 1000000 int main() { unsigned long Premiers[MAX_NUM/2]; unsigned long Multiples[MAX_NUM/2]; unsigned long premiers_trouves = 0; unsigned long nombre; unsigned long i; Multiples[premiers_trouves]=5*5; Premiers[premiers_trouves++]=5; int pasNombre = 2; for(nombre=7; nombre<MAX_NUM; nombre+=(pasNombre=6-pasNombre)) { unsigned long pointeurCourant=0; //Test de primalité direct : for(i=Premiers[pointeurCourant]; (nombre!=Multiples[pointeurCourant]) && i*i < nombre;i=Premiers[++pointeurCourant]) { if(nombre>Multiples[pointeurCourant]) Multiples[pointeurCourant] +=2*i; } if(nombre!=Multiples[pointeurCourant]) { Multiples[premiers_trouves] = nombre * nombre; Premiers[premiers_trouves++] = nombre; } } Premiers[premiers_trouves++]=2; Premiers[premiers_trouves++]=3; printf("Total : %lu\n",premiers_trouves); return 0; }
Temps du code : 0m0,120s.
En conclusion
Cet article vous aura, je l'espère interessé. Il montre
l'importance du choix de l'algorithme, mais aussi la nécessité
de comprendre les mécanismes internes du langage (je pense par
exemple aux appels de fonction, ou aux possibilités offertes
par les boucles for). Soyons franc cependant : améliorer
de 445 000 % son code n'est pas un événement
courant. Mais l'important est de comprendre le
principe…
Pour ceux qui souhaiteraient aller plus loin sachez que l'on
utilise aussi des tests dits probabilistes qui sont éminemment
plus rapides. Je n'en traite pas ici, et par méconnaissance du
sujet, et pour ne pas embrouiller le lecteur en apportant des
idées autrement plus complexes que les algorithmes "de base"
que je présente ici.
Je redirigerai donc le lecteur curieux vers les tests de
primalité de
Solovay-Strassen
et de
Miller-Rabin.
Télécharger les codes de tests.
Première approche | Deuxième approche | ||
---|---|---|---|
Algo | Temps | Algo | Temps |
1 | 8m54.778s | 1 | 8m6.025s |
2 | 0m1.013s | ||
3 | 0m0.504s | ||
4 | 0m0.361s | 2 | 0m0.348s |
3 | 0m0.211s | ||
4 | 0m0.120s |
- 1 ↑ En fait, on l'initalise à 4 pour qu'il soit à deux au moment où on veut l'utiliser.
- 2 ↑ À moins de jouer avec les options du compilateur, mais on se l'est interdit dans le cahier des charges ! À titre d'exemple, avec l'option -O3, on gagne une cinquantaine de millisecondes.
- 3 ↑ Notez que dans cette approche, 2 et 3 vont nous gêner si on les met directement dans la liste. On les ajoutera à la fin.