English version(Document écrit en août 1995.)
Armpi 3 est un programme de calcul de π en décimal pour les processeurs ARM (en particulier pour le Risc PC d'Acorn). On utilisera la formule suivante: π/4 = arctan(1/2) + arctan(1/3) avec arctan x = x - x3/3 + x5/5 - x7/7 + ... car elle est simple et converge assez rapidement, bien que l'on préfère souvent d'autres formules plus rapides. C'est cette formule que j'ai aussi utilisée sur l'Apple II (6502) et l'Atari ST (68000); vous pouvez voir les comparaisons.
Les résultats seront donnés à la fin.
Le programme donné ici est une routine standard dont l'adresse de retour doit se trouver dans le registre LR (Link Register). Lors de l'appel, le registre R0 contient le nombre de chiffres voulu aux erreurs d'arrondi près, et le registre R1 contient l'adresse de la zone mémoire réservée pour les calculs et le stockage du résultat. Cette zone mémoire doit contenir au moins 2N + B + 1544 octets (cf ci-dessous), où B est le plus petit multiple de 4 supérieur ou égal à 5N/12. R0 et R1 doivent être des multiples de 4.
Au retour, les registres R0, R1 et R13 ne sont pas modifiés, et le résultat se trouve à l'adresse donnée par R1, poids faible en premier, avec un chiffre décimal par octet.
Note: l'ARM doit être en configuration little-endian.
Les calculs vont être effectués en binaire. Le résultat final sera ensuite converti en décimal. Pour N chiffres décimaux, il faut un peu moins de 5N/12 octets (1 octet = 8 chiffres binaires). On va donc faire les calculs sur B octets, où B est le plus petit multiple de 4 supérieur ou égal à 5N/12 (il vaut mieux choisir trop grand que trop petit).
On va regrouper les divisions par un même nombre, i.e. au lieu de calculer et additionner (1/2n)/n et (1/3n)/n, on va calculer (1/2n+1/3n)/n. Le calcul de 1/2n consiste en un décalage; en fait, on ajoutera un 1 à 1/3n, ce qui prendra un temps constant négligeable. Le calcul de 1/3n consiste en une division par 9 à chaque itération; on choisira un algorithme spécial de division par 9, très rapide: l'idée est de multiplier par l'inverse de 9 dans Z/232Z. La division par n se fera en calculant des tables et en lisant dans ces tables. Les additions et les soustractions se feront de manière classique. Enfin, il reste la conversion en base 10, qui prend un temps relativement long; l'algorithme général est classique, mais on fera les calculs sur 4 chiffres décimaux en parallèle (on a 4 chiffres par mot de 32 bits), avec représentation redondante (un chiffre peut être supérieur à 9).
Les sections suivantes contiennent les explications des différentes parties du programme. Le source assembleur est donné dans un fichier séparé. Voici les différentes parties:
Les routines données avant la partie principale seront utilisées
comme des macros, mais la syntaxe est ici différente: pour plus de
clarté, on ne mettra pas le $
pour les paramètres. Des scripts
Perl permettent d'extraire du fichier donné
ici le source assembleur conforme à la syntaxe
standard ou le source au format extASM.
Un autre script Perl crée un
programme BASIC (pour le Risc PC) à partir du
source assembleur renvoyé par le premier script (syntaxe standard).
On veut effectuer une division d'un grand nombre par 9, le résultat étant stocké à la place de ce nombre. On va effectuer une succession de divisions de nombres N = (NA,NB) de 36 bits (dont les 4 bits de poids fort NA forment un nombre inférieur à 9) par 9, avec quotient sur 32 bits et reste sur 4 bits.
D'abord, on se débarrasse des 4 bits de poids fort en soustrayant NA × &FFFFFFFC, puis en ajoutant 36 si le résultat est négatif. On obtient ainsi un quotient partiel, égal à NA.a (+ 4 éventuellement) avec a = &FFFFFFFC / 9 = &1C71C71C, le reste étant la nouvelle valeur de NB. Maintenant, il reste à diviser un nombre de 32 bits par 9.
L'inverse de 9 dans Z/232Z est 2a+1. Soit R le résultat sur 32 bits de N.(2a+1). La table suivante donne la valeur de N modulo 9 en fonction de l'intervalle contenant R:
[0,a]......... N % 9 = 0 [a+1,2a]...... N % 9 = 5 [2a+1,3a+1]... N % 9 = 1 [3a+2,4a+1]... N % 9 = 6 [4a+2,5a+2]... N % 9 = 2 [5a+3,6a+2]... N % 9 = 7 [6a+3,7a+3]... N % 9 = 3 [7a+4,8a+3]... N % 9 = 8 [8a+4,9a+3]... N % 9 = 4
On multiplie NB par 2a+1, qui s'écrit en base 2: 00111000111000111000111000111001 et qui peut aussi être écrit en chiffres signés: 0100m00100m00100m00100m00100m001 où m = -1. Ainsi la multiplication peut être effectuée en utilisant 1 SUB et 3 ADDs. Ensuite, on détermine le reste de manière à ramener le résultat de la multiplication dans [0,a].
Il s'agit ici de diviser un grand nombre par un nombre de 32 bits. Le quotient sera stocké dans un autre bloc mémoire (contrairement à la division par 9). L'algorithme donné ici ne marche qu'avec un diviseur dont le MSB (bit de poids fort) est nul; le diviseur est donc en fait sur 31 bits (ce qui correspond aussi aux diviseurs positifs en arithmétique signée), mais cela est largement suffisant pour les calculs que l'on veut effectuer.
La division longue s'effectue de manière classique en base 232. À chaque itération, on connaît le reste de la division effectuée à l'itération précédente (mot de 32 bits inférieur au diviseur), initialement nul. On doit diviser le nombre de 64 bits, appelé dividende partiel, formé du reste et du chiffre suivant du dividende (mot de 32 bits). On obtient un chiffre du quotient (sur 32 bits) et le reste, qui sera utilisé à l'itération suivante. Le dividende partiel et le diviseur seront normalisés, i.e. décalés d'un même nombre de bits vers la gauche de telle sorte que le MSB du diviseur soit à 1; comme le dividende et le diviseur sont multipliés par un même nombre, le quotient est inchangé et le reste est automatiquement normalisé lors du calcul.
Maintenant, expliquons un peu en détail la division du dividende partiel X = (NA,NB) sur 64 bits (2 mots) par le diviseur DV sur 32 bits (dont le MSB est à 1). La division va se faire en 4 étapes identiques, à l'aide de tables calculées avant la division globale (ces tables dépendent uniquement du diviseur). On a par hypothèse NA < DV.
En considérant les 9 bits de poids fort du dividende partiel normalisé, on obtient, en lisant dans une table à 512 entrées (au maximum), 8 bits du quotient, formant un nombre NQ appelé quotient partiel. En fait, il peut généralement y avoir 2 quotients partiels possibles: NQ et NQ+1; pour avoir un quotient unique, il faudrait considérer plus de bits (jusqu'à 38), ce qui est trop pour une table. On considère donc que c'est NQ (comme si les bits suivants les 9 bits de poids fort du dividende étaient tous nuls), et on fera une correction plus tard si c'était NQ+1. Ensuite, il faut soustraire (NQ<<25).DV de X, en lisant NQ.DV dans une table à 256 entrées (chaque entrée étant un quotient partiel possible), et décaler le résultat de 8 bits vers la gauche. NQ.DV est un nombre de 39 bits (NQ a 8 bits et DV a 31 bits). Mais si NQ est le quotient réel, on sait que les 8 bits de poids fort de X-NQ.DV sont nuls, et dans ce cas, on n'a pas besoin de représenter les 8 bits de poids fort de NQ.DV; on a seulement besoin du bit de poids faible de ces 8 bits pour déterminer le quotient. On ne stockera donc que 32 bits du produit NQ.DV, ce qui forme juste un mot P (cet algorithme n'est donc pas valable pour un diviseur de 32 bits). On calcule ainsi X' = ((X<<7)-P)<<1. Si le bit perdu lors du dernier décalage est non nul ou si X' ≥ DV, alors le quotient réel est NQ+1, et non NQ; on fait la correction éventuelle du quotient et de X' maintenant. On recommence ces calculs 3 fois. Au bout de 4 étapes, on obtient les 32 bits du quotient (que l'on stocke en mémoire) et le reste normalisé.
Cette table est calculée de la manière suivante, où DV est le diviseur normalisé:
NA contient les 9 bits de poids fort des produits NQ.DV. On stocke entre 256 et 512 valeurs.
Cette table est calculée de la manière suivante, où DV est le diviseur normalisé:
Note: puisque le diviseur est sur 31 bits, le LSB de DV est nul. Dans les additions, les bits de poids fort (jusqu'à 7 bits) sont perdus (cf Division longue par n).
On a un nombre binaire, défini à une puissance de 2 près, à convertir en base 10: π (on a en fait calculé π/4, mais 4 est une puissance de 2). 3 variables (nombres longs) vont être utilisées dans la conversion: le nombre binaire, un nombre décimal T égal à une puissance de 2 (l'exposant est décrémenté à chaque itération, et devient négatif), et le résultat temporaire R en base 10 (de plus en plus précis). Les deux nombres décimaux ont le même nombre de chiffres (zéros non significatifs compris) et la valeur du chiffre de poids faible doit être (légèrement) supérieure à celle du chiffre de poids faible du nombre binaire.
Comme le bit de poids fort du nombre binaire π vaut 2, la première puissance de 2 considérée est 2. La valeur initiale de R est évidemment 0. On parcourt successivement les bits du nombre binaire. À chaque itération, on ajoute T à R si le bit vaut 1; puis on divise T par 2. On itère jusqu'à ce que T soit nul (par manque de précision).
Pour les nombres décimaux, on a un chiffre par octet. Comme les nombres binaires, les nombres décimaux sont stockés poids faible en premier: ceci est dû au fait que l'ARM est little-endian par défaut.
Pour optimiser, les nombres T et R sont entrelacés par mots de 32 bits; les mots de T sont situés avant les mots de R correspondants. Il s'agit d'une optimisation spécifique à l'ARM, qui permet de lire (ou d'écrire) plusieurs mots consécutifs à l'aide d'une seule instruction. Aussi, on a besoin de moins de registres.
On pourrait effectuer les opérations de division par 2 et d'addition comme on peut le faire à la main (i.e. chiffre par chiffre, en s'assurant que chaque chiffre reste entre 0 et 9). Mais on peut accélérer très fortement le calcul, d'une part en faisant des opérations élémentaires sur 32 bits au lieu de 8, i.e. en considérant 4 chiffres à la fois, d'autre part en utilisant une représentation redondante, i.e. des chiffres pouvant être supérieurs à 9 (mais ne dépassant jamais 255), ce qui évite de calculer les retenues à chaque addition.
La division par 2 s'effectue mot par mot en partant du poids fort. À chaque itération, on effectue une rotation du mot vers la droite avec retenue; on obtient un mot M et une nouvelle retenue. Les bits 7, 15, 23 et 31 contiennent les retenues affectant chacun des 4 chiffres. On écrit donc M = N+C, où C contient les retenues et N les autres bits. On ajoute C>>5 + C>>7 = 5(C>>7) à N pour obtenir le résultat voulu. Note: ces calculs en parallèle sont possibles grâce au fait que la base (10) est divisible par le diviseur (2).
Pour qu'elle soit très rapide, l'addition s'effectue mot par mot
en partant des mots de poids faible, et sans retenue. Mais il faut
régulièrement diminuer la valeur des chiffres du résultat pour éviter
des dépassements de capacité. Ce nettoyage
se fait toutes les
8 itérations, i.e. à chaque fois qu'un octet du nombre
binaire a été entièrement lu, après l'addition éventuelle. On parcourt
alors un par un les chiffres du résultat, en partant du chiffre de
poids faible; à chaque fois qu'un chiffre + la retenue est supérieur
ou égal à 128, on le diminue de 80 et on génère une retenue égale à 8.
Ainsi les chiffres sont toujours compris entre 0 et 199, et il n'y a
jamais de dépassement de capacité.
À la fin de la conversion, les chiffres sont compris entre 0 et 190.
Alors on normalise
le résultat: on obtient un nombre ayant des
chiffres compris entre 0 et 9. On recopie en même temps le nombre
normalisé où il faut (les mots ne sont plus entrelacés).
Dans la première partie (avant la conversion), 3 zones mémoires sont utilisées: une contenant les nombres de la forme 4/32n+1, une contenant le quotient dans la division longue par 2n+1, et une contenant le résultat provisoire.
D'abord, on calcule B de la façon suivante: B = ((N + N/4)/3 + 3) BIC 3. On calcule les adresses des différentes zones utilisées, on initialise le diviseur et le décalage pour le normaliser. On stocke 4 (1/2 + 1/3) = 1101010101... dans la zone qui contiendra le résultat en binaire (le bit de poids fort valant 2), et 4/3 = 01010101... dans la zone qui contient les puissances négatives de 3. Ensuite, les calculs peuvent commencer.
À chaque itération, on effectue la division par 9, on calcule la table des quotients partiels et la table des produits pour la division par 2n+1; on met un bit à 1, correspondant au 4/22n+1, et on passe à la conversion si ce bit se situe en dehors de la zone, i.e. si la précision demandée a été atteinte. Ensuite, on effectue la division par 2n+1, et on ajoute ou soustrait le quotient au/du résultat suivant la parité de n. On remet le bit à 0. On ajoute 2 au diviseur, et on enlève 1 au décalage si le diviseur a dépassé une nouvelle puissance de 2. Puis on boucle.
Pour la conversion, on considère de nouvelles zones mémoires. On initialise les nombres décimaux et on effectue la conversion comme indiqué ci-dessus.
200 000 décimales ont été calculées en 7 heures 40 minutes 7 secondes: le calcul en binaire a pris 3 heures 47 minutes 15 secondes et la conversion en base 10 a pris 3 heures 52 minutes 51 secondes.
Les 3 temps de calcul sont à peu près proportionnels au carré du nombre de décimales calculées. Avec 20 000 décimales, les temps sont les suivants: 4 minutes 36 secondes au total, 2 minutes 22 secondes pour le calcul en binaire et 2 minutes 14 secondes pour la conversion. Les temps ont été respectivement divisés par 100, 96, 104.
Le calcul en binaire et la conversion en base 10 prennent à peu près
le même temps. Il peut donc être intéressant d'optimiser l'un ou l'autre,
mais si on n'optimise pas les deux, l'accélération restera limitée. Pour le
calcul en binaire, on peut prendre une autre formule (toujours sous forme
d'une somme d'atan). Pour la conversion, c'est plus dur: on peut effectuer
le nettoyage
un peu moins souvent (toutes les 8 additions au lieu de
tous les 8 bits lus) ou faire une division par 4 au lieu de 2 divisions par
2 quand un bit est nul et qu'il n'y a pas d'addition à faire; cf
Armpi 4.
Il vaudrait peut-être mieux effectuer tous les calculs en une base
de la forme 10n, où n est bien choisi.
Les comparaisons avec les autres ordinateurs et autres programmes de calcul de π que j'ai écrits sont données dans un fichier séparé; tous les programmes jusqu'à Armpi 4 utilisent la même formule, mais les algorithmes de calculs sont assez différents.
