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 n^o1

 
#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;i<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 n^o2

 
#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 n^o3

 
#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 à 2¹, 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 n^o4

 
#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 mieux².

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 n^o4, mais plus rapide que le test n^o1 avec lequel il doit être comparé.

Implémentation : code n^o2-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 n^o2-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 à x^1/2.

Implémentation : code n^o2-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 n² (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 7²=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 : 25³.
On lance ensuite la boucle for (qui commence à 7, 5 étant manuellement ajouté).

Implémentation : code n^o2-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.

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.