« [Il est] raisonnable de spéculer que calculer le n-ième digit de pi n'est pas vraiment plus facile que calculer tous les digits jusqu'au n-ième » (Les frères Borwein, 1989)

Avec un peu de retard, il est temps de terminer ma tétralogie d'articles sur les décimales de pi (1-2-3), avec sans conteste la formule qui, le 19 septembre 1995, à 0h29 (heure locale), a révolutionné le monde des décimales de π. Cette formule ressemble à ça :

BBP
La formule BBP

Avant les années 90, calculer des décimales de π ressemblait toujours à peu près à la même chose : on  trouve une formule donnant π, on l'implémente dans son ordinateur (ou, avant les années 50, sur son papier) et, si la formule converge bien, on découvre par approximations successives plus ou moins rapidement des décimales exactes de π. L'arrivée de l'algorithme compte-gouttes, qui donne une à une les décimales de π a eu son petit effet, mais c'est 4 ans plus tard que la révolution est arrivée : il est possible de calculer isolément la 400 milliardième décimale, sans calculer toutes les précédentes !

Un petit mot sur le mathématicien canadien Simon Plouffe, le P de BBP (les deux autres étant David Bailey et Peter Borwein - le même qui disait 6 ans plus tôt que c'était impossible) : il est le créateur de l'inverseur de Plouffe (Après un calcul issu d'une modélisation bien compliquée  de la cinétique de l'oxydation du sulfite de cuivre (par exemple), vous trouvez la valeur 2.638794927335... Il suffit simplement de rentrer ces décimales dans l'inverseur qui, en donnera une écriture bien plus commode) ; il est également co-auteur, avec N. Sloane, de l'OEIS (L'encyclopédie des suites entières qui permet, par exemple, de trouver à quoi correspond la suite 1,11,21,1211 si on ne la connaît pas déjà), et il est également entré dans le livre des records pour avoir récité de tête  en 1977 les 4096 premières décimales de π (1)...

 

Calculer isolément les décimales de π...
La formule BBP permet donc, pour quelqu'un qui sait bien s'en servir, de calculer isolément des décimales de π. Pour être honnête, ce ne sont pas à proprement parler les décimales que l'on calcule, mais des digits. (en base 16, ou en base 2). Les décimales du nombre 23 sont "2-3". Écrit en binaire, ce nombre devient 10111, ses digits sont alors "1-0-1-1-1". On parle aussi de digits pour la base 16, qui peuvent se déduire facilement de ceux de la base 2 par groupement des digits par paquets de 4. En base 16, le nombre 23 se note "17".
Ainsi, on sait, et c'est Bellard (encore lui) qui l'a calculé en octobre 1996, que le 400 milliardième digit (en base 2) de π est un 0, suivi immédiatement par 1001110000111000. Le record a été battu depuis, et on connaît depuis 1999 quelle est le 40,000 milliardième digit de π.
Pour ce qui est du calcul isolé véritable des décimales (en base 10) de π, une formule similaire à celle en base 16 existe, mais l'utiliser pour calculer (par exemple) le millionième digit prend plus de temps que de calculer les 999 999 décimales précédentes...

 

La formule en elle-même n'est pas bien compliquée (elle demande de savoir intégrer des fonctions rationnelles, et quelque petites interversions séries/intégrales...), elle aurait parfaitement pu être découverte 250 ans plus tôt par Euler. La seule difficulté résidant dans l'idée qu'une telle formule puisse exister ! Elle a depuis été détournée pour donner des formules du même genre permettant le calcul des digits de π² ou de π√2...

Mais au fait... Comment utiliser cette formule pour calculer le 100(..)000e digit de π ? Cela repose sur 4 idées :

1ère idée : il est facile d'additionner localement deux grands entiers
(et de multiplier localement un grand entier par un petit)

Étant donné deux grands nombres, possédant par exemple 200 chiffres chacun (on peut ici raisonner avec des décimales), on est capable de dire avec précision quelle va être la 100e décimale de leur somme, sans avoir à calculer entièrement la somme. Il suffit de poser l'addition :

addition_locale

Voici localement ce que donne la somme de deux nombres N et N' de 200 chiffres. Les chiffres a,b,c représentent les 100e, 101e et 102e chiffres de N, respectivement pour N'. Ce qui nous intéresse ici, c'est de connaître le 100e chiffre, a'', de la somme. Pour cela, on peut simplement calculer la somme abc+a'b'c' (Une simple somme de nombres à 3 chiffres). On trouvera un nombre à 3 chiffres (éventuellement à 4 chiffres si a+a' donne une retenue), que l'on écrit abc+a'b'c'=(1?)def. Seuls ces 3 derniers chiffres nous intéressent.
Le 100e chiffre a'' de N+N' a alors toutes les chances de correspondre au chiffre d de la somme abc+a'b'c'. Il ne peut arriver qu'une mauvaise surprise, c'est que b+b'=9 et c+c'=9 : on ne sait pas si la somme des 104e chiffres de N et N' amènera une retenue qui changera d. Un tel évènement est rare, et si il arrive, on peut le détecter, et corriger la valeur de a'' si besoin est.

Finalement, pour calculer le n-ième chiffre d'une somme de deux grands entier, il suffit simplement de faire la somme à partir de la n-ième en poussant un peu vers la droite pour éviter les erreurs de retenues. On peut étendre l'idée à la somme de plusieurs grands entiers et à la multiplication d'un grand entier par un plus petit.

2eme idée : On peut calculer le n-ième digit d'un nombre de la forme 1/(k.16i) sans trop de calculs.

Pour simplifier, on va raisonner en base 10 : l'énoncé dit qu'il est facile de calculer la n-ième décimale de 1/(k.10i). Prenons par exemple n=1000, i=42 et k=49 : on cherche la 1000e décimale de 1/(49.1042).
En fait, la 100e décimale de  1/(49.1042) n'est autre que la 999e décimale de  1/(49.1041) ou la 99e décimale de 1/(49.1040) (Multiplier par 10 revient à décaler les décimales). Au final, on est amené à chercher la 958e décimale de 1/49.
En continuant sur le principe de décalage par multiplication par 10, on peut dire que la 958e décimale de 1/49 est la 957e décimale de 10/49, et ainsi de suite : c'est la 1ère décimale (juste après la virgule) de P=10957/49.

Maintenant, supposons que l'on puisse facilement faire la division euclidienne (celle que l'on pose à l'école, avec le quotient et le reste) de 10957 par 49 : on a 10957=49q+r, avec r<49. On aura alors  P=10957/49 = q+r/49. Le chiffre juste après la virgule de P est le même que celui juste après la virgule de r/49. Puisque r est petit, la division r/49 sera facile à faire, et ce chiffre facile à trouver.

Faire la division euclidienne de 10957 par 49 est l'objet de l'idée suivante :

3eme idée : On peut facilement trouver le reste dans la division euclidienne de 16n-i-1 par k sans manipuler de grands nombres.

Ramenons-nous encore une fois dans le cas décimal. On cherche donc le reste r dans la division euclidienne de 10957 par 49. Autrement dit, on cherche combien vaut 10957 modulo 49 (c'est simplement une affaire de vocabulaire, ça veut dire la même chose). On raisonnera alors dans l'arithmétique modulo 49, où les seuls nombres permis sont les nombres entre 0 et 48. Si un nombre est plus grand, on retranche 49 autant de fois qu'il sera nécessaire pour tomber dans l'intervalle [0,49]. Par exemple, dans cette arithmétique, 104=6, puisque 104-49-49=6. Dans cette arithmétique, les additions, soustractions et multiplication se comportent comme dans l'arithmétique habituelle.

On commence par calculer 10s dans l'arithmétique modulo 49 pour les puissances de 2, en procédant de proche en proche :

101 = 10         102 = 100 = 2         104 = 22 = 4         108 = 42 = 16
1016 = 162 = 256 = 11         1032 = 112 = 121 = 23         1064 = 232 = 529 = 39
10128 = 392 = 1521 = 2         10256 = 22 = 4         10512 = 42 = 16

On peut alors écrire 957 comme somme de puissance de 2 :
957 = 512+256+128+32+16+8+4+1
On a alors le résultat recherché :
10957 = 10512.10256.10128.1032.1016.108.104.101
        =16.4.2.23.11.16.4.10 = 20725760 = 34

Cet algorithme, permettant de calculer le reste de la division de AN  (N grand) par k (petit) en passant par la décomposition en puissance de 2, ne se limite pas à son utilisation pour la formule BBP, mais est utilisé dans tous les algorithmes d'arithmétique (par exemple, dans l'algorithme RSA, qui s'occupe de crypter vos coordonnées bancaires). Même si N est très grand, le calcul sera tout de même très rapide.

4eme idée : la formule BBP permet de calculer le n-ième digit de π en ne travaillant que sur des petits nombres

La formule BBP, en développant un peu, s'écrit :

BBpd

Grâce à l'idée 2, trouver le n-ième digit (et les quelques uns qui suivent) de chaque terme de la somme se fait sans manipuler de trop grands entiers.
On va faire ça pour tous les termes de la somme entre k=0 et N un peu plus grand que n. En choisissant bien le "un peu plus grand que", les termes de la série entre N et ∞ seront trop petits pour apporter des erreurs de retenues (Cela vient du fait que 1/16k décroit très vite).
L'idée 1 permettra alors de trouver le n-ième digit tant convoité, et les quelques qui suivent, en additionnant tous les digits de 0 à N...

Deux autres formules du même genre, signé Simon Plouffe :

Plouffe1

Plouffe2


(1) À propos des records de mémorisations des décimale du nombre π, il faut tout de même dire que Plouffe s'est arrêté volontairement à la 4096e décimale (alors qu'il en avait appris 4400, mais 4096 est un nombre plus joli).
Ce record a été battu deux ans plus tard par l'anglais M. Poultney, qui en a mémorisé 5050.
Le record actuel est détenu par le Chinois Chao Lu, qui en a mémorisé 67890. Sa récitation a duré un peu plus de 24h...
Parmi les records, on peut également citer L'anglais D. Turner, qui a récité 1250 décimales de π à l'endroit, puis à l'envers. Il y a aussi le suédois Mats Bergsten qui a donné les 9778 premières décimales tout en jonglant avec 3 balles, sans oublier la canadienne Grace Hare, qui en a récité 31, mais qui n'avait alors que 3 ans... (voir là-bas)


Sources :
Jean-Paul Delahaye, Le fascinant nombre π, Éditions Belin, Pour la Science