Mail/M'écrire



Implémentation d'un algorithme simple.
(André Brouty janvier 1996)
Cet exposé est sous licence libre LLDD.

La formule choisie est classique

PI = 16arctg(1/5) - 4arctg(1/239))
il convient de calculer les arcs tangentes. Ceux-ci se calculent a l'aide de la série de Taylor bien connue de terme générale:
Ui(X)=(-1)i/((2*i+1)*X(2*i+1))
qui elle-même s'obtient pas récurence suivant les formules:
V0 = 1/X
Vn+1 = [(-1)*(2*n+1)/((2*n+3)*X*X)]*Vn

La difficulté de l'implémentation de ces formules vient qu'il faut faire les calculs en multiprécision. On aura besoin de faire des additions et des divisions, éventuellement des multiplications.

Les additions en multiprécision ne posent pas de problème particulier, car la somme de deux nombres de p et q chiffres est un nombre d'au plus max(p,q) + 1 chiffres. Par contre pour la multiplication il en va autrement car le produit de deux nombres de p et q chiffres donne un nombre de p+q-2, p+q-1 ou p+q chiffres et la complexité de l'algorithme n'est plus linéaire. Comme on a décidé de faire simple on ne multipliera qu'un grand nombre par un petit, le petit tenant dans un mot mémoire de la machine où se fera le calcul. De même pour la division. Ce choix a l'avantage de la simplicité et de la performance mais va imposer des limites au nombre de décimales que l'on pourra calculer. On va donc représenter un grand nombre par un tableau d'entiers, chacun ayant une représentation interne sur la machine de b bits. Par exemple un tableau de cent entiers permet de représenter un nombre sur 100b bits. En théorie seulement car de nouvelles contraintes vont se poser. On sera donc amené à multiplier un nombre sur 100b bits par un nombre sur b bits en multipliant le nombre de b bits par chaque case du tableau. Le problème est qu'en faisant cette opération il ne faut pas obtenir un nombre supérieur à b bits qu'on ne saura pas stocker. Le nombre placé dans chaque case du tableau ne devra pas dépasser une valeur maximale afin d'éviter ce problème. Evaluons la valeur maximale autorisée.

Examinons le cas de la multiplication:

On fait comme à l'école primaire mais sur chaque case du tableau. La case du tableau représentant un chiffre sur b bits, dans une base dont la valeur sera notée BASE. Chaque case du tableau C est inférieure à BASE

C[i] <= BASE - 1
C[i-1] = C[i]*N + R
(R la retenue de la multiplication de la case précédente).
R = C[i+1]*N/BASE donc
R <= N*(BASE - 1)/BASE et ainsi R < N donc
C[i]*N + R < N*(BASE - 1) + N <= 2b - 1
On obtient donc la relation:
N <= (2b - 1)/BASE (1)
Il faudra pour la multiplication choisir l'entier N vérifiant cette contrainte.

Examinons le cas de la division:

De même on divisera un grand nombre représenté par un tableau par un nombre ordinaire N. Alors comme à l'école primaire on a

C[i] = q*N + R ( R le reste, R < N)
puis on abaisse la tranche suivante et on obtient
R*BASE + C[i+1] à diviser par N. Or
R <= N - 1 et C[i+1] < (N - 1)*BASE + BASE d'où:
N*BASE < 2b - 1
N < (2b - 1)/BASE
on retrouve la même inégalité que pour la multiplication.

le nombre de bits b étant fixé par le type de machine, il nous reste à choisir la BASE de numération qu'on a intéret à prendre la plus grande possible afin d'avoir le moins de calculs à faire. Plus cette valeur de BASE est grande plus le nombre N sera petit.

La valeur de BASE étant choisie, il nous faut déterminer combien de décimales on peut calculer avec ces contraintes.

Un nombre X > 1 étant donné, le nombre de zéros après la virgule de la représentation en décimale de 1/X est donné par E[log(X)] si X n'est pas une puissance de 10 et par log(X) - 1 sinon, le log étant en base 10.
Cette remarque est importante car le calcul des décimales de PI s'arrêtera quand toutes les cases du tableau représentant ces décimales seront nulles.
Compte tenues des contraintes de la formule (1) combien de décimales pourra-t-on obtenir de chaque arctg(1/X) ? Autrement dit combien aura-t-on de zéros pour le terme de la série

Ui(X)=(-1)i/((2*i+1)*X(2*i+1))
avec l'indice i maximal vérifiant (1) ? Le nombre de zéro de Ui(X) est donné par
E[log(2*i+1) + (2*i+1)*log(X)];
X peut valoir, 5 ou 239 pour la formule de Machin. Les tableaux ci-dessous donnent les résultats pour les différentes valeurs de b et de BASE.

Les machines à mots de 32 bits

Nous avons choisi d'utiliser des entiers signés, donc notre mot fait 31 bits et la relation (1) nous impose:
N < 2147483647/BASE
BASE ET INDICESDECIMALES ATTEINTES POUR arctg(1/X)
BASE N MAX i=(N-3)/2X=5X=8X=239
1 000 00021471072150219405104
100 0002147410735150121939551070
10 000214748107372150106193940510756

Les machines à mots de 64 bits

Nous avons choisi d'utiliser des entiers signés, donc notre mot fait 63 bits et la relation (1) nous impose:
N < 9223372036854775807/BASE
BASE ET INDICESDECIMALES ATTEINTES POUR arctg(1/X)
BASE N MAX i=(N-3)/2X=5X=8X=239
10159 2234 6106 4498 33121 935
101492 23346 11564 47183 297219 366
1013922 337461 167644 690832 9572 193 685

On est donc amené à changer de base suivant le nombre de décimales désirées, ce qui conduit à une baisse de performance pour les grands calculs.

Performances

ce programme a les performances suivantes:

SUR SUN SPARC 10, 32 BITSSUR DEC ALPHA OSF, 64 BITS
DECIMALESTEMPSTEMPS
1 5001.8 s0.5s
10 0001mn 32s 27.3s
500 000 - 19h 40mn

Le temps de calcul pour 500 000 décimales sur ce type de machine n'est pas très bon.

Les machines évoluant très vite on trouvera ici les temps de calcul sur les diverses machines du marchés permettant ainsi de suivre l'évolution de leur performance.