Dernière modification:
- Source XML: 2003-06-18, 23:39:56 (20)
- Feuille de styles XSL: 2003-08-06, 22:11:54 (488)
Note: article envoyé le 1997-01-27, paru dans ARMada News n°0C. Il s'agit de la version 6 de l'algorithme.
Le problème consiste à trouver le plus rapidement possible tous les nombres premiers inférieurs à un entier M donné. L'algorithme le plus rapide est le crible d'Ératosthène (le principe est rappelé ci-dessous). Nous cherchons à l'implémenter sur un Risc PC actuel muni d'un StrongARM à 202 MHz et avec assez de mémoire RAM; nous n'avons pas de cache de niveau 2. Le résultat aura la forme d'un tableau de booléens: la valeur d'un élément t[i] sera 1 si et seulement si l'indice i est premier. Nous choisirons une représentation par un tableau d'octets: l'octet vaut 0 si l'indice correspondant n'est pas premier, l'octet vaut 1 si l'indice est premier. Pour les timings, nous prendrons M = 1 000 000 (il faut prendre une grande valeur: la taille du tableau doit faire plusieurs fois celle du cache pour prendre en compte la taille limitée de ce dernier).
L'algorithme consiste à initialiser le tableau à 1 (nombres premiers
par défaut). On va barrer
les multiples stricts de 2, 3, 5, 7, ...,
p, ... où p parcourt les nombres premiers inférieurs
à la racine carrée de M. Les nombres non barrés seront les
nombres premiers. Plus précisément, on prend le plus petit nombre premier:
2. On barre alors les multiples de 2 différents de 2 en mettant l'octet
correspondant à 0. On prend ensuite l'entier non barré suivant 2: c'est 3.
On barre alors les multiples de 3 différents de 3. On prend l'entier non
barré suivant 3: 4 est déjà barré (multiple de 2), mais 5 ne l'est pas. Et
ainsi de suite... Bon, si vous ne connaissiez pas, j'espère que vous avez
compris.
Commençons (déjà!) par deux petites optimisations. Nous voulons barrer les multiples stricts de p: 2p, 3p, 4p, ..., kp, ... Si 1 < k < p, alors kp a déjà été barré. On n'a pas besoin de le barrer une deuxième fois! Donc les multiples à barrer sont p.p, (p+1).p, (p+2).p, ... (on a juste une multiplication à effectuer en plus: carré de p). Supposons le nombre premier p impair, le cas p = 2 étant considéré séparément (lors de l'initialisation du tableau). Alors (p+1).p, (p+3).p, (p+5).p, ... sont déjà barrés, car ce sont des nombres pairs. Les multiples à barrer sont maintenant p.p, (p+2).p, (p+4).p, ... en ajoutant 2p à chaque fois.
On peut implémenter cet algorithme; vous pouvez le faire comme exercice. Une implémentation correcte en assembleur devrait prendre environ 140 ms sur notre machine préférée. J'ai aussi implémenté cet algorithme en C et j'ai compilé avec gcc sur une Sun SparcStation4 à 110 MHz, et les calculs se font en... 111 ms! soit moins que sur Risc PC. C'est normal: le crible d'Ératosthène fait des accès intensifs à la mémoire, et ceux-ci ne sont pas du tout locaux. Par conséquent, le cache de niveau 1 et le write buffer (*) du StrongARM ne sont pas suffisants pour empêcher les accès à la mémoire principale de tout ralentir sur le Risc PC.
Nous avons trois méthodes pour être plus rapide sur Risc PC:
Je pense que (3) n'est pas possible, ou alors c'est une conséquence de (1) ou (2): calculs supplémentaires... En revanche, nous allons voir que (1) et (2) sont possibles et peuvent s'appliquer simultanément.
Commençons par (1). Rappelons que nous voulons barrer les multiples de p suivants: p.p, (p+2).p, (p+4).p, ... mais certains sont déjà barrés, et nous voulons éviter de les barrer encore une fois. Avec la deuxième optimisation, nous avons évité de barrer les multiples de 2 une seconde fois; maintenant, nous pouvons éviter de barrer les multiples de 3, mais c'est un peu plus compliqué: en effet, nous devons alternativement ajouter 2p et 4p; nous commencerons par ajouter 4p si p = 3k+1, et 2p si p = 3k+2. On gagne ainsi une écriture en mémoire sur 3. Notons qu'au bout de chaque période, on aura ajouté 2p+4p = 6p. On peut faire de même avec les multiples de 5, ce qui est encore un peu plus compliqué, et cette fois-ci la période est de 2.3.5.p = 30p. L'inconvénient est que l'on ne peut pas aller bien loin: la période augmente très rapidement et on y gagne de moins en moins! Si on va trop loin, on y perd: le code risque de ne plus tenir dans le cache, cela complique les calculs, etc. Nous nous contenterons d'aller jusqu'à 5, même s'il est certainement un peu plus rapide d'éviter aussi de barrer les multiples de 7 une seconde fois. Voilà pour (1).
Passons maintenant à (2), mais faisons d'abord un petit rappel sur le cache du StrongARM (sans entrer dans les détails qui ne nous seront pas utiles), en particulier le cache données. Remarque: ceci ne s'applique pas aux ARM (610, 710, 810...). Le cache est divisé en blocs: les lignes de cache. Pour le StrongARM, une ligne fait 32 octets et le cache fait 16 Ko. La division du cache en lignes induit une division (logique) de la RAM en lignes (de même longueur). Voyons ce qui se passe lors d'une lecture en mémoire (cachable). Si la donnée est déjà dans le cache (cache hit), elle y est directement lue et c'est rapide. Sinon (cache miss), une ligne de cache est éventuellement écrite en RAM pour libérer une ligne et la ligne de la RAM contenant la donnée est entièrement chargée dans le cache, et là c'est très lent, bien plus que s'il n'y avait pas de cache. Prenons maintenant le cas de l'écriture; c'est principalement celui-ci qui nous intéresse. Si l'emplacement mémoire est dans le cache, l'écriture se fait uniquement dans le cache (write back), et la ligne (ou seulement la moitié: il y a un dirty bit pour chaque moitié) sera recopiée en RAM ultérieurement; c'est donc rapide. Dans le cas contraire, l'écriture se fait directement en RAM, comme s'il n'y avait pas de cache; le write buffer est utilisé, mais s'il est plein, c'est lent. Conséquence: quand on veut faire de nombreux accès en écriture au même emplacement (et assez rapprochés dans le temps), il vaut mieux tout d'abord charger la ligne dans le cache en faisant une lecture en mémoire, même si on n'y avait rien écrit auparavant!
Voici maintenant mon algorithme pour que les écritures se fassent dans le cache; il s'agit d'un crible d'Ératosthène par blocs: nous découpons le tableau en blocs de taille un peu inférieure à celle du cache (14 Ko est optimal pour M = 1 000 000). Nous allons d'abord remplir le premier bloc en le chargeant auparavant dans le cache, puis le suivant, et ainsi de suite... Le seul problème est que, pour chaque p, il faut se souvenir où on s'était arrêté dans le bloc précédent; toutes les données utiles sont stockées dans une liste, qui ne prend heureusement pas trop de place; autrement les lectures dans cette liste pourraient enlever du cache une ligne appartenant au bloc, ce qui est mauvais car la ligne ne serait plus du tout rechargée dans le cache. Remarquons que l'algorithme de remplacement de ligne de cache du StrongARM (round robin) évite qu'une telle chose se produise trop souvent, ce qui ne serait pas le cas d'un algorithme aléatoire (ARM 610...). Les interruptions pourraient aussi poser un problème, mais elles sont relativement peu fréquentes.
Enfin, voici mon code assembleur. Je ne sais pas si vous le comprendrez comme ça, surtout qu'il y a quelques petites optimisations et astuces en plus! (et peut-être quelques bugs.) Il ne faut maintenant que 46 ms pour trouver tous les nombres premiers jusqu'à M = 1 000 000. Note: si on ne fait pas de préchargement du bloc dans le cache, cela prend 80 ms.
On pose N = √M, soit ici N = 1000. C'est une routine relogeable. On note π(N) le nombre de nombres premiers inférieurs ou égaux à N. En entrée, R0 doit contenir l'entier N ≥ 5, R1 l'adresse d'un tableau de 8×⌊N²/8⌋+132+8π(N) octets (R1 doit être un multiple de 32), et R2 la taille d'un bloc (multiple de 32 supérieur ou égal à 64, par exemple 14×1024). En sortie, pour 0 ≤ i < N², R1[i] indique si i est premier ou non.
erat STMFD SP!, {R0-R12,LR}
MLA R3, R0, R0, R1 ;R3 pointe sur R1[N^2].
MOV R5, #0
ADR R14, data
LDR R11, [R14], #offsets-data
BIC R4, R3, #7
ADD R4, R4, #120 ;R4 = 8*FLOOR(R3/8)+120.
BIC R3, R3, #1 ;R3 = 2*FLOOR(R3/2).
MOV R7, R1
STMIB R4!, {R1,R5}
; R3: fin du tableau de booléens: R1+N^2 ou R1+N^2-1 si N est impair.
; R4: données temporaires: couples (nombre premier, adresse).
; R5: 0.
; R6: nombre premier p courant.
; R7: fin du bloc courant.
; R11: &11111111.
; R14: tableau d'offsets (longueur: 15 mots).
next_block MOV R8, R7 ;Début du bloc suivant.
ADD R7, R7, R2 ;Fin du bloc...
CMP R7, R3
MOVHI R7, R3 ;ou fin du tableau de booléens.
; Charge le nouveau bloc dans le cache, et supprime (avec un petit temps
; de retard) les multiples de 2, 3 et 5 de ce bloc: stocke 0100 0001 0001
; 0100 0101 0001 0000 0101 0000 0100 0101 0001 0100 0100 0001.
LDR R10, [R4, #-4]
LDR R12, [R10]
LDR R12, [R10, #32]
MOV R9, #&01000000 ;0001.
MOV R6, #&00000100 ;0100.
ORR R8, R9, R6 ;0101.
suppr_235 LDR R12, [R10, #64]
STMIA R10!, {R6,R9}
STR R9, [R10], #4
STMIA R10!, {R6,R8,R9}
STMIA R10!, {R5,R8}
STMIA R10!, {R5,R6,R8,R9}
STR R6, [R10], #4
STMIA R10!, {R6,R9}
LDR R12, [R10, #36]
CMP R10, R7
BCC suppr_235
STR R10, [R4, #-4]
MOV R8, R4
B old_prime
; Supprime les multiples de p dans le bloc courant.
old_loop STRB R5, [R9], R12, LSL #1 ;Supprime p(30k+1).
CMP R9, R7
BCS old_store1
STRB R5, [R9], R6, LSL #2 ;Supprime p(30k+7).
CMP R9, R7
BCS old_store2
STRB R5, [R9], R6, LSL #1 ;Supprime p(30k+11).
CMP R9, R7
BCS old_store3
STRB R5, [R9], R6, LSL #2 ;Supprime p(30k+13).
CMP R9, R7
BCS old_store4
STRB R5, [R9], R6, LSL #1 ;Supprime p(30k+17).
CMP R9, R7
BCS old_store5
STRB R5, [R9], R6, LSL #2 ;Supprime p(30k+19).
CMP R9, R7
BCS old_store6
STRB R5, [R9], R12, LSL #1 ;Supprime p(30k+23).
CMP R9, R7
BCS old_store7
STRB R5, [R9], R6, LSL #1 ;Supprime p(30k+29).
CMP R9, R7
BCC old_loop
old_store0 STMDB R8, {R6,R9}
old_prime LDMIA R8!, {R6,R9} ;Prochain (p,adr) ou R6 = 0.
ADD R10, PC, R6, LSR #25
BICS R6, R6, #&F8000000 ;C = 1.
CMPNE R9, R7 ;Note: R9 (impair) != R7 (pair).
ADDCC R12, R6, R6, LSL #1 ;R12 = 3p.
SUBCC PC, R10, #4*28
BNE old_prime
; Nouveaux nombres premiers...
SUB R8, R8, #8
TEQ R8, R4
LDRNE R6, [R8, #-8] ;R6: dernier nombre premier dans la
BNE last_prime ;liste si celle-ci est non vide.
MOV R6, #7 ;Premier nombre premier: 7.
B new_prime
new_loop STRB R5, [R9], R12, LSL #1 ;Supprime p(30k+1).
CMP R9, R7
BCS new_store1
STRB R5, [R9], R6, LSL #2 ;Supprime p(30k+7).
CMP R9, R7
BCS new_store2
STRB R5, [R9], R6, LSL #1 ;Supprime p(30k+11).
CMP R9, R7
BCS new_store3
STRB R5, [R9], R6, LSL #2 ;Supprime p(30k+13).
CMP R9, R7
BCS new_store4
STRB R5, [R9], R6, LSL #1 ;Supprime p(30k+17).
CMP R9, R7
BCS new_store5
STRB R5, [R9], R6, LSL #2 ;Supprime p(30k+19).
CMP R9, R7
BCS new_store6
STRB R5, [R9], R12, LSL #1 ;Supprime p(30k+23).
CMP R9, R7
BCS new_store7
STRB R5, [R9], R6, LSL #1 ;Supprime p(30k+29).
CMP R9, R7
BCC new_loop
new_store0 STMIA R8!,{R6,R9}
last_prime BIC R6, R6, #&F8000000
next_odd ADD R6, R6, #2 ;Nombre impair suivant.
LDRB R10, [R1, R6]
TEQ R10, #0 ;Boucle tant que ce n'est pas
BEQ next_odd ;un nombre premier p.
new_prime MLA R9, R6, R6, R1 ;R9: pointeur sur le premier
MUL R10, R11, R6 ;multiple à supprimer: p^2.
CMP R9, R7
LDRCCB R10, [R14, R10, LSR #28]
ADDCC R12, R6, R6, LSL #1 ;R12 = 3p.
SUBCC PC, PC, R10, LSL #2
STR R5, [R8] ;Fin de la liste des (p,adr).
TEQ R7, R3
BNE next_block
LDMDB R14, {R11,R12} ;Correction pour les
STMIA R1, {R11,R12} ;entiers de 0 à 7.
LDMFD SP!, {R0-R12,PC}
old_store1 ORR R6, R6, #&18000000
B old_store0
old_store2 ORR R6, R6, #&30000000
B old_store0
old_store3 ORR R6, R6, #&48000000
B old_store0
old_store4 ORR R6, R6, #&60000000
B old_store0
old_store5 ORR R6, R6, #&78000000
B old_store0
old_store6 ORR R6, R6, #&90000000
B old_store0
old_store7 ORR R6, R6, #&A8000000
B old_store0
new_store1 ORR R6, R6, #&18000000
B new_store0
new_store2 ORR R6, R6, #&30000000
B new_store0
new_store3 ORR R6, R6, #&48000000
B new_store0
new_store4 ORR R6, R6, #&60000000
B new_store0
new_store5 ORR R6, R6, #&78000000
B new_store0
new_store6 ORR R6, R6, #&90000000
B new_store0
new_store7 ORR R6, R6, #&A8000000
B new_store0
data DCD &11111111
DCB 0,0,1,1,0,1,0,1
offsets DCB 0,37,25,0,22,0,0,34,19,0,0,31,0,28,16
Timings complets:
Note: article envoyé le 1997-03-01, paru dans ARMada News n°0E. Il s'agit de la version 8 de l'algorithme.
Nous avons déjà une version très rapide du crible d'Ératosthène. La question qui se pose naturellement est: pouvons-nous aller plus vite? Cela revient à peu près à se demander ce qui prend le plus de temps et dans quelles proportions. Rappelons que le crible d'Ératosthène prend actuellement 46 ms.
Considérons maintenant le problème plus simple, consistant à écrire une constante en mémoire (sur un million d'octets). Si nous découpons la mémoire en blocs et préchargeons chaque bloc dans le cache, comme dans le crible d'Ératosthène, cela prend environ 39 ms. Par conséquent, ce sont les transferts entre le cache et la RAM qui prennent la plus grande partie du temps (plus de 80%) dans le crible d'Ératosthène. L'algorithme donné dans l'article précédent peut être encore un peu amélioré en préchargeant convenablement le bloc (accès RAM) en parallèle avec la boucle principale; cela permet de gagner encore 1 ou 2 ms, mais pas beaucoup plus car nous sommes déjà très proches des 39 ms.
Maintenant cherchons la méthode la plus rapide d'écrire une constante en mémoire, ce qui nous donnera une borne inférieure sur la durée du crible d'Ératosthène. L'écriture par STM avec 4 ou 8 registres (avec alignement sur une adresse multiple de 16 ou 32) semble optimale; elle prend 24 ms. Notons qu'en remplaçant les STM par des STR, i.e. en écrivant 4 octets à la fois au lieu de 16, le temps passe à 40 ms; ceci doit être dû au fait que le bus du Risc PC n'est pas en enhanced mode.
Contrairement à la version avec préchargement du cache, nous ne faisons pas de lecture en mémoire principale. C'est pour cela que c'est plus rapide. La seule façon d'appliquer cette méthode au crible d'Ératosthène est de stocker tous les résultats intermédiaires dans un bloc fixe en mémoire, appelé bloc image car il représente un bloc du tableau situé en RAM: le bloc image est préchargé (accès RAM) tout au début, et ensuite tous les STR se font dans ce bloc, donc dans le cache. Lorsque le bloc a fini d'être calculé, il est recopié en RAM par LDM/STM; les LDM se font dans le cache, et les STM dans la RAM via le write buffer. Il faut faire attention en choisissant la taille du cache, car si une partie du bloc image est retirée du cache (notamment à cause des interruptions), elle ne pourra plus être rechargée dans le cache par la suite. Pour éviter ce problème, nous pourrions recharger de temps en temps le bloc image dans le cache, mais cela n'est apparemment pas nécessaire même avec de plus grandes valeurs de N.
Pour N = 1000 (recherche des nombres premiers inférieurs à 1 000 000), avec une taille de blocs de 9600 octets, cela prend 30,8 ms... à comparer avec les 24 ms de l'écriture d'une constante en mémoire par STM. Notons aussi que cette implémentation du crible d'Ératosthène est plus rapide qu'une simple copie de bloc d'un million d'octets en mémoire.
Le code se trouve sur http://www.ens-lyon.fr/~vlefevre/acorn/erat8.fr.html [maintenant http://www.vinc17.org/acorn/erat8.fr.html].