Suivez-nous sur X
|
|
|
0,
A,
B,
C,
D,
E,
F,
G,
H,
I,
J,
K,
L,
M,
N,
O,
P,
Q,
R,
S,
T,
U,
V,
W,
X,
Y,
Z,
ALL
|
|
0,
A,
B,
C,
D,
E,
F,
G,
H,
I,
J,
K,
L,
M,
N,
O,
P,
Q,
R,
S,
T,
U,
V,
W,
X,
Y,
Z
|
|
0,
A,
B,
C,
D,
E,
F,
G,
H,
I,
J,
K,
L,
M,
N,
O,
P,
Q,
R,
S,
T,
U,
V,
W,
X,
Y,
Z
|
|
A propos d'Obligement
|
|
David Brunet
|
|
|
|
Programmation : Assembleur - Calculer précisément Pi
(Article écrit par Franck Chatton et Alain Mignon et extrait d'Amiga News Tech - juillet 1991)
|
|
Notre concours lecteur nous entraîne ce mois-ci au pays des mathématiques. Franck Chatton et Alain Mignon
proposent en effet un programme qui sort des sentiers battus : il calcule le nombre Pi
avec une précision pouvant atteindre 45 000 décimales !
La première tentative de calcul de Pi vraiment raisonnable, après les pifométrages à base de fraction de
l'antiquité (Lambert n'ayant prouvé l'irrationalité de Pi qu'en 1761), semble revenir à
Leibniz qui utilise le développement en série de l'arctangente : sachant que Tan(Pi/4)=1,
Pi/4=1-1/3+1/5-1/7+1/9. Mais on a trouvé largement mieux depuis. Signalons que Schanks
a calculé en 1874 les 707 décimales (ne rigolez pas, il n'avait pas d'ordinateur ! N'empêche qu'elles
pourraient quand même être justes...) qui sont gravées dans une salle du palais de la découverte.
Un peu de mathématiques
Nous n'allons surtout pas utiliser la formule de Leibniz qui est, en matière de convergence,
encore plus un veau qu'une Harley, ni même celle d'Euler (Pi/4=Atan(1/2)+Atan(1/3))
qui se traîne encore passablement, mais celle de Machin : Pi/4=4Atan(1/5)-Atan(1/239).
Pour calculer l'arctangente, nous allons utiliser son développement en série entière,
sachant que Atan(x)=SIGMA(-1)^i/(2i+1)*x^(2i+1). Enfin, pour calculer les arctangentes,
nous n'allons pas calculer bêtement la somme telle que je viens de l'écrire, mais nous
allons utiliser l'algorithme de Horner, par exemple pour calculer aX^2+bX+c, nous ne
calculons pas a*X*X+b*X+c, mais ((a)*X+b)*X+c, ce qui est la manière la plus économique en
nombre d'opérations donc en temps, d'évaluer un polynôme.
Petit historique
Initialement, nous disposions d'un programme de calcul de Pi (merci à Laurent Charbonnel de
nous avoir fourni le source), qui bien qu'un peu lent (ah, les subtilités de l'euphémisme :
plus de 10 heures de calcul pour... 5000 décimales !), ridiculisait déjà
notre polytechnicien de prof de maths prétentieux, avec son super logiciel "c'est môâââ qui
l'ai fait" sur (âmes sensibles s'abstenir) XT en Turbo Pascal et qui pendant ce temps calculait
poussivement ses 2000 décimales, et encore... Ce programme laissait une large place à
l'optimisation, et la version actuelle n'a plus beaucoup (et même plus du tout) de points
communs avec son ancêtre éloigné, hormis bien sûr la fonction.
Un peu de pratique
L'utilisation de la méthode de Horner nous oblige à calculer le rang jusqu'auquel nous devons
pousser le développement de la série pour obtenir la précision requise avant de commencer le
calcul. Le rang maximal a été calculé en tenant compte des erreurs diverses et vous garantit une
précision de 1E-p quand vous demandez p décimales.
"Oh l'autre eh, comment y va les stocker ses décimales, ça passe même pas en double précision !"
pensez-vous. Remarque éminemment pertinente au demeurant. Nous allons regrouper les décimales
par groupes de quatre dans des mots, eux-mêmes situées dans des zones de mémoire que nous aurons
allouées.
Attention, ce n'est pas du DCB. Mais alors pourquoi quatre et pas cinq ? Parce que 65 535<10^5.
Les heureux possesseurs de 68020 et plus pourront encore améliorer l'algorithme
car les mots longs permettent de stocker neuf décimales, mais le 68000 ne peut malheureusement
multiplier et diviser que des mots, d'où la limitation à 45 000 décimales de ce programme
(il est possible de faire plus, mais il faut déjà environ neuf heures avec ce programme et il
serait notablement ralenti par cette nouvelle modification).
L'arctangente étant impaire et le signe des coefficient alterné, nous allons améliorer la
méthode de Horner en n'écrivant pas aX-bX^3+cX^5+...=X*(a+X*(0+X*(-b+X*(0+X*(c+...)))))
mais aX-bX^3+cX^5+...=X*(a-X2*(b-X2(c-...))))), avec X2=X^2 (étonnant, non ?) calculé une
fois pour toutes au début. Ruse subtile comme nous allons calculer deux arctagentes,
pourquoi calculer deux fois les mêmes coefficients. Une partie des coefficients 1/(2i+1)
étant commune a deux développements, leur calcul sera mené simultanément à partir d'u certain rang.
Attention, pour les possesseurs d'A500 sans horloge, le chronomètre ne
fonctionnera pas, mais il ne provoque cependant pas de plantage. Enfin, le but initial de
notre programme étant la performance pure, nous avons encore gratté un peu de temps de calcul
au système à coups de Forbid et de Disable, le tout assaisonné d'une suppression bestiale
de DMA d'écran qui piquait quelques cycles au 68000 en haute résolution, avec pour effet
secondaire de fournir un économiseur d'écran pour les distraits et un témoin d'exécution du programme.
|
