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

Introduction

Voici une implémentation en ARM du crible d'Ératosthène, spécialement optimisée pour le Risc PC avec un processeur StrongARM (202 MHz). Avec cette configuration, la bande passante est très faible par rapport à la rapidité du processeur. Nous cherchons donc à:

diminuer le nombre d'écritures en mémoire (cache ou RAM) en évitant des écritures redondantes,
utiliser le cache au maximum en découpant le tableau de booléens en blocs (tenant dans le cache) et en traitant ces blocs les uns après les autres,
utiliser le write buffer au maximum en évitant de faire des accès RAM trop rapprochés dans le temps.

Un article avec plus de détails paraîtra dans le numéro 0C d'ARMada News (revue éditée par l'ARMada, une association française des utilisateurs de machines Acorn).

Routines

La routine assembleur se trouve ci-dessous. Vous pouvez l'assembler avec !ObjAsm, après avoir éventuellement remplacé le nom de fichier h:RegNames. Pour l'utiliser, vous pouvez prendre le petit programme C qui se trouve plus loin. Sinon voici comment l'appel se fait pour rechercher les nombres premiers de 0 à N²-1: le registre R0 doit contenir l'entier N (supérieur ou égal à 5), le registre R1 doit contenir l'adresse d'un tableau de 4.⌊N²/4⌋ + 128 + π(N) octets, où π(N) est le nombre de nombres premiers inférieurs ou égaux à N, et le registre R2 doit contenir la taille des blocs (multiple de 32 supérieur ou égal à 64). Au retour, les registres ne sont pas modifiés et les N² premiers octets du tableau indiquent si le nombre correspondant est premier (1) ou non (0).

Note: le format de sortie (tableau d'octets) a été choisi et fixé à l'avance. On aurait pu avoir une routine plus rapide si un autre format de sortie, utilisant moins de mémoire (tableau de bits et/ou non-stockage des nombres pairs), avait été choisi. Cependant, cela change le problème, qui est de trouver une optimisation pour un format de sortie fixé. De toute façon, les mêmes techniques seraient utilisées.

; Crible d'Eratosthene, version 7 (1997-02-09)
; Recherche des nombres premiers jusqu'à N^2-1.

; Entrée:
; R0: entier N >= 5.
; R1: tableau de 4*[N^2/4]+128+8*pi(N) octets. Adresse multiple de 32.
; R2: taille d'un bloc (multiple de 32, >= 64).

; Sortie:
; Registres non modifiés.
; Résultats dans R1[0..N^2-1]: 0 (non premier) ou 1 (premier).

		GET	h:RegNames	;Définit R0, ..., R14, SP, LR.

		AREA	area_erat, READONLY, CODE, PIC
		EXPORT	erat

erat		STMFD	SP!, {R0,R3-R12,LR}
		MLA	R3, R0, R0, R1	;R3 = R1+N^2.
		MOV	R5, #0		;R5 = 0.
		ADR	R14, offsets
		BIC	R3, R3, #1	;R3 = 2*FLOOR(R3/2).
		ADD	R4, R3, #116	;R4 = 4*FLOOR(R3/4)+116.
		MOV	R7, R1
		STMIB	R4!, {R1,R5}

; Charge le premier bloc.

		ADD	R8, R7, R2	;Fin du premier bloc...
		CMP	R8, R3
		MOVHI	R8, R3		;ou fin du tableau de booléens.
		MOV	R0, R1
load_first	LDRB	R11, [R0], #32
		CMP	R0, R8
		BCC	load_first

; Dans la suite:
; R1: tableau.
; R2: taille d'un bloc.
; R3: R1+N^2 (N pair) ou R1+N^2-1 (N impair).
; R4: données temporaires: couples (nombre premier, adresse).
; R5: 0.
; R6: nombre premier p courant.
; R7: fin du bloc courant.
; 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.

; Initialisation du bloc et suppression des multiples de 2, 3 et 5
; (on dépasse au plus de 116 octets): on stocke 0100 0001 0001 0100
; 0101 0001 0000 0101 0000 0100 0101 0001 0100 0100 0001.

		LDR	R10, [R4, #-4]
		MOV	R9, #&01000000	;0001.
		MOV	R6, #&00000100	;0100.
		ORR	R8, R9, R6	;0101.
suppr_235	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}
		CMP	R10, R7
		BCC	suppr_235
		STR	R10, [R4, #-4]
		MOV	R8, R4
		ADD	R0, R7, R2	;Fin du bloc...
		CMP	R0, R3		;ou fin du tableau de booléens
		BICHI	R0, R3, #&1F	;(multiple de 32 inférieur).
		B	old1_prime

; Charge le bloc suivant et supprime les multiples de p
; dans le bloc courant.

old1_loop	ADDS	R11, R11, #&80000000	;LDR: une fois sur 2.
		CMPCS	R0, R7
		LDRCSB	R11, [R0], #-32		;Charge le bloc suivant...
		STRB	R5, [R9], R12, LSL #1	;Supprime p(30k+1).
		CMP	R9, R7
		BCS	old1_store1
		STRB	R5, [R9], R6, LSL #2	;Supprime p(30k+7).
		CMP	R9, R7
		BCS	old1_store2
		STRB	R5, [R9], R6, LSL #1	;Supprime p(30k+11).
		CMP	R9, R7
		BCS	old1_store3
		STRB	R5, [R9], R6, LSL #2	;Supprime p(30k+13).
		CMP	R9, R7
		BCS	old1_store4
		STRB	R5, [R9], R6, LSL #1	;Supprime p(30k+17).
		CMP	R9, R7
		BCS	old1_store5
		STRB	R5, [R9], R6, LSL #2	;Supprime p(30k+19).
		CMP	R9, R7
		BCS	old1_store6
		STRB	R5, [R9], R12, LSL #1	;Supprime p(30k+23).
		CMP	R9, R7
		BCS	old1_store7
		STRB	R5, [R9], R6, LSL #1	;Supprime p(30k+29).
		CMP	R9, R7
		BCC	old1_loop		;Boucle.
old1_store0	STMDB	R8, {R6,R9}	;Stocke (p,adr).
		CMP	R0, R7		;Branchement si le chargement
		BCC	old2_prime	;du bloc suivant est terminé.
old1_prime	LDMIA	R8!, {R6,R9}	;Prochain (p,adr) ou R6 = 0.
		ADD	R10, PC, R6, LSR #25
		BICS	R6, R6, #&F8000000	;Note: C = 1.
		CMPNE	R9, R7		;Note: R9 (impair) != R7 (pair).
		ADDCC	R12, R6, R6, LSL #1	;R12 = 3p.
		SUBCC	PC, R10, #4*30	;Branchement au bon STRB.
		BNE	old1_prime	;Branchement si R6 != 0.
load_block	CMP	R0, R7		;Finit de charger le bloc suivant.
		LDRCSB	R11, [R0], #-32
		BCS	load_block
		B	first_new

old2_loop	STRB	R5, [R9], R12, LSL #1	;Supprime p(30k+1).
		CMP	R9, R7
		BCS	old2_store1
		STRB	R5, [R9], R6, LSL #2	;Supprime p(30k+7).
		CMP	R9, R7
		BCS	old2_store2
		STRB	R5, [R9], R6, LSL #1	;Supprime p(30k+11).
		CMP	R9, R7
		BCS	old2_store3
		STRB	R5, [R9], R6, LSL #2	;Supprime p(30k+13).
		CMP	R9, R7
		BCS	old2_store4
		STRB	R5, [R9], R6, LSL #1	;Supprime p(30k+17).
		CMP	R9, R7
		BCS	old2_store5
		STRB	R5, [R9], R6, LSL #2	;Supprime p(30k+19).
		CMP	R9, R7
		BCS	old2_store6
		STRB	R5, [R9], R12, LSL #1	;Supprime p(30k+23).
		CMP	R9, R7
		BCS	old2_store7
		STRB	R5, [R9], R6, LSL #1	;Supprime p(30k+29).
		CMP	R9, R7
		BCC	old2_loop		;Boucle.
old2_store0	STMDB	R8, {R6,R9}	;Stocke (p,adr).
old2_prime	LDMIA	R8!, {R6,R9}	;Prochain (p,adr) ou R6 = 0.
		ADD	R10, PC, R6, LSR #25
		BICS	R6, R6, #&F8000000	;Note: C = 1.
		CMPNE	R9, R7		;Note: R9 (impair) != R7 (pair).
		ADDCC	R12, R6, R6, LSL #1	;R12 = 3p.
		SUBCC	PC, R10, #4*28	;Branchement au bon STRB.
		BNE	old2_prime	;Branchement si R6 != 0.

; Nouveaux nombres premiers.

first_new	LDR	R11, [R14, #data-offsets] ;&11111111.
		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		;Boucle.
new_store0	STMIA	R8!,{R6,R9}	;Stocke (p,adr).

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,R3-R12,PC}

old1_store1	ORR	R6, R6, #&18000000
		B	old1_store0
old1_store2	ORR	R6, R6, #&30000000
		B	old1_store0
old1_store3	ORR	R6, R6, #&48000000
		B	old1_store0
old1_store4	ORR	R6, R6, #&60000000
		B	old1_store0
old1_store5	ORR	R6, R6, #&78000000
		B	old1_store0
old1_store6	ORR	R6, R6, #&90000000
		B	old1_store0
old1_store7	ORR	R6, R6, #&A8000000
		B	old1_store0

old2_store1	ORR	R6, R6, #&18000000
		B	old2_store0
old2_store2	ORR	R6, R6, #&30000000
		B	old2_store0
old2_store3	ORR	R6, R6, #&48000000
		B	old2_store0
old2_store4	ORR	R6, R6, #&60000000
		B	old2_store0
old2_store5	ORR	R6, R6, #&78000000
		B	old2_store0
old2_store6	ORR	R6, R6, #&90000000
		B	old2_store0
old2_store7	ORR	R6, R6, #&A8000000
		B	old2_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

		END

Le programme suivant prend 3 arguments:

le nombre N (pour rechercher les nombres premiers de 0 à N²-1),
la taille des blocs (multiple de 32),
un entier nul pour afficher les nombres premiers, ou un entier positif pour avoir le temps d'exécution (c'est alors le nombre d'itérations, plus il est grand, plus c'est précis).

La taille des blocs optimale devrait être environ un peu inférieure à C/2 - 4N/Ln(N), où C est la taille du cache de données. Lorsque N devient très grand, les blocs peuvent devenir trop petits; mais la valeur de N est d'abord limitée par la quantité de RAM.

/* Crible d'Eratosthene - test et timing
 *
 * Usage: erat <N> <taille_bloc> <p>
 * --> Recherche des nombres premiers jusqu'à N^2-1
 * p = 0: affichage
 * p > 0: benchmark (nombre d'itérations)
 */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

void erat(int, char *, int);

int main(int argc, char **argv)
{
  int i,n,s,p;
  char *t;
  clock_t start,stop;

  if (argc != 4) exit(1);
  if ((n = atoi(argv[1])) < 5) exit(2);
  if ((s = atoi(argv[2])) < 64 || s%32) exit(3);
  if ((p = atoi(argv[3])) < 0) exit(4);
  if (!(t = malloc(n*(n+8)+160))) exit(5);
  t = (char *) (((int) t + 31) & ~31);
  if (p) {
      start = clock();
      for(i=0;i<p;i++) erat(n,t,s);
      stop = clock();
      printf("Temps = %.3f ms\n",
             (double) (stop-start) * (1000.0/CLK_TCK) / p);
    }
  else {
      erat(n,t,s);
      for(i=0;i<n*n;i++)
        if (t[i]) printf("%d\n", i);
    }
  return 0;
}

Résultats

Pour N = 1000 (recherche des nombres premiers inférieurs à 1 000 000), avec une taille de blocs de 7168 octets (7 Ko), cela prend 43,3 ms.

Allez voir aussi une version plus rapide: la version 8.

Dernière modification:

Source XML: 2003-06-18, 23:39:56 (20)
Feuille de styles XSL: 2003-08-06, 22:11:54 (488)

webmaster@vinc17.org