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[j] = true;
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.
MAX | Pi(MAX) | Temps (sec) |
1 000 000 | 78 498 | 0.031 |
10 000 000 | 664 579 | 0.657 |
100 000 000 | 5 761 455 | 7.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 = (int) Math.sqrt(MAX) + 1;
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 * i) / 2;
while (j < MEMORY_SIZE) {
array[j] = true;
j += i;
}
}
}
Globalement nous divisons le temps par 2.
MAX | Pi(MAX) | Temps (sec) |
1 000 000 | 78 497 | 0.016 |
10 000 000 | 664 578 | 0.281 |
100 000 000 | 5 761 454 | 3.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 = (long) Math.sqrt(MAX) + 1;
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 += (2 * i);
}
}
}
//--//
public static boolean getBit(long i) {
byte block = array[(int) (i >> 4)];
byte mask = (byte) (1 << ((i >> 1) & 7));
return ((block & mask) != 0);
}
public static void setBit(long i) {
int index = (int) (i >> 4);
byte block = array[index];
byte mask = (byte) (1 << ((i >> 1) & 7));
array[index] = (byte) (block | mask);
}
Et non atteignons le milliard....
MAX | Pi(MAX) | Temps (sec) |
1 000 000 | 78 497 | 0.031 |
10 000 000 | 664 578 | 0.203 |
100 000 000 | 5 761 454 | 3.047 |
1 000 000 000 | 50 847 533 | 35.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, (byte) 0);
long debut = numSegment * MAX_SEGMENT;
long fin = debut + MAX_SEGMENT;
long finSqrt = (long) Math.sqrt(fin) + 1;
int p = -1;
int i = 1;
while (p < finSqrt) {
p = amorce[i++];
long j = p - (debut % p);
j += ((j & 1) == 0 ? p : 0);
long p2 = p * 2;
while (j < MAX_SEGMENT) {
segment[(int) (j >> 4)] |= (1 << ((j >> 1) & 7));
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.
MAX | Pi(MAX) | Temps |
10 000 000 000 | 455 052 511 | 5mn 24s |
100 000 000 000 | 4 118 054 813 | 1h06mn |
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%.
MAX | Pi(MAX) | Temps |
10 000 890 900 | 455 091 136 | 4mn 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.
MAX | Pi(MAX) | Temps |
10 000 890 900 | 455 091 136 | 3mn 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...