Crible d'Ératosthène sur Risc PC avec StrongARM

Partie 1

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:

  1. Améliorer l'algorithme en diminuant le nombre d'accès mémoire.
  2. Essayer d'utiliser le cache de niveau 1 en localisant les accès mémoire.
  3. Essayer d'utiliser le write buffer en évitant de le remplir.

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:

(*) Write buffer
c'est un système permettant au processeur de ne pas se bloquer sur un STR (écriture en mémoire), qui peut prendre du temps, et de continuer à exécuter les instructions suivantes pendant l'écriture en mémoire. Plusieurs écritures peuvent être mises en attentes, suivant la taille du write buffer (pour le StrongARM: 8 blocs de 1 à 16 octets).

Partie 2

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].



Valid XHTML 1.0!
Dernière modification:
webmaster@vinc17.org