PI: implémentation

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:

U_i(X)=(-1)ⁱ/((2*i+1)*X^(2*i+1))

qui elle-même s'obtient pas récurence suivant les formules:

V₀ = 1/X
V_n+1 = [(-1)*(2*n+1)/((2*n+3)*X*X)]*V_n

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 <= 2^b - 1

On obtient donc la relation:

N <= (2^b - 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 < 2^b - 1
N < (2^b - 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)ⁱ/((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 INDICES			DECIMALES ATTEINTES POUR arctg(1/X)
BASE	N MAX	i=(N-3)/2	X=5	X=8	X=239
1 000 000	2147	1072	1502	1940	5104
100 000	21474	10735	15012	19395	51070
10 000	214748	107372	150106	193940	510756

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 INDICES			DECIMALES ATTEINTES POUR arctg(1/X)
BASE	N MAX	i=(N-3)/2	X=5	X=8	X=239
10¹⁵	9 223	4 610	6 449	8 331	21 935
10¹⁴	92 233	46 115	64 471	83 297	219 366
10¹³	922 337	461 167	644 690	832 957	2 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 BITS SUR DEC ALPHA OSF, 64 BITS

DECIMALES TEMPS TEMPS

1 500 1.8 s 0.5s

10 000 1mn 32s 27.3s

500 000 - 19h 40mn

	SUR SUN SPARC 10, 32 BITS	SUR DEC ALPHA OSF, 64 BITS
DECIMALES	TEMPS	TEMPS
1 500	1.8 s	0.5s
10 000	1mn 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.