Formule BBP

Un article de Wikipédia, l'encyclopédie libre.

La formule BBP (ou Bailey-Borwein-Plouffe) permet de calculer l'énième décimale de π en base 2 (ou 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 19 septembre 1995 par Simon Plouffe en collaboration avec David Bailey et Peter Borwein

Sommaire

1 La formule
2 Exploitation de la formule pour calculer les décimales de π
3 Complexité de cette méthode
- 3.1 Complexité temporelle
- 3.2 Complexité spatiale
4 Formules dérivées
5 Les records
6 Et pour le calcul des décimales ?
7 Annexe : démonstration de la formule BBP

[modifier] La formule

$\pi = \sum_{k=0}^\infty \frac{1}{16^k}\left(\frac{4}{8k+1}-\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6}\right)$

[modifier] Exploitation de la formule pour calculer les décimales de π

Le but est de calculer le N^e digit de π en base 16.

Déjà, on remarque que la (N+1)^e décimale de π en base 16 est le même que la 1^re décimale de 16^Nπ. En effet, comme en base dix, multiplier un chiffre en base 16 par 16 permet de décaler la virgule d'un rang vers la droite. Donc en multipliant un nombre par 16^N, la virgule est décalée de N rang vers la droite. Il « suffit » donc de calculer le premier digit de 16^Nπ, égal par la formule BBP à:

$16^{N}\pi = \sum_{k=0}^\infty 16^{N-k}\left(\frac{4}{8k+1}-\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6}\right)$

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 $S_N(a)=\sum_{k=0}^\infty \frac{16^{N-k}}{8k+a}$ . Le calcul des premiers digits de S_N(a) permettra d'obtenir ceux de 16^Nπ, par la relation :

$16^{N}\pi = 4\ S_N(1) - 2 S_N(4) - S_N(5) - S_N(6)$

Découpons la somme S_N(a) en deux :

$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)$

et calculons A_N(a) et B_N(a) indépendamment.

[modifier] Calcul de B_N(a)

$B_N(a)=\sum_{k=N}^ \infty \frac{16^{N-k}}{8k+a}$

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 digits.

En effet, le premier terme de la somme est : $b_N=\frac{1}{8N+a}$ . Comme on cherche la N^e décimale de π (N = 1 000 000 000 par exemple), le premier terme b_N 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_k > 16 b_k+1 :

$\frac{b_k}{b_{k+1}}=\frac{16^{N-k}}{16^{N-(k+1)}}\frac{8(k+1)+a}{8k+a}=16\left( 1+\frac{8}{8k+a}\right) \longrightarrow 16^{+}$ .

Finalement, la somme B_N(a) est de la forme (au pire) :

$B_N\ =\ \ 0,**********.\ .\ .\ .\$

$+\ 0,0*********.\ .\ .\$

$+\ 0,00********.\ .\ .\$

$+\ 0,000********.\ .\ .\$

Donc pour obtenir B_N(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 : $B_N'(a)=\sum_{k=N}^{N+P+10} \frac{16^{N-k}}{8k+a}$

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.

[modifier] Calcul de A_N(a)

$A_N(a)=\sum_{k=0}^{N-1} \frac{16^{N-k}}{8k+a}$

Le problème pour calculer A_N(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, aussi grande qu'elle soit. On peut donc s'en « débarrasser » en utilisant arithmétique modulo.

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

Pour cela, on effectue la division euclidienne de 16^N-k par 8k+a :

$\exists q\in \mathbb{Z}, \exists r < 8k+a,\ 16^{N-k}=q(8k+a)+r$

Donc $\frac{16^{N-k}}{8k+a}=q+\frac{r}{8k+a}$

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

Et $\frac{r}{8k+a}=\frac{16^{N-k}[8k+a]}{8k+a}$

Il suffit donc de calculer : $A_N'(a)=\sum_{k=0}^{N-1} \frac{16^{N-k}[8k+a]}{8k+a}$ .

En utilisant la méthode d'exponentiation rapide, 16^N-k[8k+a] se calcule rapidement (temps d'exécution en O(log₂(N-k)).

[modifier] Conclusion

Au final, pour obtenir les premiers digits de π en base 16 (ou 2), il faut calculer les premiers digits de :

$\pi_N=4\ S_N'(1)-2\ S_N'(4)-S_N'(5)-S_N'(6)$

avec $S_N'(a)=\sum_{k=0}^{N-1} \frac{16^{N-k}[8k+a]}{8k+a}+\sum_{k=N}^{N+P+10} \frac{16^{N-k}}{8k+a}$ .

[modifier] Complexité de cette méthode

Pour calculer le n^e digit de π en base 16 (et donc le 4n^e digit en base 2):

[modifier] Complexité temporelle

B_n'(a) se calcule en complexité linéaire (O(1))
A_n'(a) : en utilisant la méthode d'exponentiation rapide, ses termes se calculent en O(log₂(n)). Donc la somme des n termes, A_n'(a), se calcule en O(n*log₂(n)).

Donc S_n'(a) se calcule en O(1)+O(n*log₂(n))=O(n*log₂(n)). Donc finalement, π_n se calcule en 4*O(n*log₂(n))=O(n*log₂(n)).

D'où un temps des calcul proportionnel à n*log₂(n), soit quasi-linéaire.

[modifier] Complexité spatiale

La complexité en utilisation mémoire est constante : seules des sommes successives sur de petits nombres sont effectuées (une précision d'une dizaine de décimales est nécessaire).

[modifier] Formules dérivées

[modifier] Simon Plouffe

Formule originale : $\pi\ =\ \sum_{i=0}^\infty \frac{1}{16^i}\left(\frac{4}{8i+1}-\frac{2}{8i+4}-\frac{1}{8i+5}-\frac{1}{8i+6}\right)$

$\forall r \in \mathbb{C}\ \ \ \ \pi\ =\ \sum_{i=0}^\infty \frac{1}{16^i}\left(\frac{4+8r}{8i+1}-\frac{8r}{8i+2}-\frac{4r}{8i+3}-\frac{2+8i}{8i+4}-\frac{1+2r}{8i+5}-\frac{1+2r}{8i+6}+\frac{r}{8i+7}\right)$

$\pi\sqrt{2}\ =\ \sum_{i=0}^\infty \frac{(-1)^i}{8^i}\left(\frac{4}{6i+1}+\frac{1}{6i+2}+\frac{1}{6i+3}\right)$

$\pi^2\ =\ \sum_{i=0}^\infty \frac{1}{16^i}\left(\frac{16}{(8i+1)^2}-\frac{16}{(8i+2)^2}-\frac{8}{(8i+3)^2}-\frac{16}{(8i+4)^2}-\frac{4}{(8i+5)^2}-\frac{4}{(8i+6)^2}-\frac{2}{(8i+7))^2}\right)$

$\pi^2\ =\ \frac{9}{8}\ \sum_{i=0}^\infty \frac{1}{64^i}\left(\frac{16}{(6i+1)^2}-\frac{24}{(6i+2)^2}-\frac{8}{(6i+3)^2}-\frac{6}{(6i+4)^2}-\frac{1}{(6i+5)^2}\right)$

$\pi^2\ =\ \frac{2}{27}\ \sum_{i=0}^\infty \frac{1}{729^i}\left(\frac{243}{(12i+1)^2}-\frac{405}{(12i+2)^2}-\frac{81}{(12i+4)^2}-\frac{27}{(12i+5)^2}-\frac{72}{(12i+6)^2}-\frac{9}{(12i+7)^2}-\frac{9}{(12i+8)^2}-\frac{5}{(12i+10)^2}+\frac{1}{(12i+11)^2}\right)$

[modifier] Viktor Adamchick et Stan Wagon (1997)

$\pi\ =\ \sum_{i=0}^\infty \frac{(-1)^i}{4^{i}}\left(\frac{2}{4i+1}+\frac{2}{4i+2}+\frac{1}{4i+3}\right)$

[modifier] Fabrice Bellard

$\pi\ =\ \frac{1}{64}\ \sum_{i=0}^\infty \frac{(-1)^i}{2^{10i}}\left(-\frac{32}{4i+1}-\frac{1}{4i+3}+\frac{256}{10i+1}-\frac{64}{10i+3}-\frac{4}{10i+5}-\frac{4}{10i+7}+\frac{1}{10i+9}\right)$

[modifier] Les records

Pour comparaison, le record actuel de calcul de toutes les décimales de π est de 1 241 milliards de décimales (soit environ 4 123 milliard de chiffres binaires).

7 octobre 1996 (Fabrice Bellard) : 400 milliardième digit en base 2
septembre 1997 (Fabrice Bellard) : 1 000 milliardième digit en base 2
février 1999 (Colin Percival) : 40 000 milliardième digit en base 2
2001 : 4 000 000 milliardième digit en base 2

[modifier] Et pour le calcul des décimales ?

Actuellement, aucune formule réellement efficace n'a été découverte pour calculer le n^e digit de π en base 10. Simon Plouffe a mis au point en décembre 1996, à partir d'une très ancienne série de calcul de π basée sur les coefficients du binôme de Newton, une méthode pour calculer les chiffres en base 10, mais sa complexité en O(n³*log₂(n)) la rendait en pratique inutilisable. Fabrice Bellard a bien amélioré l'algorithme pour atteindre une complexité en O(n²), mais cela n'est pas suffisant pour concurrencer les méthodes classiques de calcul de toutes les décimales.

[modifier] Annexe : démonstration de la formule BBP

Démontrons déjà un résultat général :

$\forall n \in \mathbb{N}: \sum_{k=0}^\infty \frac{1}{16^k(8k+n)}=\sqrt{2}^{n}\sum_{k=0}^\infty \frac{\left(\frac{1}{\sqrt{2}}\right) ^{8k+n}}{8k+n}$

$=\sqrt{2}^{n}\sum_{k=0}^\infty \left[\frac{x^{8k+n}}{8k+n}\right]_0^{\frac{1}{\sqrt{2}}}$

$=\sqrt{2}^{n}\sum_{k=0}^\infty \left(\int_{0}^{\frac{1}{\sqrt{2}}} x^{n+8k-1}dx \right)$

$=\sqrt{2}^{n}\int_{0}^{\frac{1}{\sqrt{2}}} \left(\sum_{k=0}^\infty x^{n+8k-1}\right)dx$

$=\sqrt{2}^{n}\int_{0}^{\frac{1}{\sqrt{2}}} \bigg(x^{n-1}\underbrace{\sum_{k=0}^\infty \left(x^8\right)^k} \bigg)dx$

$=\sqrt{2}^{n}\int_{0}^{\frac{1}{\sqrt{2}}} \bigg(x^{n-1} \frac{1}{1-x^8}\bigg)dx$

d'où : $\sum_{k=0}^\infty \frac{1}{16^k(8k+n)}=\sqrt{2}^{n}\int_{0}^{\frac{1}{\sqrt{2}}} \frac{x^{n-1}}{1-x^8}dx$

Appliquons ce résultat à la formule BBP:

$\sum_{k=0}^\infty \frac{1}{16^k}\left(\frac{4}{8k+1}-\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6}\right)$

$= 4\sum_{k=0}^\infty \frac{1}{16^k(8k+1)}-2\sum_{k=0}^\infty \frac{1}{16^k(8k+4)}-\sum_{k=0}^\infty \frac{1}{16^k(8k+5)}-\sum_{k=0}^\infty \frac{1}{16^k(8k+6)}$

$= 4\left(\sqrt{2}^{1}\int_{0}^{\frac{1}{\sqrt{2}}} \frac{x^{1-1}}{1-x^8}dx \right)-2\left(\sqrt{2}^{4}\int_{0}^{\frac{1}{\sqrt{2}}} \frac{x^{4-1}}{1-x^8}dx \right)-\left(\sqrt{2}^{5}\int_{0}^{\frac{1}{\sqrt{2}}} \frac{x^{5-1}}{1-x^8}dx \right)-\left(\sqrt{2}^{6}\int_{0}^{\frac{1}{\sqrt{2}}} \frac{x^{6-1}}{1-x^8}dx \right)$

$= \int_{0}^{\frac{1}{\sqrt{2}}} \frac{4\sqrt{2}-8x^3-4\sqrt{2}x^4-8x^5}{1-x^8}dx$

On pose $y=\sqrt{2}x$ :

$= \int_{0}^{1} \left(\frac{4\sqrt{2}-\frac{8}{\sqrt{2}^3}y^3-\frac{4\sqrt{2}}{\sqrt{2}^4}y^4-\frac{8}{\sqrt{2}^5}y^5}{1-\frac{1}{\sqrt{2}^8}y^8}\right)\frac{dy}{\sqrt{2}}$