Voir la source de Formule BBP

{{à recycler|date=octobre 2020}}
{{Sources|date=septembre 2013}}
{{Voir homonymes|BBP}}

La '''formule BBP''' (ou '''formule de Bailey-Borwein-Plouffe''') permet de calculer le ''n''-ième chiffre après la virgule du [[Pi|nombre {{math|π}}]] en [[base 2]] (ou [[base 16|16]]) sans avoir à en calculer les précédents, et en utilisant très peu de mémoire et de temps. Elle a été obtenue le {{Date|19|septembre|1995}} par [[Simon Plouffe]] en collaboration avec [[David H. Bailey]] et [[Peter Borwein]]<ref>{{article|langue=en
 | nom1 = Bailey | prénom1 = David H. | nom2 = Borwein | prénom2 = Peter B. | nom3 = Plouffe | prénom3 = Simon
 | année = 1997
 | titre = On the Rapid Computation of Various Polylogarithmic Constants
 | journal = [[Liste des journaux scientifiques en mathématiques #M|Math. Comp.]]
 | volume = 66 | numéro = 218 | pages = 903-913 <!-- | url = http://crd.lbl.gov/~dhbailey/dhbpapers/digits.pdf   PDF is free at doi link, and is better formatted. -->
 | doi = 10.1090/S0025-5718-97-00856-9
 | math reviews= 1415794
 }}.</ref>.

== La formule ==
Dans sa forme originelle, la formule BBP est donnée par

:<math>\pi=\sum_{k=0}^\infty\frac1{16^k}\left(\frac4{8k+1}-\frac2{8k+4}-\frac1{8k+5}-\frac1{8k+6}\right)</math>.

{{démonstration|
Notons <math>S_n=\sum_{k=0}^\infty\frac1{16^k(8k+n)}</math> et démontrons la formule de Plouffe généralisée :
:<math>(0)\quad\forall r\in\Complex\qquad\pi\ =\ (4+8r)S_1-8rS_2 -4rS_3 -(2+8r)S_4-(1+2r)S_5-(1+2r)S_6+rS_7</math>
(le cas {{math|''r'' {{=}} 0}} est sa formule originale ; le cas {{math|''r'' {{=}} –1/4}} est, sous une forme plus détaillée, celle d'Adamchick et Wagon).

Posons <math>\alpha=1-i=\sqrt 2\mathrm e^{-\mathrm i\pi/4}</math> et calculons de deux façons l'intégrale suivante :
<center><math>I=\int_0^1\frac{\mathrm dy}{\alpha-y}.</math></center>
Elle est d'une part reliée aux {{mvar|S{{ind|n}}}} par
:<math>\begin{align}I&=\frac 1\alpha\int_0^1\frac{\mathrm dy}{1-y/\alpha}=
\frac1\alpha\int_0^1\sum_{m=0}^\infty\frac{y^m}{\alpha^m}\mathrm dy=\sum_{m=0}^\infty\frac{\mathrm e^{\mathrm i(m+1)\pi/4}}{(m+1)\sqrt 2^{m+1}}\\&=
\frac{1+\mathrm i}2S_1+\frac {\mathrm i}2S_2+\frac{-1+\mathrm i}4S_3-\frac14S_4-\frac{1+\mathrm i}8S_5-\frac {\mathrm i}8S_6+\frac{1-\mathrm i}{16}S_7+\frac1{16}S_8,\end{align}</math>
et d'autre part calculable par des méthodes élémentaires (en calculant séparément sa partie réelle et sa partie imaginaire), ou de façon plus synthétique via le [[logarithme complexe]] :
:<math>\begin{align}I&=-\left[\ln(\alpha-y)\right]_0^1=\ln\left(\frac\alpha{\alpha-1}\right)=\ln(1+\mathrm i)=\ln(\sqrt 2\mathrm e^{\mathrm i\pi/4})=\frac{\ln 2}2+\mathrm i\frac\pi4.\end{align}</math>

L'égalité entre ces deux expressions de {{mvar|I}} équivaut à :
:<math>(1)\qquad\pi=4\mathrm{Im}(I)=2S_1+2S_2+S_3-\frac 12S_5-\frac 12S_6-\frac 14S_7,</math>
:<math>(2)\qquad\ln 2=2\mathrm{Re}(I)=S_1-\frac 12S_3-\frac 12S_4-\frac 14S_5+\frac 18S_7+\frac 18S_8.</math>

Mais {{math|ln 2}} s'exprime par ailleurs plus directement en fonction des {{mvar|S{{ind|n}}}} :
:<math>\begin{align}(3)\quad\ln 2&=[-\ln(2-y^2)]_0^1=
\int_0^1\frac{y}{1-y^2/2}\mathrm dy=\sum_{k\ge 0}\frac 1{(2k+2)2^k}\\&=S_2+\frac 1 2 S_4+\frac 1 4 S_6+\frac 1 8 S_8.\end{align}</math>
(2) et (3) donnent donc, par soustraction des membres de droite, une relation entre les {{mvar|S{{ind|n}}}} :
:<math>(4)\qquad 0=2S_1-2S_2-S_3-2S_4-\frac 1 2 S_5-\frac 1 2 S_6+\frac 1 4 S_7.</math>
En multipliant le membre de droite de (4) par {{math|1+4''r''}} et en ajoutant ce produit au membre de droite de (1), on obtient l'égalité (0) annoncée.
}}

== Exploitation de la formule pour calculer les chiffres après la virgule de {{math|π}} ==
{{Article détaillé|Approximation de π{{!}}Approximation de {{math|π}}}}

Le but est de calculer le ''N''-ième chiffre après la virgule de {{math|π}} en [[base 16]].

Déjà, on remarque que le (''N'' + 1)-ième chiffre après la virgule de {{math|π}} en base 16 est le même que le {{1er|chiffre}} après la virgule de 16<sup>''N''</sup>{{math|π}}. En effet, comme en [[base 10]], multiplier un nombre en base 16 par 16 permet de décaler la virgule d'un rang vers la droite. En multipliant un nombre par 16<sup>''N''</sup>, la virgule est donc décalée de ''N'' rangs vers la droite. Ainsi, il suffit de calculer le premier chiffre de 16<sup>''N''</sup>{{math|π}}, égal par la formule BBP à :

<center><math>16^{N}\pi = \sum_{k=0}^\infty 16^{N-k}\left(\frac4{8k+1}-\frac2{8k+4}-\frac1{8k+5}-\frac1{8k+6}\right)</math>.</center>

Mais calculer les premiers chiffres derrière la virgule de ce nombre n'est pas si simple, pour deux raisons :
*d'abord, ce nombre étant très grand, cela demande d'effectuer des calculs sur des nombres très grands ;
*ensuite, parce que cette somme est infinie.

Posons <math>S_N(a)=\sum_{k=0}^\infty \frac{16^{N-k}}{8k+a}</math>. Le calcul des premiers chiffres de {{math|''S<sub>N</sub>''(''a'')}} permettra d'obtenir ceux de 16<sup>''N''</sup>{{math|π}}, par la relation :

<center><math>16^N\pi = 4\ S_N(1) - 2 S_N(4) - S_N(5) - S_N(6)</math>.</center>

Découpons la somme {{math|''S<sub>N</sub>''(''a'')}} en deux :

<center><math>S_N(a)=\sum_{k=0}^\infty \frac{16^{N-k}}{8k+a}=\sum_{k=0}^{N-1} \frac{16^{N-k}}{8k+a}+\sum_{k=N}^ \infty \frac{16^{N-k}}{8k+a} = A_N(a) + B_N(a)</math></center>

et calculons {{math|''A<sub>N</sub>''(''a'')}} et {{math|''B<sub>N</sub>''(''a'')}} indépendamment.

=== Calcul de {{math|''B<sub>N</sub>''(''a'')}} ===
<math>B_N(a)=\sum_{k=N}^ \infty \frac{16^{N-k}}{8k+a}</math>

Bien que ce soit une somme infinie, ce terme est très simple à calculer, car on remarque que ses termes deviennent vite très petits et on ne cherche que les premiers chiffres.

* En effet, le premier terme de la somme est : <math>b_N=\frac1{8N+a}</math>. Comme on cherche le ''N''-ième chiffre derrière la virgule de {{math|π}} (''N'' = {{formatnum:1000000000}} par exemple), le premier terme ''b<sub>N</sub>'' est très inférieur à 1.
* De plus, chaque terme suivant a un zéro de plus derrière la virgule que le précédent, car pour ''k'' ≥ ''N'', ''b<sub>k</sub>'' > 16 ''b''<sub>''k''+1</sub> :
<center><math>\frac{b_k}{b_{k+1}}=\frac{16^{N-k}}{16^{N-(k+1)}}\frac{8(k+1)+a}{8k+a}=16\left(1+\frac8{8k+a}\right)\longrightarrow16^+</math>.</center>

Finalement, la somme {{math|''B<sub>N</sub>''(''a'')}} est de la forme (au pire) :

<math>B_N\ =\ \ 0,**********.\ .\ .\ .\ </math>

::<math>+\ 0,0*********.\ .\ .\ </math>

::<math>+\ 0,00********.\ .\ .\ </math>

::<math>+\ 0,000********.\ .\ .\ </math>

Donc pour obtenir {{math|''B<sub>N</sub>''(''a'')}} avec une précision de P chiffres derrière la virgule, il suffit de calculer les ''P'' premiers termes de la somme, plus les quelques suivants pour éviter les problèmes de retenues qui peuvent éventuellement apparaître.

Il suffit donc de calculer : <math>B_N'(a)=\sum_{k=N}^{N+P+10} \frac{16^{N-k}}{8k+a}</math>

Cette somme n'étant composée que d'un petit nombre de termes (de nombre constant), son temps de calcul est négligeable pour un [[ordinateur]].

=== Calcul de {{math|''A<sub>N</sub>''(''a'')}} ===
<math>A_N(a)=\sum_{k=0}^{N-1} \frac{16^{N-k}}{8k+a}</math>

Le problème pour calculer {{math|''A<sub>N</sub>''(''a'')}} est que les premiers termes sont extrêmement grands (''N'' chiffres en base 16 devant la virgule). Néanmoins, comme on ne cherche que les premiers chiffres derrière la virgule, peu importe la [[Partie entière et partie fractionnaire|partie entière]], aussi grande qu'elle soit. On peut donc s'en « débarrasser » en utilisant l'[[arithmétique modulaire]].

Toute la difficulté se réduit donc à trouver la partie fractionnelle de <math>\frac{16^{N-k}}{8k+a}</math>.

Pour cela, on effectue la [[division euclidienne]] de 16''<sup>N-k</sup>'' par 8''k''+''a'' :

:<math>\exists q\in\Z,\exists r < 8k+a,\ 16^{N-k}=q(8k+a)+r</math>

Donc <math>\frac{16^{N-k}}{8k+a}=q+\frac r{8k+a}</math>

:<math>\frac r{8k+a}</math> est inférieur à 1, donc c'est la partie fractionnelle de <math>\frac{16^{N-k}}{8k+a}</math>.

Et <math>\frac r{8k+a}=\frac{16^{N-k}\bmod(8k+a)}{8k+a}</math>

Il suffit donc de calculer : <math>A_N'(a)=\sum_{k=0}^{N-1} \frac{16^{N-k}\bmod(8k+a)}{8k+a}</math>.

En utilisant la méthode d'[[exponentiation rapide]], 16<sup>''N''-''k''</sup>mod (8''k''+''a'') se calcule rapidement (temps d'exécution en [[Comparaison asymptotique|O]](log{{ind|2}}(''N''-''k'')).

=== Conclusion ===
Finalement, pour obtenir les premiers chiffres de {{math|π}} en base 16 (ou 2), il faut calculer les premiers chiffres de :

<math>\pi_N=4\ S_N'(1)-2\ S_N'(4)-S_N'(5)-S_N'(6)</math>

avec <math>S_N'(a)=\sum_{k=0}^{N-1} \frac{16^{N-k}\bmod(8k+a)}{8k+a}+\sum_{k=N}^{N+P+10} \frac{16^{N-k}}{8k+a}</math>.

== Complexité de cette méthode ==
Pour calculer le ''n''-ième chiffre après la virgule de {{math|π}} en base 16 (et donc le 4''n''-ième chiffre en base 2) :

=== Complexité temporelle ===
* {{math|''B<sub>N</sub>{{'}}''(''a'')}} se calcule en temps constant ([[Comparaison asymptotique|O]](1)). La complexité du calcul de {{math|''S<sub>N</sub>{{'}}''}} est donc la même que la complexité du calcul de {{math|''A<sub>N</sub>{{'}}''(''a'')}}.
* {{math|''A<sub>N</sub>{{'}}''(''a'')}} : en utilisant la méthode d'exponentiation rapide, ses termes se calculent en O(log<sub>2</sub>(''n'')) multiplications sur des entiers de taille log<sub>2</sub>(''n''). En notant ''M''(''k'') la complexité de la multiplication de deux entiers de taille ''k'', la complexité est donc O(log(''n'')''M''(log(''n''))). Finalement, la somme des ''n'' termes, {{math|''A<sub>N</sub>{{'}}''(''a'')}}, se calcule en temps O(''n'' log(''n'')''M''(log(''n''))). Même en utilisant [[Algorithme de multiplication|l'algorithme de multiplication naïf]] plutôt que l'[[algorithme de Karatsuba]] ou la [[Transformation de Fourier rapide|transformée de Fourier rapide]], on obtient une complexité quasi linéaire de O(''n'' log(''n''){{3}}).

=== Complexité spatiale ===
Le calcul de {{math|''B<sub>N</sub>{{'}}''(''a'')}} s'effectue en [[complexité en espace|espace]] constant (somme d'un nombre fixé de termes, avec un nombre fixé de chiffres significatifs). Le calcul de {{math|''A<sub>N</sub>{{'}}''(''a'')}} nécessite d'effectuer des calculs modulo 8''k''+''a'', c'est-à-dire de manipuler des nombres de taille log(''k'') avec ''k'' ≤ ''N''. À chaque étape de l'algorithme, on manipule un nombre constant de tels nombres : la complexité en espace du calcul de {{math|''A<sub>N</sub>{{'}}''(''a'')}} est donc O(log(''n'')). L'algorithme total utilise donc un espace logarithmique.

== Formules dérivées ==
=== Simon Plouffe ===
Formule originale : <math>\pi\ =\ \sum_{k=0}^\infty \frac1{16^k}\left(\frac4{8k+1}-\frac2{8k+4}-\frac1{8k+5}-\frac1{8k+6}\right)</math>

:<math>\forall r \in\Complex\ \ \ \ \pi\ =\ \sum_{k=0}^\infty \frac1{16^k}\left(\frac{4+8r}{8k+1}-\frac{8r}{8k+2}-\frac{4r}{8k+3}-\frac{2+8r}{8k+4}-\frac{1+2r}{8k+5}-\frac{1+2r}{8k+6}+\frac r{8k+7}\right)</math>

:<math>\pi\sqrt2\ =\ \sum_{k=0}^\infty\frac{(-1)^k}{8^k}\left(\frac4{6k+1}+\frac1{6k+2}+\frac1{6k+3}\right)</math>

:<math>\pi^2\ =\ \sum_{k=0}^\infty\frac1{16^k}\left(\frac{16}{(8k+1)^2}-\frac{16}{(8k+2)^2}-\frac8{(8k+3)^2}-\frac{16}{(8k+4)^2}-\frac4{(8k+5)^2}-\frac4{(8k+6)^2}-\frac2{(8k+7))^2}\right)</math>

:<math>\pi^2\ =\ \frac98\ \sum_{k=0}^\infty\frac1{64^k}\left(\frac{16}{(6k+1)^2}-\frac{24}{(6k+2)^2}-\frac8{(6k+3)^2}-\frac6{(6k+4)^2}-\frac1{(6k+5)^2}\right)</math>

:<math>\pi^2\ =\ \frac2{27}\ \sum_{k=0}^\infty\frac1{729^k}\left(\frac{243}{(12k+1)^2}-\frac{405}{(12k+2)^2}-\frac{81}{(12k+4)^2}-\frac{27}{(12k+5)^2}-\frac{72}{(12k+6)^2}-\frac9{(12k+7)^2}-\frac9{(12k+8)^2}-\frac5{(12k+10)^2}+\frac1{(12k+11)^2}\right)</math>

=== Viktor Adamchick et Stan Wagon (1997) ===
:<math>\pi\ =\ \sum_{i=0}^\infty\frac{(-1)^i}{4^{i}}\left(\frac2{4i+1}+\frac2{4i+2}+\frac1{4i+3}\right)</math>.

=== [[Fabrice Bellard]] ===
:<math>\pi\ =\ \frac1{64}\ \sum_{i=0}^\infty \frac{(-1)^i}{2^{10i}}\left(-\frac{32}{4i+1}-\frac1{4i+3}+\frac{256}{10i+1}-\frac{64}{10i+3}-\frac4{10i+5}-\frac4{10i+7}+\frac1{10i+9}\right)</math>.

=== Géry Huvent (2001) ===
:<math> \pi\ =\ \frac1{128}\ \sum_{i=0}^\infty\frac1{2^{12i}}\left(\frac{768}{24i+3}+\frac{512}{24i+4}+\frac{128}{24i+6}-\frac{16}{24i+12}-\frac{16}{24i+14}-\frac{12}{24i+15}+\frac2{24i+20}-\frac1{24i+22}\right)</math>
:<math> \pi\ =\ \frac1{256}\ \sum_{i=0}^\infty\frac1{2^{12i}}\left(\frac{2048}{24i+2}-\frac{1024}{24i+4}+\frac{192}{24i+9}+\frac{128}{24i+10}+\frac{32}{24i+12}-\frac4{24i+18}-\frac4{24i+20}-\frac3{24i+21}\right)</math>
:<math> \pi^3\ =\ \frac1{16}\ \sum_{i=0}^\infty\frac{(-1)^i}{2^{10i}}\left(\frac{32}{(4i+1)^3}+\frac8{(4i+2)^3}+\frac1{(4i+3)^3}\right)</math> 
::<math>\ +\ \frac52\ \sum_{i=0}^\infty\frac{(-1)^i}{2^{6i}}\left(\frac{32}{(12i+1)^3}-\frac{192}{(12i+2)^3}+\frac{88}{(12i+3)^3}-\frac8{(12i+5)^3}+\frac{84}{(12i+6)^3}-\frac4{(12i+7)^3}+\frac{11}{(12i+9)^3}-\frac{12}{(12i+10)^3}+\frac1{(12i+11)^3}\right)</math>

:<math> \pi^4\ =\ \frac{27}{164}\ \sum_{i=0}^\infty\frac1{2^{12i}}\left(\frac{2048}{(24i+1)^4}-\frac{38912}{(24i+2)^4}+\frac{81920}{(24i+3)^4}-\frac{2048}{(24i+4)^4}-\frac{512}{(24i+5)^4}-\frac{23552}{(24i+6)^4}+\frac{256}{(24i+7)^4}-\frac{27648}{(24i+8)^4}-\frac{10240}{(24i+9)^4}\right.</math>
::<math>-\frac{2432}{(24i+10)^4}-\frac{64}{(24i+11)^4}-\frac{3584}{(24i+12)^4}-\frac{32}{(24i+13)^4}-\frac{608}{(24i+14)^4}-\frac{1280}{(24i+15)^4}-\frac{1728}{(24i+16)^4}+\frac8{(24i+17)^4}-\frac{368}{(24i+18)^4}-\frac4{(24i+19)^4}</math>
::<math>-\left.\frac8{(24i+20)^4}+\frac{160}{(24i+21)^4}-\frac{38}{(24i+22)^4}+\frac1{(24i+23)^4}\right)</math><ref>{{Lien web |auteur=Géry Huvent |titre=Formules BBP |url=https://c8cf35d7-ad1f-4031-84c6-bb15d656c101.filesusr.com/ugd/a97b48_b22c37a8b999457f85fe6d8528582dfa.pdf?index=true |année=2001}}.</ref>.

== Les records ==
Pour comparaison, le record de calcul de toutes les décimales de {{math|π}} est, en 2016, de 22 600 milliards de décimales (soit environ 70 000 milliards de chiffres binaires).

* {{Date|7|octobre|1996}} ([[Fabrice Bellard]]) : 400 milliardième chiffre en base 2
* septembre [[1997]] (Fabrice Bellard) : {{formatnum:1000}} milliardième chiffre en base 2
* février [[1999]] (Colin Percival) : {{formatnum:40000}} milliardième chiffre en base 2
* [[2001]] : {{formatnum:4000000}} milliardième chiffre en base 2

== Calcul en base 10 ==
Actuellement, aucune formule réellement efficace n'a été découverte pour calculer le ''n''-ième chiffre de {{math|π}} en base 10. [[Simon Plouffe]] a mis au point en décembre [[1996]], à partir d'une très ancienne série de calcul de {{math|π}} basée sur les coefficients du [[binôme de Newton]], une méthode pour calculer les chiffres en base 10, mais sa complexité en [[Comparaison asymptotique|O]](''n''{{3}} × log<sub>2</sub>(''n'')) la rendait en pratique inutilisable. [[Fabrice Bellard]] a bien amélioré l'algorithme pour atteindre une complexité en [[Comparaison asymptotique|O]](''n''{{2}}), mais cela n'est pas suffisant pour concurrencer les méthodes classiques de calcul de toutes les décimales.

== Références ==
{{Références}}

== Bibliographie ==
* {{Ouvrage |langue=fr |langue originale=de |auteur1=Jörg Arndt |auteur2=Christoph Haenel |traducteur=Henri Lemberg et François Guénard |titre=À la poursuite de π |titre original=Pi. Algorithmen, Computer, Arithmetik |éditeur=[[Vuibert]] |lieu=Paris |date=mars 2006 |année première édition=1998 |pages totales=273 |isbn=978-2-7117-7170-7 |numéro chapitre=10 |titre chapitre=L'algorithme BBP}}.

== Lien externe ==
{{MathWorld|nom_url=BBPFormula|titre=BBP Formula}}

{{Portail|mathématiques}}

[[Catégorie:Pi]]