Pi et Plouffe

[Gilles59]

En 1995 Simon Plouffe découvre une nouvelle formule pour calculer la valeur de Pi, c'est la formule de Bailey-Borwein-Plouffe (BBP).
Cette formule est étonnante pour au moins deux raisons:

1/ Sa simplicité ....



A tel point qu'il est surprenant qu'elle n'ait pas été découverte avant !

2/ Elle permet (ce qui semblait impossible alors, y compris pour les plus grands 'spécialistes' de PI) de calculer le n-ième chiffre de PI sans avoir besoin de connaitre les chiffres précédents... mais attention, cela fonctionne efficacement en base 16 (ou 8 ou 4 ou 2) seulement

http://fr.wikipedia.org/wiki/Simon_Plouffe

J'avoue ne pas avoir bien compris (en tout cas réussi à écrire un programme HP50G) comment calculer le n-ième chiffre hexadécimal (ou binaire) de PI mais ci-après quelques variations sur ce thème sur la 50G

Le code HP50 pour calculer la formule de Plouf en "précision infinie" (mode EXACT)

Code : Tout sélectionner

« 0 DUP ROT FOR k 
     '(1/(16^k))*((((4/((8*k)+1))-(2/((8*k)+4)))-(1/((8*k)+5)))-(1/((8*k)+6)))' EVAL + 
  NEXT
»
'~PI' STO

Par exemple
5 ~PI donne :

'47/15+53/6552+829/5026560+79/15590400+857/4561108992+1901/244729774080'

EVAL

Code : Tout sélectionner

'40413742330349316707/12864093722915635200'

->NUM

Code : Tout sélectionner

3.14159265323
(les 2 derniers chiffres son faux dans ce cas)

On peut aller plus loin en précision et on dépasse alors la précision des 12 chiffes affichables en mode réel de la 50

Voilà un petit prog pour aller plus loin, vers une précision 'infinie' ( ou plutot limitée par la mémoire et le temps de calcul). ici par exemple on calcule a/b avec 240 décimales :

Code : Tout sélectionner

« ""  -> Pi
 «
  FXND                       @ Transforme a/b en  a b
  SWAP OVER IDIV2 
  0 240 FOR n
   SWAP STR 'Pi'SWAP STO+
   10 * OVER IDIV2 
  NEXT
  3 DROPN Pi
 »
»
'Calc' STO

100 ~PI CALC

donne alors :

Code : Tout sélectionner

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628
03482534211706798214808651328230664709384460955058223172535940812848111745028410270193
8521105559644622948954930381964428810975665933446128475648233786783165

Arrivé là deux questions :

1/Comment savoir _à priori_ le nombre de chiffres corrects de PI calculé de cette façon ?

2/ Mais ce qui m'intéresse vraiment serait l'algo qui donne le n-ieme chiffre hexa de PI.... Là,je sèche,Wikipédia me parait peu clair

[zpalm]

Une très belle et remarquable formule. J'étais tombé il y a quelques temps sur l'histoire de sa découverte par Simon Plouffe lui-même: The story behind a formula for Pi .
Une histoire un peu triste sur les rapports humains ....

Edit: pour répondre à ta question sur l'algorithme :

• The BBP Algorithm for Pi
• BBP Code Directory

(google est ton ami ne l'oublie pas )

[rogeroge]

Bonsoir...
Bel exposé sur Pi et aussi Plouffe, personnalité scientifique que je découvre.
Je ne suis pas un matheux de haut niveau mais la formule de Leibniz n'est-elle pas plus facile à mettre en oeuvre sur une calculatrice pour répondre au seul critère de la simplicité ?
Pi/4 = 1 - 1/3 + 1/5 - 1/7 + ...

[Gilles59]

'Pi/4 = 1 - 1/3 + 1/5 - 1/7 + ...'
J'ai testé et la formule de leibnitz converge beaucoup moins vite

Avec une suite de 100 fractions j'obtiens 'seulement' 3.15149340106 contre 240 décimales exactes pour l'autre formule

Mais il est vrai que le gros avantage de la formule BBP est de pouvoir calculer assez facilement (bien que je n'ai pas réussi) en hexadécimal un chiffre de PI sans calculer tous les précédents.

[C.Ret]

Très bien présenter ce sujet sur PI et la fomule de Plouffe !
Il faut que je prenne le temps de lire les référence indiquée, cela semble très instructif.


Sinon, je suis d'accord les formules comme celles de Libtniz ne convergent pas vite, et tout cas pas aussi vite que les formules plus récentes.
Par contre, l'une d'entre-elles permet de composer un algorithme très simple car il se compose de la succession de deux opérations simples (une division et une multiplication). Ce qui permet de programmer cela afin d'utiliser un grand nombre de rgistre pour calculer (en prenant le temps) un nombre conséquent de décimales.



Cela me rappèle le Misez P'tit, Optimisez n°14 - une approximation de PI et le seul code qui permet de calculer (en théorie) jusqu'à 1044 décimales de PI sur un SHARP PC-1211:

Par chance, elle (la méthode) est bien expliquée sur ce site : http://serge.mehl.free.fr/anx/pi_100deci.html

J'ai donc adapté l'algorithme donné pour mon pauvre SHARP PC-1211.
Cela donne un programme d'environ 217 STEPs qui laisse donc disponible 175 MEMORIES pour stocker les dicimales de PI, soit un maximum théorique d'environ 1044 décimales. Théorique, car personne n'aura la patience (et les piles) pour arriver à un tel résultat. Il faut quelque(s) fraction(s) de seconde sur un PC pour avoir aussi peu de décimale !

Il doit être possible d'optimiser ce code ou de gagner quelque mémoire, mais est-ce bien utile ?

L'avantage de cette version, qui utilise le même algorithme que le programme de SVI (*) , est que l'on ne pert pas de temps à calculer les décimales au-dela de la limite demandée. Le nombre de registre utilisé A(6 à E) dépend du nombre demandé. Les calculs sont donc rapides pour un nombre faible de décimales.

Code : Tout sélectionner

  REM ---------------  saisie nombre decimales

1:CLEAR :INPUT D

  REM --------------- initialisation - D = nmax, E nombre de block, T(0 à ...) = A(7 à E)=2n_max/(2n_max+1)  

2:E=7+INT .2D,D=INT (3+D/LOG 2),G=2D,F=G+1:GOSUB 8

  REM --------------- boucle principale  To=To+2, divise T() par 2n-1, multiplie T() par n-1

3:PAUSE D:G=G+2:IF D<2 BEEP 1:FOR A=7 TO E:PRINT A(A):NEXT A:END
4:F=2D-1:GOSUB 8:F=D-1:GOSUB 9:D=D-1:GOTO 3

  REM - - - - - - - - ss-prg division  A(7 à E) par F

8:C=0:FOR A=7 TO E:A(A)=A(A)+€6C,B=INT (A(A)/F),C=A(A)-FB,A(A)=B:NEXT A:RETURN

  REM - - - - - - - - ss-prg multiplication A(7 à E) par F

9:C=0:FOR A=E TO 7 STEP -1:B=FA(A)+C,C=INT €-6B,A(A)=B-€6C:NEXT A:RETURN

( € correspond en fait à la touche Exp qui permet de saisir les exposants de 10).

Variables:
A: compteur i
B: quotient/produit/valeur intermédiaire
C: retenue
D: nbr décimales de PI souhaitée/ nmax/ n
E: nombre/limite registre stokage des décimales: T0=A(7)=G .
F: facteur ou dividente pour respectivement division et multiplication.
G à A(177) mémorisation de PI avec jusqu'à presque 1000 décimales. En fait 6 décimales par registre mémoire (valeur absolue pour simplifier expression pour des calculs plus rapides).


(*) M. Labat (de PARIS) et publié dans Science&Vie n°759 (pdf)

[Marge]

'Pi/4 = 1 - 1/3 + 1/5 - 1/7 + ...'
J'ai testé et la formule de leibnitz converge beaucoup moins vite

Ah oui, je confirme, c'est la plus facile à retenir mais la plus lente. Je crois qu'elle n'est pas de Leibniz, d'ailleurs, mais je n'ai plus son auteur en tête.

[yvesbolo]

Aaah Plouffe... Que de bons souvenirs !

J'ai fait un projet de processeur pour calculer Pi avec cette formule (sur FPGA). C'est peut-être ce que j'ai préféré durant toutes mes études

[rogeroge]

'Pi/4 = 1 - 1/3 + 1/5 - 1/7 + ...'
J'ai testé et la formule de leibnitz converge beaucoup moins vite

Ah oui, je confirme, c'est la plus facile à retenir mais la plus lente. Je crois qu'elle n'est pas de Leibniz, d'ailleurs, mais je n'ai plus son auteur en tête.

J'ai regardé dans mes vieux cours qui datent de 1967 à 1969 :
Voici les formules simples des valeurs approchées de Pi :

Leibniz : Pi / 4 = 1 - 1/3 + 1/5 - 1/7 + ...+ ((-1)^n) / (2n + 1) + ...

Wallis : Pi / 2 = (2^2) / (1x3) x (4^2) / (3x5) x (6^2) / (5x7) ... ((2n)^2) / ((2n-1) x (2n+1)) ...

Je pense aussi que celle de Wallis est aussi lente.

[Marge]

"Le mathématicien écossais James Gregory (1638 - 1675), professeur à l'Université de Saint Andrews puis à celle d'Edimbourg, fut notamment l'inventeur d'un télescope à miroir secondaire concave. Il tenta en vain de démontrer que le problème de la quadrature du cercle était impossible à résoudre ; il prétendit même avoir réussi et publia une démonstration qui ne persuada ni Huygens, ni Leibniz. [...]

James Gregory a aussi découvert la formule suivante :

arctan (x) = x - x³/3 + x⁵/5 - x⁷/7 + ... [...]

par une méthode que l'on interprète aujourd'hui comme le calcul de arctan(x) en tant que primitive de 1/(1+x²) = 1 - x² + x⁴ - x⁶ + x⁸ -...
[...]

En prenant x=1, on obtient l'extraordinaire somme infinie :
Pi = 4 (1 - 1/3 + 1/5 - 1/7 +...) [...]

Malheureusement, Gregory ne l'a jamais écrite explicitement. Peut-être est-ce parce qu'il avait compris qu'elle n'était guère utile pour calculer pi. [...] La convergence est exécrable ! On parle de convergence logarithmique : le nombre d'étapes pour gagner un chiffre est de plus en plus grand. [...]

La formule de arctan(x) de Gregory sera ultérieurement utilisée pour le calcul de Pi, mais avec des valeurs inférieures à 1, car plus x est proche de 0, plus elle converge.
En fait, la formule pour Pi avait déjà été proposée vers 1410 par le mathématicien indien Madhava, mais elle était restée inconnue en occident."

Jean-Paul Delahaye, Le fascinant Nombre Pi, Belin, 2000, p. 69.

Plus loin (à la page suivante), l'auteur affirme qu' "en 1674, [Leibniz] propose la formule de Pi que Gregory était sur le point d'écrire [ce qui aura précipité la fin de ce dernier ?]. Cette formule, que certains nomment formule de Gregory, d'autres formule de Leibniz, (et que je choisis de nommer formule de Madhava-Gregory-Leibniz) ne fut pas la seule rencontre de l'inventeur du calcul différentiel avec Pi."

L'histoire des sciences est emplie de ces injustices, de ces orgueils et amours-propres blessés !
(à voir la convergence de celle de Wallis, je crois que tu as raison, Rogeroge, elle paraît aussi logarithmique).

Excellent livre que celui-là, et on y parle aussi de Simon Plouffe (et ce dernier parle aussi de ce livre dans le lien de zpalm : le monde est presque rond !).

[Gilles59]

L'histoire de PI est vraiment passionnante.

Ci après un lien qui contient plusieurs programmes pour le calcul de PI sur une HP50G
http://sense.net/~egan/hpgcc/

Un des programme trouve 5000 décimales en 46 sec (le plus rapide les trouve en 38s)


Bien sur, c'est en C. Les programmes sont bien commentés et comparés

Si c'est programmes sont très rapides, à noter que le lien S&V donné par C.RET donne un prog 41C qui en calcule 3600 décimales soit pas beaucoup moins mais en plusieurs mois :O ...

Lu plus en détail l'article....On peut monter à "200,001 correct digits in 17.44 hours." ...

Il y aussi une comparaison entre ARM/C et Saturn. Le program en C est 15 fois plus rapide,mais 20 fois plus volumineux ( 20 ko au lieu de 1 Ko)

Le site est aussi intéressant car il explique vraiment en détail comment programmer en C sur la 50G.

[gege]

Pour savoir quelle est l'erreur en ignorant les termes à partir du rang n, on peut commencer par calculer le terme de rang n justement.
L'erreur sera forcément supérieure (car les termes sont tous positifs).
Vu la décroissance très rapide des termes, la vraie erreur ne sera pas beaucoup plus grande que le terme de rang n.
G.E.

[Gilles59]

J'ai largement simplifié le programme RPL du calcul de PI basé sur BBP (mais il y a semble-t-il plus rapide pour calculer les décimales de PI):

Code : Tout sélectionner

 
«
 2000 ALOG
 0 DUP 1700 FOR k
  k 1. DISP
  k DUPDUP 151 * SWAP SQ DUP 120 * ROT + 47 + 5 PICK *
  UNROT
  DUP SQ 512 * UNROT DUP2 * 1024 * UNROT 712 * SWAP 194 * 15 + + + + 16 k ^ *
  IQUOT +
 NEXT
»
'PI' STO

207 octets, chaque décimale de Pi prends un nibble (demi-octet)... Bref çà n'utilise quasi pas de RAM (quelques kilo)

Ce programme calcule les 2000 première décimales PI en quelques minutes sur l'émulateur (258 sec. sur mon portable assezancien) ,pas testé sur la vraie calc).

J'ai essayé de monter à 10.000 décimales (ce programme est tres économe en mémoire et avec 200 ko on pourrait en calculer en théorie des centaines de milliers), mais je me suis aperçu que les enitiers 'exacts' ne sont pas infinis : 10^10000 génére une erreur :entier trop grand ! 5000 est faisable

Pour déterminer le nombre de boucles nécessaire j'utilise :
2+log(8n)+n*log(16)
trouvé ici : http://www.pi314.net/fr/plouffe.php

Code : Tout sélectionner

31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066
47093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783
16527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259
036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575
27248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694
05132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129
02196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334
46850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818
57780532171226806613001927876611195909216420198938095257201065485863278865936153381827968230301952035301852968995773622
59941389124972177528347913151557485724245415069595082953311686172785588907509838175463746493931925506040092770167113900
98488240128583616035637076601047101819429555961989467678374494482553797747268471040475346462080466842590694912933136770
28989152104752162056966024058038150193511253382430035587640247496473263914199272604269922796782354781636009341721641219
92458631503028618297455570674983850549458858692699569092721079750930295532116534498720275596023648066549911988183479775
35663698074265425278625518184175746728909777727938000816470600161452491921732172147723501414419735685481613611573525521
33475741849468438523323907394143334547762416862518983569485562099219222184272550254256887671790494601653466804988627232
79178608578438382796797668145410095388378636095068006422512520511739298489608412848862694560424196528502221066118630674
427862203919494504712371378696095636437191728746776465757396241389086583264599581339047802758185

[/size]

PS: il manque une itération pour que les 4 dernières décimales soient bonnes

[gege]

Ben voilà, ce nombre de termes correspond au premier terme :
log(1/16^n*4/8n)=1.66+log(8n)+n*log(16)
G.E.

[C.Ret]

Bonjour,

J'ai un peu de mal avec cette formule:

Code : Tout sélectionner

n   :   1     2     3     4      10      82     100       165       247
log :  -1.51 -3.01 -4.39 -5.72  -13.3  -100.9  -122.7    -201.2    -300.1

Faut-il comprendre que pour obtenir 200 décimales correctes de PI, il faut utiliser n au moins égal à 247 ?

[Gilles59]

Bonjour,



Faut-il comprendre que pour obtenir 200 décimales correctes de PI, il faut utiliser n au moins égal à 247 ?

Oui c'est bien çà !