Calcul de Pi - forth et multiprécision

Le problème

Quand on utilise bc, si l'on veut que le résultat d'un calcul soit effectué avec 1000 décimales, il suffit d'entrer la bonne valeur avec le paramètre scale.
Ensuite, le programme s'occupe de tout.

Mais comment faire avec un langage informatique ?
Parce que, là, le programmeur a la responsabilité de tout.

La réponse d'Astrobe

Rappel : le site d'Astrobe a pour adresse : http://membres.lycos.fr/astrobe/

Je reproduis ici, une réponse faite par Astrobe, à une question posée un peu en l'air.
Elle a le mérite d'être à la fois simple (donc compréhensible) et pédagogique.
Ce serait dommage que d'autres ne puissent en profiter.

Ma question :

> Parmi les formules que l'on trouve sur internet, il en est
> une qui ne fait appel qu'à des calculs simples :
> pi = 4 (1 - 1/3 + 1/5 - 1/7 +  1/9 ...)
> Il paraît qu'elle converge lentement...
> Question : comment réaliser, le plus simplement possible le calcul
> 1/3  ou 1/5  ou 1/7 ... avec 4IM en étant le plus rapide
> possible (et le plus exact).

La réponse :

Le souci est que 4IM ne connait pas les nombres décimaux, et encore moins
les flottants. Une qualité de l'exemple en GForth que tu as trouvé est qu'il
travaille en nombres entiers (à propos, il faudrait se poser la question de
savoir s'il a besoin d'entiers 32 bits ou non ). Si on fait 1 3 /, on
obtient 0.

Qu'à cela ne tienne, on va calculer 10*pi :
10pi= 4 * (10 - 10/3+10/5-10/7+10/9 ... )
    = 4 * ( 10 - 3 + 2 - 1 + 1 ) en faisant des divisions entières
    = 36 c'est a dire 3.6

Pas terrible. Il faudrait passer a 100pi pour avoir la première décimale
correcte. Vu que nos nombres sont limités à 65535, on ne peut espérer guère
mieux qu'obtenir 31416... ( 4 décimales, snif! ).
Mais une idée serait de calculer (pi-3)*10000, ce qui nous permettrait
d'obtenir une décimale de plus. En appliquant cette idée plusieurs fois, on
pourrait calculer les décimales de proche en proche.

Exploration d'une solution avec Forth

Forth peut être utilisé en mode interprété ou en mode compilé. Le mode interprété permet de manipuler, sans avoir le souci de trop structurer la solution.

Version choisie et raisons

Ce sera gforth (qui calcule sous 32 bits), pour les raisons suivantes :

Choix du calcul à effectuer

On calculera 1/17 avec 20 décimales.
La valeur de comparaison est donnée par bc :

scale=20
1/17
.05882352941176470588

Déroulement du calcul avec gforth

Commençons pas "voler", à Augustin Vidovic, la routine d'affichage d'un nombre de 4 chiffres
: 4. 0 <# # # # # #> ;
Les 4 dièses centraux correspondent chacun à l'affichage d'un chiffre.

Entrer 1 et le multiplier par 10000 (voir texte Astrobe). Le mot .s sera utilisé pour afficher le contenu de la pile sans y provoquer de modification.
1 10000 *
Affichage de la pile
.s
<1> 10000 ok
La pile contient un seul élément (c'est le sens du <1> et cet élément est 10000). Forth nous dit qu'il est prêt pour la suite (c'est le sens du ok).

Effectuons la division entière du contenu de la pile par 17 en empilant à la fois le quotient entier en la valeur entière du reste. Le mot gforth utilisé est /mod On utilisera à nouveau .s pour vérifier ce qui se passe.
17 /mod .s <2> 4 588 ok
Le quotient 588 est sur le sommet de la pile. Le reste 4 est en dessous.

Affichons le quotient entier, en l'enlevant du sommet de la pile (rappelons que le mot 4. permet d'afficher 4 chiffres, en "mangeant" le sommet de la pile) :
4. 0588 ok
La première série de 4 chiffres est 0588
Vérifions le contenu de la pile
.s <1> 4 ok

La suite : l'hypothèse de départ est-elle la bonne ?

Si oui, il suffit de multiplier le reste par 10000 et de répéter les manipulations.
Copie de la console gforth :
10000 * ok
17 /mod .s <2> 16 2352 ok
4. 2352 ok
.s <1> 16 ok

La nouvelle série de 4 chiffres est 2352

Poursuivons sur notre console gforth :
10000 * ok
17 /mod .s <2> 13 9411 ok
4. 9411 ok

10000 * ok
17 /mod .s <2> 1 7647 ok
4. 7647 ok

10000 * ok
17 /mod .s <2> 4 588 ok
4. 0588 ok

Résultat

En réunissant les chiffres affichés 4 par 4 avec le forth, on obtient la valeur
05882352941176470588
qui correspond aux décimales obtenues avec bc.
Remarquons qu'il nous manque la partie le "début" du résultat, qui, dans la réalité est : 0.05882352941176470588

Utilisation de boucles

Afin de s'éviter les saisies répétitives au clavier, il est possible d'utiliser une boucle contrôlée par un compteur.
En gforth, cela peut s'écrire :
maxi mini do (code...) loop

La contrainte est ici qu'il faut placer la boucle dans un mot compilé.

Décomposons ce dont nous avons besoin.

Le mot d'affichage

Pas de surprise ici. C'est celui qui a déjà été utilisé :
: 4. 0 <# # # # # #> type ;

Le mot pour la partie calcul et affichage

Il découle des manipulations précédentes :
: div17 10000 * 17 /mod 4. ;

Le mot qui effectue la boucle et appelle le calcul

La boucle sera répétée 5 fois, ce qui se traduit par 5 0 do
: inverse17 5 0 do div17 loop ;

Application

Vérifions ce que cela donne. Pour cela, il faut placer 1 sur la pile (le numérateur de la fraction 1/17) :
1 inverse17 05882352941176470588 ok

Remarques et critiques

Il importe de tempérer l'enthousiame sur ce qui précède (remarques inspirées d'Astrobe).

Le risque de dépassement de valeur

Dans le cas du calcul des décimales de Pi, la multiplication par 10000 du reste de la division ne provoque pas de dépassement de la valeur maximale sur 32 bits. Mais dans le cas de nombre trop grand au départ, il serait prudent de vérifier les risques de débordements.

Le ralentissement dû aux boucles répétitives

Chuck Moore défend l'idée qu'à chaque fois que l'on évite une boucle on gagne sur le temps d'exécution. Dans le cas d'un calcul répété 100 fois, il préférerait
: inverse17 20 0 do div17 div17 div17 div17 div17 loop ;
à
: inverse17 100 0 do div17 loop ;

Dans ce cas, le programmeur doit choisir entre l'efficacité et le respect de l'orthodoxie...

Adaptation à 4IM

Le Forth 4IM est présenté sur d'autres pages de ce site. Rappelons qu'il calcule sur 16 bits et non sur 32 bits comme gforth.
En conséquence, nous allons multiplier par 100 (et non par 10000) et progresser pa deux chiffres à la fois (au lieu de 4).

Le mot d'affichage

Il devient :
: 2? <# # # DROP #> "? .

Le mot pour la partie calcul et affichage

Il devient :
: DIV17 100 * 17 /MOD 2? .

Le mot qui effectue la boucle et appelle le calcul

Nous allons tricher et sauter la boucle :
: INVERSE17 DIV17 DIV17 DIV17 DIV17 DIV17 .

Que nous appliquerons deux fois de suite (en utilisant le sommet de pile, la seconde fois... :
1 INVERSE17
0588235294 Ready
INVERSE17
1176470588 Ready

Ce qui donne successivement
0588235294 puis 1176470588 (voir plus haut).
Remarque :
1 INVERSE17 INVERSE17
donne directement :
05882352941176470588 Ready

Utilisation d'une boucle

4IM adopte une syntaxe particulière pour la réalisation de boucles.
Le mieux est ici de lire la domcumentation mise en ligne par Astrobe, ainsi que les notes d'expérimentation de mon propre site (elles progressent lentement).

Voici un mot qui fonctionne bien :
: MOT 5 :: DUP IF >R DIV17 DIV17 R> 1- IDEM ; THEN DROP .
Ready
1 MOT
05882352941176470588Ready

Remarque :
La valeur de contrôle de la boucle (initialement 5, puis 4 ... jusqu'à 0) étant sur la pile, on utilise la pile de retour (par les mots >R et R>) de façon à la sauvegarder momentanément.
Cette façon de faire est très classique en forth et est expliquée sur les documents permettant l'apprentissage du langage (mais pas sur mon site :-))

Conclusion

Le calcul en multiprécision est possible en Forth :
- Avec un Forth 32 bits tel que gforth (sous Linux ou Windows),
- Avec un forth 16 bits tel que 4IM (soit en version autonome sans aucun système, soit sous Dos).
4IM peut fonctionner sous Linux, via l'émulateur dosemu.


Réagir :lerautal A suivre....