logo prime i.t. Prime I.T.
Soutenir Prime I.T. : Bitcoin Litecoin Dodgecoin English version

Le crible d'Eratosthène

Le principe

Le crible d'Eratosthène est une méthode permettant de trouver tous les nombres premiers jusqu'à une borne fixée N.
Pour cela, il faut prendre les nombres de 2 à N et éliminer successivement tous ceux qui sont multiples d'un entier et donc pas premier (un nombre non premier est dit composite).

L'élimination fonctionne de la façon suivante :

  • Prendre le premier élément du tableau, c'est à dire 2 et éliminer tous les multiples de 2 du tableau.
  • Puis prendre le nombre suivant non éliminé, ici 3 et éliminer les multiples de 3.
  • Puis prendre le nombre suivant non éliminé, ici 5 (car 4 a été éliminé à la première étape) et éliminer les multiples de 5.
  • etc...

A la fin les non-éliminés sont les nombres premiers

Première implémentation

D'un point de vue algorithmique, il nous faut donc plusieurs choses :

  • Un tableau de booléen de 2 à N permettant de savoir si un nombre est éliminé ou non.
  • Une boucle parcourant ce tableau.
  • Une deuxième boucle imbriquée éliminant les multiples pour un nombre sélectionné.


Ce qui donne cette première version relativement naïve (écrite en java) :

private static final int MAX = 100000000;
private static boolean[] array = new boolean[MAX];
//--//
for (int i = 2; i < MAX; i++) {
  if (!array[i]) {
    int j = i + i;
    while (j < MAX) {
      array[jtrue;
      j += i;
    }
  }
}

Note : Tout au long de la présentation des algorithmes, nous garderons le même formalisme pour les variables "clés", à savoir :

  • array représente le tableau des booléens avec la valeur true signifiant éliminé (c'est un nombre composite) et false non éliminé (c'est potentiellement un premier).
  • MAX représente la taille maximum de la liste des nombres premiers que l'on souhaite trouver (c'est la borne fixée N).


Sur une vielle machine AMD 3200+ 2Ghz et en java (pas forcément réputé pour sa vélocité), nous arrivons déjà à une vitesse intéressante. Pi(x) représente le nombre de nombres premiers jusqu'à x.

MAXPi(MAX)Temps (sec)
1 000 00078 4980.031
10 000 000664 5790.657
100 000 0005 761 4557.578

C'est a dire, nous avons trouvé les 5.7 millions nombres premiers compris entre 2 et 100 millions en 7 secondes. source
C'est bon mais pas suffisant...

Deuxième implémentation

Cette première version de l'algorithme est largement optimisable. Effectivement il est possible de travailler sur deux axes d'améliorations à savoir trouver plus vite les nombres en accélérant la vitesse du traitement (la double boucle doit être plus rapide) et/ou agrandir la plage de nombre premier à trouver (pouvoir agrandir la borne N)

Pour la vitesse de traitement, il "suffit simplement" que les 2 boucles ne parcours pas les nombres qui sont déjà trivialement éliminés...Ok mais comment ?
Si nous déroulons l'algorithme à la main, nous avons :

  • Prendre le nombre 2 et éliminer les multiples : 2*2, 2*3, 2*4, 2*5, 2*6...
  • Prendre le nombre 3 et éliminer les multiples : 3*2, 3*3, 3*4, 3*5, 3*6...
  • Prendre le nombre 5 et éliminer les multiples : 5*2, 5*3, 5*4, 5*5, 5*6...
  • Prendre le nombre n et éliminer les multiples : n*2, n*3, n*4, n*5, n*6...


Nous voyons que l'algorithme tente d'éliminer trop de multiples, par exemple 3*2 de l'étape 2 a déjà été supprimé dans l'étape 1 (2*3), de même avec 5*2, 5*3 de l'étape 3 avec 2*5 de l'étape 1 et 3*5 de l'étape 2. En conclusion, pour un nombre n donnée, il est inutile de chercher à éliminer les multiples en dessous de son carré (n*n) car ceux-ci ont déjà été éliminés dans les tours précédents. Cela se traduit donc par une augmentation de la borne inférieure de la deuxième boucle qui passe de "i+i" à "i*i" et à un abaissement de la borne supérieure de la première boucle qui passe à racine carré de N.

Pour le deuxième point d'optimisation sur l'agrandissement de la borne N, nous n'allons pas sauvegarder les nombres pairs. Effectivement, hormis le chiffre 2 tous les autres nombres pairs sont composites. Avec cette astuce évidente nous divisons le stockage par deux.

Ce qui donne cette deuxième version :

private static final int MAX = 100000000;
private static final int SQRT_MAX = (intMath.sqrt(MAX1;
private static final int MEMORY_SIZE = (int) (MAX / 2);
private static boolean[] array = new boolean[MEMORY_SIZE];
//--//
for (int i = 3; i < SQRT_MAX; i += 2) {
  if (!array[(i / 2)]) {
    int j = (i * i2;
    while (j < MEMORY_SIZE) {
      array[jtrue;
      j += i;
    }
  }
}

Globalement nous divisons le temps par 2.

MAXPi(MAX)Temps (sec)
1 000 00078 4970.016
10 000 000664 5780.281
100 000 0005 761 4543.063

Les plus assidus remarqueront que nous perdons 1 nombre premier par rapport à la première implémentation (colonne Pi). Tout simplement parce que cet algorithme ne stocke plus les nombres pairs et donc 2 est omis. source

Troisième implémentation

Avec cette deuxième implémentation, nous voyons que nous arrivons à une limite fondamentale du crible d'Eratosthène. Sans toucher aux paramètres de la jvm java, si nous essayons de passer à 1 milliard, nous obtenons :

java.lang.OutOfMemoryError: Java heap space
  at org.cpas.prime.Eratos2.<clinit>(Eratos2.java:12)
Exception in thread "main" 

L'algorithme a besoin de beaucoup de mémoire pour s'exécuter car nous devons stocker les N booléens (ou les N/2 pour la version 2). La complexité en mémoire du crible d'Eratosthène est en O(n).
Essayons donc de stocker plus intelligemment. Au lieu d'avoir un tableau de booléens nous pouvons passer en binaire où chaque bit représentera un nombre éliminé ou non éliminé. Le support binaire sera le byte java égal à 8 bits.

private static final long MAX = 1000000000L;
private static final long SQRT_MAX = (longMath.sqrt(MAX1;
private static final int MEMORY_SIZE = (int) (MAX >> 4);
private static byte[] array = new byte[MEMORY_SIZE];
//--//
for (long i = 3; i < SQRT_MAX; i += 2) {
  if (!getBit(i)) {
    long j = (i * i);
    while (j < MAX) {
      setBit(j);
      j += (* i);
    }
  }
}
//--//
public static boolean getBit(long i) {
  byte block = array[(int) (i >> 4)];
  byte mask = (byte) (<< ((i >> 17));
  return ((block & mask!= 0);
}

public static void setBit(long i) {
  int index = (int) (i >> 4);
  byte block = array[index];
  byte mask = (byte) (<< ((i >> 17));
  array[index(byte) (block | mask);
}

Et non atteignons le milliard....

MAXPi(MAX)Temps (sec)
1 000 00078 4970.031
10 000 000664 5780.203
100 000 0005 761 4543.047
1 000 000 00050 847 53335.625

source

Quatrième implémentation

C'est bien beau tout ça mais cette limite d'espace est toujours présente. Aussi astucieux que nous pouvons l'être, nous aurons toujours à allouer un tableau de 1 à N éléments qui pour un N suffisamment grand dépassera la taille mémoire de l'ordinateur.
Ne partez pas chez Surcouf acheter de la ram, la solution est le crible d'Eratosthène segmenté.

L'idée est de découper le tableau en X morceaux (segments) et de lancer l'algorithme du crible sur chacun de ces morceaux. Non seulement c'est possible mais en plus ça ne "plombe" pas trop les performances de l'application (le temps de traitement de chaque segment étant à peu près équivalent au temps de calcul de la troisième implémentation).
Suite à l'optimisation de la deuxième implémentation, nous savons que pour trouver tous les nombres premiers de 2 à N, nous pouvons nous contenter d'appliquer le crible de 2 à racine carrée de N. Ce qui n'est pas anodin puisque pour trouver tous les nombres premiers entre 2 et 10 milliards, nous avons besoin de parcourir uniquement jusqu'à 100 mille. Le but d'Eratosthène segmenté est donc de calculer cette liste de nombre premier jusqu'à racine carrée de N (l'amorce) et de "cribler les différents segments" avec cette amorce.
Par contre l'optimisation sur le parcours des nombres dans la deuxième boucle change. Dans le cas d'un Eratosthène segmenté c'est plus délicat à mettre en place car le début du "criblage" commence sur le début du segment qui est un nombre choisi arbitrairement. Nous repassons donc d'un "i*i" à un "i+i" mais il doit être possible de faire mieux. Trop de neurones étant en jeu pour optimiser cette partie, je laisse le lecteur trouver la solution...

Remarque intéressante les segments sont indépendants, il n'est donc pas obligatoire de calculer d'abord le premier segment pour en déduire le second. Cela signifie donc qu'il est possible de calculer les nombres premiers entre 2 bornes quelconques (du moment que l'amorce est suffisamment grande pour atteindre la racine carrée de la borne sup du segment), cela signifie aussi qu'il est tout à fait envisageable de créer un Eratosthène distribué entre plusieurs ordinateurs (chacun cherchant son segment).

Arrays.fill(segment, (byte0);

long debut = numSegment * MAX_SEGMENT;
long fin = debut + MAX_SEGMENT;
long finSqrt = (longMath.sqrt(fin1;

int p = -1;
int i = 1;
while (p < finSqrt) {
  p = amorce[i++];
  long j = p - (debut % p);
  j += ((j & 1== ? p : 0);

  long p2 = p * 2;
  while (j < MAX_SEGMENT) {
    segment[(int) (j >> 4)] |= (<< ((j >> 17));
    j += p2;
  }
}

numSegment représente le numéro du segment et MAX_SEGMENT la taille du segment, avec nous en déduisons le début et la fin du segment.

MAXPi(MAX)Temps
10 000 000 000455 052 5115mn 24s
100 000 000 0004 118 054 8131h06mn

Enfin il est à noter que le calcul des nombres premiers de l'amorce dans Eratos4 est fait sous la forme de trial division. Il aurait été plus rapide de le faire avec un Eratosthène (une poigné de millisecondes contre 500ms) mais il faut bien faire varier les plaisirs... source

Cette implémentation fait un dump des segments et les écrit dans des fichiers sous la forme prime_{N° segment}_{Debut segment}.bin
Pour lire ces fichiers, voici un exemple de comment faire : source

Quatrième implémentation (bis)

En analysant le code du crible d'Eratosthène segmenté nous voyons que celui-ci est couteux en temps de calcul au début du criblage d'un segment, lorsque nous appliquons l'élimination avec les petits chiffres de l'amorce. Ceci est dû au fait que... justement les chiffres sont trop petits et qu'il faut donc beaucoup de "i+i" pour atteindre la fin du segment.
Voyez :

  • Sur un segment d'1 milliard, il faut 333 millions de tour de boucle pour éliminer tous les multiples de 3.
  • Pour 5, il faut 200 millions de tour de boucle.
  • Pour 17, nous passons à 58 millions.


Nous allons donc faire une petite optimisation en calculant une fois pour toute ces "i+i" pour les petits chiffres de l'amorce (nommé le mask) puis nous l'appliquerons tel quel lors du criblage d'un segment. A noter que pour "l'appliquer tel quel" il faudra que la taille du segment soit divisible par tous ces petits chiffres. Pour Eratos4bis nous allons pré-calculer un mask jusqu'à 17 et nous allons positionner la taille du segment à 1000089090 = 3 * 5 * 7 * 11 * 13 * 17 * 3918.
Avec cette optimisation l'algorithme perd 4 secondes au démarrage (le temps de calculer le mask) puis regagne ce temps à chaque calcul de segment. Au final le gain est de performance est de l'ordre de 15%.

MAXPi(MAX)Temps
10 000 890 900455 091 1364mn 35s
100 008 909 000-A vérifier-54mn

source

Quatrième implémentation (ter)

Puisque java est un peu lent, essayons de retranscrire Eratos4bis en C.

MAXPi(MAX)Temps
10 000 890 900455 091 1363mn 35s
100 008 909 000-A vérifier-35mn

Effectivement c'est nettement plus rapide !
Testé sous environnement Windows avec une machine 32 bits simple coeur, apparemment ce gain n'est pas si important sur les machines récentes. Peut-être que java arrive à exploiter pleinement les ressources machines contrairement au C et donc à améliorer les performances ? source

Et maintenant ?

C'est bien gentil d'avoir calculé des milliards de nombre premier mais maintenant qu’en faisons-nous ? Un segment d' 1 milliard représente à peu près un fichier de 60Mo (26Mo en zippé). Donc en moins d'une heure de temps de calcul, nous trouvons tous les nombres premiers de 1 à 100 milliards (il y en a 4 milliards) mais surtout avec 6Go de données (ou 2Go compressé).
Cet algorithme nous montre bien qu'il est relativement rapide de générer linéairement tous les nombres premiers entre 2 bornes, par contre les stocker et/ou les utiliser est un autre défi...


 
Un site réalisé par Les composants associés (cpas) logo cpas