Formule BBP

La formule BBP (ou Bailey-Borwein-Plouffe) permet de calculer le n^ème chiffre après la virgule 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 H. Bailey (en) et Peter Borwein (en).

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

Exploitation de la formule pour calculer les chiffres après la virgule de π

Le but est de calculer le N^e chiffre après la virgule de π en base 16.

Déjà, on remarque que le (N+1)^e chiffre après la virgule de π en base 16 est la même que le 1^re chiffre après la virgule de 16^Nπ. En effet, comme en base dix, multiplier un nombre 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 chiffre 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 :

1. D'abord, ce nombre étant très grand, cela demande d'effectuer des calculs sur des nombres très grands ;
2. 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 chiffres 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.

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

• En effet, le premier terme de la somme est : $b_N=\frac{1}{8N+a}$. Comme on cherche le N^e chiffre derrière la virgule 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.

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

Conclusion

Au final, pour obtenir les premiers chiffres de π en base 16 (ou 2), il faut calculer les premiers chiffres 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}$.

Complexité de cette méthode

Pour calculer le n^e chiffre après la virgule de π en base 16 (et donc le 4n^e chiffre en base 2):

Complexité temporelle

• B_n'(a) se calcule en temps constant (O(1)). La complexité du calcul de S_n' est donc la même que la complexité du calcul de A_n'(a).
• A_n'(a) : en utilisant la méthode d'exponentiation rapide, ses termes se calculent en O(log₂(n)) multiplications sur des entiers de taille log₂(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, A_n'(a), se calcule en temps O(n log(n)M(log(n))). Même en utilisant l'algorithme de multiplication naïf, on obtient une complexité quasi-linéaire de O(n log(n)³).

Complexité spatiale

Le calcul de B_n'(a) s'effectue en espace constant (somme d'un nombre fixé de termes, avec une nombre fixé de chiffres significatifs). Le calcul de A_n'(a) nécessite d'effectuer des calculs modulo 8k+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 A_n'(a) est donc O(log(n)). L'algorithme total utilise donc un espace logarithmique.

Formules dérivées

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+8r}{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)$

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

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

Géry Huvent (2001)

$\pi^3\ =\ \frac{1}{16}\ \sum_{i=0}^\infty\frac{(-1)^i}{2^{10i}}\left(\frac{32}{(4i+1)^3}+\frac{8}{(4i+2)^3}+\frac{1}{(4i+3)^3}\right)$ $\ +\ \frac{5}{2}\ \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}-\frac{8}{(12i+5)^3}+\frac{84}{(12i+6)^3}-\frac{4}{(12i+7)^3}+\frac{11}{(12i+9)^3}-\frac{12}{(12i+10)^3}+\frac{1}{(12i+11)^3}\right)$

$\pi^4\ =\ \frac{27}{164}\ \sum_{i=0}^\infty\frac{1}{2^{12i}}(\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}$ $-\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}+\frac{8}{(24i+17)^4}-\frac{368}{(24i+18)^4}-\frac{4}{(24i+19)^4}$ $-\frac{8}{(24i+20)^4}+\frac{160}{(24i+21)^4}-\frac{38}{(24i+22)^4}+\frac{1}{(24i+23)^4})$

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 chiffre en base 2
• septembre 1997 (Fabrice Bellard) : 1 000 milliardième chiffre en base 2
• février 1999 (Colin Percival) : 40 000 milliardième chiffre en base 2
• 2001 : 4 000 000 milliardième chiffre en base 2

Et pour le calcul des décimales ?

Actuellement, aucune formule réellement efficace n'a été découverte pour calculer le n^e chiffre 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.

Annexe : démonstration de la formule BBP

Notons $S_n=\sum_{k=0}^\infty \frac{1}{16^k(8k+n)}$ et démontrons la formule de Plouffe généralisée :

$(0)\quad \forall r \in \mathbb{C}\qquad\pi\ =\ (4+8r)S_1-8rS_2 -4rS_3 -(2+8r)S_4-(1+2r)S_5-(1+2r)S_6+rS_7$

(le cas r = 0 est sa formule originale ; le cas r = − 1 / 4 est, sous une forme plus détaillée, celle d'Adamchick et Wagon).

Posons :

$\alpha=1-i=\sqrt 2e^{-i\pi/4}$ et calculons de deux façons l'intégrale suivante :

$I=\int_0^1\frac{dy}{\alpha-y}.$

Elle est d'une part reliée aux S_n par

$\begin{align}I&=\frac 1\alpha\int_0^1\frac{dy}{1-y/\alpha}= \frac1\alpha\int_0^1\sum_{m=0}^\infty\frac{y^m}{\alpha^m}dy=\sum_{m=0}^\infty\frac{e^{i(m+1)\pi/4}}{(m+1)\sqrt 2^{m+1}}\\&= \frac{1+i}2S_1+\frac i2S_2+\frac{-1+i}4S_3-\frac14S_4-\frac{1+i}8S_5-\frac i8S_6+\frac{1-i}{16}S_7+\frac1{16}S_8,\end{align}$

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 :

$\begin{align}I&=-\left[\ln(\alpha-y)\right]_0^1=\ln\left(\frac\alpha{\alpha-1}\right)=\ln(1+i)=\ln(\sqrt 2e^{i\pi/4})=\frac{\ln 2}2+i\frac\pi4.\end{align}$

L'égalité entre ces deux expressions de I équivaut à :

$(1)\qquad\pi=4\mathop{Im}(I)=2S_1+2S_2+S_3-\frac 12S_5-\frac 12S_6-\frac 14S_7,$

$(2)\qquad\ln 2=2\mathop{Re}(I)= S_1-\frac 12S_3-\frac 12S_4-\frac 14S_5+\frac 18S_7+\frac 18S_8.$

Mais ln 2 s'exprime par ailleurs plus directement en fonction des S_n :

$\begin{align}(3)\quad\ln 2&=[-\ln(2-y^2)]_0^1= \int_0^1\frac{y}{1-y^2/2}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}$

(2) et (3) donnent donc, par soustraction des membres de droite, une relation entre les S_n :

$(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.$

En multipliant le membre de droite de (4) par 1 + 4r et en ajoutant ce produit au membre de droite de (1), on obtient l'égalité (0) annoncée.