www.pi314.net


Histoire
Mathématiciens
Toutes les formules
Approx. numériques
Programmes
Algos perso/divers
Décimales
Poèmes
Articles/vidéos
Délires !
 Pi-Day
Images/Fonds
Musique
Liens sur Pi
Bibliographie



Boris Gourévitch
L'univers de Pi - V2.57
modif. 13/04/2013



Pi-Day dans
Pi Day Countdown
Google
Accueil Historique/Actu (Pi, site, moi) Edito Livre d'or Pages en .pdf Je me présente Quelques photos Remerciements Page des nets d'or Sites qui m'indexent Derniers changements Contact

Cette page en français This page in English


pict

Simon Plouffe

Atteindre le N  -ième digit d’un nombre a  sans connaitre les précédents grâce aux formule BBP



1 La formule BBP

Petit rappel: le 19 septembre 1995 à 0h29 (!), après des mois de recherches à tâtons, Simon Plouffe, David Bailey et Peter Borwein de Vancouver découvrent la formule apparemment simple et anodine
     sum  oo  1 (   4       2       1       1  )
p =    16k   8k-+-1-  8k+-4-- 8k+-5-- 8k+-6-
    k=0
(1)

appelée depuis formule BBP (voir page Plouffe). Nos trois mathématiciens canadiens remarquent que cette expression est très proche d’une décomposition en base 16 de p  . Dans l’article [1] de 1996 devenu célèbre, ils montrent que l’on peut se servir de 1 pour atteindre le n  -ième digit de p  en base  p
2  (a fortiori en base 16) sans avoir calculé les précédents, ceci en temps presque linéaire        3
O(n.log n  ) et en espace O(logn)  !! L’espace requis étant minime, ils appliquent immédiatement ce résultat en calculant le 40  milliardième digit de p  en base 2.

Arrêtons-nous maintenant sur le principe de ce calcul du n  -ième digit de p  en considérant le cas général d’un nombre a  et d’une formule BBP associée:

2 Atteindre le N  -ième digit d’un nombre a  sans connaitre les précédents

Nous suivons ici à quelques adaptations près l’exposé remarquablement clair de Xavier Gourdon et Pascal Sebah sur le sujet [2] qui simplifiait déjà l’article original [1]. X. Gourdon est un grand connaisseur des algorithmes de calcul puisqu’il détient depuis plusieurs années le record de vitesse du calcul des décimales de p  avec le programme PiFast (plus rapide même que le programme de l’équipe de Kanada).

Pour simplifier le propos, on se restreint à la base 2  , et le N  -ième bit de a  désigne le N  -ième bit de la partie fractionnaire de a  , c’est-à-dire après la virgule. La partie fractionnaire d’un nombre x  sera notée {x}=  x- [x]  .

2.1 Décalage

La première idée repose sur le théorème suivant

Theorem 2.1 Le N + n  -ième bit de a  s’obtient en calculant le n  -ième bit de 2Na  i.e. le n  -ième bit de {    }
 2N a .

Proof. La décomposition en base 2  de a  s’écrit

a = [a]+ {a}=  sum   ak+  sum   ak,        (a   (-  {0,1})
                  2k      2k           k
              k<0     k>0
(2)

d’où { 2Na} =  sum     -ak-=  sum     aN+k
           k>N 2k- N    k>0  2k  . Le n  -ième bit de {2N a} est donc a
  N+n  qui est bien le N + n  -ième bit de a  .  __

Example 2.2 Cherchons le 13ème bit de p  . Comme on a 212p = 12867.9635...  , en particulier, {    }
 212p  = 0.9635...  . Il reste donc à calculer le 1er  bit de 0.9635  . En passant en binaire, on a

         -1   1-  -1   1-   0-  -1   1-
0.9635...= 21 + 22 + 23 + 24 + 25 + 26 + 27 + ...= [0.1111011...]base 2
(3)

donc le 13`eme  bit de p  est 1  , et le 17`eme  bit est 0  .

Chercher le N + 1  -ième bit de p  correspond donc à trouver le premier bit de {    }
 2N p . D’après la formule BBP 1, cela correspond au premier bit de la série

       oo  ({   N-4n}   {   N- 4n }  {  N- 4n }  {  N- 4n } )
pN =  sum     4.2----  -   2.2-----  -  2-----  -  2-----   = AN + BN
     n=0    8n+ 1       8n + 4      8n + 5     8n + 6
(4)

AN  et BN  sont les sommes

         sum                   sum 
AN  =       ,        BN =      .
      0<n<N/4              n>N/4
(5)

Compte tenu de la décroissance rapide de 2N -4n  pour n > N-
    4  , il suffit d’évaluer seulement les premiers termes de BN  (en virgule flottante !), plus exactement suffisamment pour que le reste de la somme soit inférieur à la précision requise (ici 1  bit après la virgule).

Concernant AN  , on peut mettre son terme général sous la forme

{    }
 p.2n
  m
(6)

avec p, n  et m  dans N  . En écrivant p.2n = qm + r  la division euclidienne de p.2n  par m  , on voit que comme r < m  , { p.2n}   r-
   m   = m  . Ceci s’écrit encore

{ p.2n }   p.2n (mod m)
  -m--  = -----m-----.
(7)

2.2 Exponentiation

Imaginons ainsi que l’on veuille calculer le 1000`eme  bit de 231549  . Grâce aux principes précédents, cela revient à calculer le 965`eme  bit de 419  et donc le 1er  bit de 249694  . En vertu de 7, il faut alors trouver r  tel que 2964 = 49q + r  donc évaluer 2964 (mod 49)  .

Le calcul de 2964 (mod 49)  s’effectue au moyen de la méthode appelée exponentiation binaire (“binary powering” ou “binary exponentiation”), dont l’idée semble avoir été utilisée dès 200 avant J.C [3], et même envisagée par les Egyptiens pour effectuer la multiplication ! Cela revient simplement à utiliser l’arithmétique modulo 49  , autrement dit à se placer dans Z/49Z  , en trois étapes :

  1. On décompose 964  en puissances de 2  : 989 = 512 +256 + 128+ 64+ 4
  2. On calcule de proche en proche chaque puissance de 2  jusqu’à 512  : 22 = 4; 24 = (22)2 = 4*4 = 16;  28 = 162 = 256 = 49 *5+ 11 = 11; 216 = 112 = 121 = 23;                                                  232 = 232 = 529 = 39;  264 = 392 = 1521 = 2;  2128 = 22 = 4;  2256 = 16;   2512 = 162 = 11
  3. On utilise les étapes 1 et 2 : 2964 = 2512+256+128+64+4 = 11 ×16 × 4× 2× 4 = 5632 = 46  .

Pour obtenir plus généralement       n
r = p.2 (mod m)  , on initialise r  à 1 et on prend t  la plus grande puissance de 2  inférieure à n  . Puis on applique l’algorithme suivant équivalent à notre exemple :

  1. Si n > t  , alors r = 2.r (mod m)  ;  n = n- t  ; FinSi
  2. t = t/2;
  3. Si t > 1  alors r = r2 (mod m)  ; Aller à l’étape 1. ; FinSi.

Enfin on termine par r = p.r (mod m);  pour prendre en compte p  . En raison de la décomposition de n  en puissance de 2  , on utilise O(logn)  (mod m)  opérations pour ce calcul, ce qui est très faible !

2.3 Etapes finales

Il reste à obtenir r-= p.2n-(modm)
m       m  et à sommer le tout sur 0 < n < N-
        4  pour obtenir AN  dans 5. Il faut ici remarquer que si le calcul de chaque terme de AN  est effectué en virgule flottante et que, disons, le dernier bit est erroné, alors au pire la propagation des retenues sur toute la somme produira une erreur sur 1+ log2 N
       4  bits. Pour des valeurs de N  non extravagantes, on pourra donc mener les calculs de -r
m  en toute sérénité avec une précision arithmétique de 64  bits seulement.

Pour cette raison et en remarquant que pour AN  , on somme O(N )  termes nécessitant chacun O(logN )  opérations grâce à l’exponentiation, l’algorithme requiert O(N.logN.M (logN ))  opérations et O(logN )  en espace, où j < M (j) < j2  est la complexité de la multiplication de deux entiers de taille j  bits.

References

[1]   D.H. Bailey, P.B. Borwein, S. Plouffe, “On the Rapid Computation of Various Polylogarithmic Constants”, Mathematics of Computation, 1997, vol. 66, pp. 903-913.

[2]   X. Gourdon, P. Sebah, “N-th digit computation”, preprint, 2003, http://numbers.computation.free.fr/Constants/Algorithms/nthdigit.pdf.

[3]   D. E. Knuth, “The Art of Computer Programming. Vol. 2: Seminumerical Algorithms”, Addison-Wesley, Reading, MA, 1981.


Retour à la page d'accueil