Chapitre 4  Pour trouver les premières décimales de π

4.1  La formule de Gregory et le développement en série de arctan(x)

La série de Gregory est le développement en série entière de arctan(x)
On a :

Le reste de cette série alternée est du signe du premier terme négligé et est majoré en valeur absolue par la valeur absolue du premier terme négligé.

4.2  La formule de Machin

De la formule d’addition des tangentes à savoir :

on en déduit la formule pour ab<1 :

Exercices :
1/ Pour a=1/5 on cherche b pour que :

On doit donc résoudre :

On trouve :

Donc :

2/ Pour a=1/5 on cherche b pour que :

On doit donc résoudre :

On trouve :

Donc :

3/ En faisant encore deux fois le même genre de substitution, montrer la formule de Machin :

Pour cela on prend successivement :
c=7/17 et on trouve b=c−a/1+ac=9/46 et
c=9/46 et on trouve b=c−a/1+ac=−1/239
4/ Écrire un programme qui prend en entrée a et n et qui renvoie la liste L des valeurs de b_k vérifiant :
karctan(a)+arctan(b_k)=π/4 pour k=1..n (L[k−1]=b_k)
On tape le programme :

machin1(a,n):={
local k,c,L;
c:=1;
L:=[0];
for (k:=1;k<n+1;k:=k+1) {
c:=(c-a)/(1+a*c);
L:=append(L,c);
}
return(L);
};

On tape :
machin1(1/5,5)
On obtient :
[0,2/3,7/17,9/46,1/-239,-122/597]
Ainsi, pour a=1/5 et k=3, b_k=9/46 donc :
3arctan(1/5)+arctan(9/46)=π/4
On tape :
machin1(1/3,5)
On obtient :
[0,1/2,1/7,-2/11,-17/31,-41/38]
Ainsi, pour a=1/3 et k=2, b_k=1/7 donc :
2arctan(1/3)+arctan(1/7)=π/4

4.3  Les décimales de π avec les formules précédentes

4.3.1  Une remarque

Si on utilise pour calculer π la formule :

Si on prend x=1, on a arctan(1)=π/4=1−1/3+1/5+...+(−1)^k1/2k+1+... mais la convergence est lente.
On remarque que la convergence est beaucoup plus rapide pour x=1/5 et encore plus rapide pour x=1/239 d’où l’utilisation de la formule :

4.3.2  Le programme avec Xcas

Pour calculer la somme de n termes de la série on va utiliser la méthode de Hörner pour faire le moins possibles de multiplications, on a :
arctan(x)=x(1−x²(1/3−x²(1/5−x²(....−x²(1/2n−1)))))
L’utilisateur doit rentrer la valeur de x (a) et le nombre n de termes de la série.

//approx de arctan(a) par sa serie :
// a-a^3/3+..+(-1)^n*a^(2n+1)/(2n+1) 
//cette valeur est calculee par la methode de Horner
gregory(a,n):={
local t,k;
t:=1/(2*n-1);
for (k:=2*n-3;k>0;k:=k-2) {
t:=1/k-a^2*t;
}
return (a*t);
};

On tape :
16*gregory(1/5,18)-4*gregory(1/239,6)
On obtient :
3.14159265359 On tape :
evalf(16*gregory(1/5,42)-4*gregory(1/239,20))
On obtient, si on a choisit 60 Chiffres dans la configuration du CAS (menu Cfg->Configuration du CAS) :
3.14159265358979323846264338327950288419716939937510582097494

4.3.3  Combien faut-il calculer de termes ?

Soit R_n(x) le reste de la série ∑_k=1^∞(−1)^k+1x^2k−1/2k−1 : R_n(x)=∑_k=n+1^∞(−1)^k+1x^2k−1/2k−1.
On sait que |R_n(x)|<|x|²ⁿ⁺¹/2n+1 Pour avoir |R_n(1/5)|=|R_n(0.2)|<10⁻⁶¹ il faut que :
2²ⁿ⁺¹<(2n+1)10²ⁿ⁻⁶⁰ et comme 2¹0≃ 10³ cela donne si on suppose 2n+1>10 :
10^(6n+3)/10<10²ⁿ⁻⁵⁹ soit 593<14n soit n≃ 42 On vérifie pour n=42 on a (1/5)⁸⁵/85<2.56e−62
On peut aussi écrire si on suppose que 2n+1>10 :
|x|²ⁿ⁺¹/2n+1<x²ⁿ⁺¹/10<10⁻⁶¹ donc on va choisir x²ⁿ⁺¹<10⁻⁶⁰ soit (2n+1)log10(x)<−60 ou encore n>((−60)/log10(x)−1)/2 Pour x=1/5 on a n>42.4202967422 et comme 2n+1>40 on peut améliorer la majoration |x|²ⁿ⁺¹/2n+1<x²ⁿ⁺¹/40<10⁻⁶¹
ce qui donne n>((−60+log10(4))/log10(1/5)−1)/2=41.9896201841 donc on prend n=42
pour x=1/293 on a n>11.66 donc on prend n=12 et on vérifie : (1/239)²⁵/25=1.38711499837e−61.
On choisit 62 Chiffres dans la configuration du CAS (menu Cfg->Configuration du CAS) et on tape :
evalf(16*gregory(1/5,42)-4*gregory(1/239,12))
On obtient :
3.1415926535897932384626433832795028841971693993751058209749446
On tape :
evalf(pi)
On obtient :
3.1415926535897932384626433832795028841971693993751058209749446
Remarque
Avec cette méthode John Machin calcula 100 décimales de π en 1706.

4.3.4  Les formules de même type que celles de Machin

En 1973, Jean Guilloud a mis une journée pour calculer 10⁶ décimales de π en utilisant une formule de même type à savoir :

en vérifiant ses calculs avec la formule analogue :

En 1999, Yasumata Kanadaa a atteint le record en calculant 12411* 10⁸ décimales de π en utilisant une formule de même type à savoir :

et la formule analogue :

On peut vérifier ces formules avec Xcas, on tape par exemple :
tsimplify(12*atan(1/49)+32*atan(1/57)-5*atan(1/239)+12*atan(1/110443)) On obtient :
π/4
Comment trouver des formules de type Machin ?
voir aussi 5.5
Montrons pour cela que si a ∈ ℕ^* et si a²+1=a₁*a₂ avec (a₁,a₂) ∈ ℤ² et (a+a₁)(a+a−2)¬ 0 alors on a :
arctan(1/a)=arctan(1/(a+a₁))+arctan(1/(a+a₂)) (si a≥ 2 alors a ne divise pas a²+1 donc a+a−1 et a+a₂ sont non nuls et différents).
On a si xy<1, arctan(x)+arctan(y)=arctan((x+y)/(1−xy)) donc
arctan(1/(a+a₁))+arctan(1/(a+a₂))=arctan((2a+a₁+a₂)/((a+a₁)(a+a₂)−1))=arctan(1/a) puisque

• si a+a₁>1 et a+a₂>1 alors 1/(a+a₁)*1/(a+a₂)<1
• si a+a₁>1 et a+a₂<−1 alors 1/(a+a₁)*1/(a+a₂)<0<1
• si a+a₁≤ −1 et a+a₂<−1 alors 1/(a+a₁)*1/(a+a₂)=1/(−a−a₁)*1/(−a−a₂)<1

Donc :
si a²+1=a₁*a₂ alors arctan(1/(a+a₁))+arctan(1/(a+a₂))=arctan(1/a)
Exemples
si a=1 on a arctan(1)=arctan(1/2)+arctan(1/3).
si a=2 on a arctan(1/2)=arctan(1/3)+arctan(1/7) et
arctan(1/2)=arctan(1/(2−1))+arctan(1/(2−5))=arctan(1)−arctan(1/3) On a donc :
a=1 a²+1=2=1*2 donc π/4=arctan(1)=arctan(1/2)+arctan(1/3)
a=2 a²+1=5=1*5 donc arctan(1/2)=arctan(1/3)+arctan(1/7) et
arctan(1/2)=arctan(1/(2−1))+arctan(1/(5−2))=arctan(1)−arctan(1/3) a=3 a²+1=10=1*10=2*5 donc arctan(1/3)=arctan(1/4)+arctan(1/13)=arctan(1/5)+arctan(1/8) et arctan(1/3)=arctan(1/(3−2))+arctan(1/(3−5))=arctan(1)−arctan(1/2)
a=5 a²+1=26=1*26=2*13 donc arctan(1/5)=arctan(1/7)+arctan(1/18)=arctan(1/6)+arctan(1/31) et
arctan(1/5)=arctan(1/(5−2))+arctan(1/(5−13))=arctan(1/3)−arctan(1/8) a=7 a²+1=50=1*50=2*25 donc arctan(1/7)=arctan(1/8)+arctan(1/57)=arctan(1/9)+arctan(1/32) et
arctan(1/7)=arctan(1/(7−2))+arctan(1/(7−25))=arctan(1/5)−arctan(1/18)
On en déduit donc que :

et en utilisant π/4=4arctan(1/5)−arctan(1/239) on retrouve facilement les 2 formules utilisées par Jean Guilloud :
6arctan(1/8)+(6−4)arctan(1/57)+arctan(1/239)=
6(π/4−2arctan(1/3))−4(π/4−3arctan(1/3)+arctan(1/5))−π/4+4arctan(1/5)=
π/4
et
12arctan(1/18)+8arctan(1/57)−5arctan(1/239)=
12(−π/4+2arctan(1/3)+arctan(1/5))+8(π/4−3arctan(1/3)+arctan(1/5))+5(π/4−4arctan(1/5))=π/4

4.4  Pour bien afficher les décimales

//approx de arctan(a) par a-a^3/3+..+(-1)^n*a^(2n+1)/(2n+1) 
//valeur de la serie est calculee par la methode de Horner
gregory(a,n):={
local t,k;
t:=1/(2*n-1);
for (k:=2*n-3;k>0;k:=k-2) {
t:=1/k-a^2*t;
}
return (a*t);
}
:;
//pour afficher selon un tableau
affiche1(fl,n):={
local s,p,q,j,k,L,LT;
L:="";
s:=string(fl);
p:=size(s);
q:=iquo(p,11*n)
si irem(p,11*n)==0 alors 
LT:=[0$q];sinon 
LT:=[0$(q+1)];
fsi;
pour j de 0 jusque q-1 faire
pour k de 1 jusque 11 faire
L:=L+mid(s,0,n)+" ";
s:=mid(s,n);
fpour;
LT[j]=<[L];
L:="";
fpour;
si r!=0 alors
pour k de 1 jusque size(s) faire
L:=L+mid(s,0,n)+" ";
s:=mid(s,n);
fpour;
LT[q]=<[L];
fsi;
return LT;
}:;
//pour afficher selon une liste de chaines
affiche2(fl,n):={
local s,p,q,j,k,L,LT;
L:="";
s:=string(fl);
p:=size(s);
q:=iquo(p,11*n);
r:=irem(p,11*n);
si r==0 alors 
LT:=[0$q];sinon 
LT:=[0$q+1];
fsi;
pour j de 0 jusque q-1 faire
pour k de 1 jusque 11 faire
L:=L+mid(s,0,n)+" ";
s:=mid(s,n);
fpour;
LT[j]=<L;
L:="";
fpour;
si r!=0 alors
pour k de 1 jusque size(s) faire
L:=L+mid(s,0,n)+" ";
s:=mid(s,n);
fpour;
LT[q]=<L;
fsi;
return LT;
}:;

On vérifie pour n=2900 on a (1/5)⁸⁵/5801<0.33e−4058
et on vérifie pour n=860 on a (1/239)¹¹⁷²¹/1721=0.35e−4096.
On tape :
(1/5.)^5801/5801,(1/239.)^1721/1721
On obtient :
0.3247147235398264366269e-4058,0.347881711410409416161e-4096
Digits:=4000
On tape :
LR1:=affiche1(evalf(16*gregory(1/5,2900)-4*gregory(1/239,860)),5)
On obtient comme réponse une matrice.
Ou on tape :
LR2:=affiche2(evalf(16*gregory(1/5,2900)-4*gregory(1/239,860)),5)
pour j de 0 jusque size(LR2)-1 faire print(LR2[j]) fpour On obtient après 0.88s:

3.141 59265 35897 93238 46264 33832 79502 88419 71693 99375 10582 
09749 44592 30781 64062 86208 99862 80348 25342 11706 79821 48086 
51328 23066 47093 84460 95505 82231 72535 94081 28481 11745 02841 
02701 93852 11055 59644 62294 89549 30381 96442 88109 75665 93344 
61284 75648 23378 67831 65271 20190 91456 48566 92346 03486 10454 
32664 82133 93607 26024 91412 73724 58700 66063 15588 17488 15209 
20962 82925 40917 15364 36789 25903 60011 33053 05488 20466 52138 
41469 51941 51160 94330 57270 36575 95919 53092 18611 73819 32611 
79310 51185 48074 46237 99627 49567 35188 57527 24891 22793 81830 
11949 12983 36733 62440 65664 30860 21394 94639 52247 37190 70217 
98609 43702 77053 92171 76293 17675 23846 74818 46766 94051 32000 
56812 71452 63560 82778 57713 42757 78960 91736 37178 72146 84409 
01224 95343 01465 49585 37105 07922 79689 25892 35420 19956 11212 
90219 60864 03441 81598 13629 77477 13099 60518 70721 13499 99998 
37297 80499 51059 73173 28160 96318 59502 44594 55346 90830 26425 
22308 25334 46850 35261 93118 81710 10003 13783 87528 86587 53320 
83814 20617 17766 91473 03598 25349 04287 55468 73115 95628 63882 
35378 75937 51957 78185 77805 32171 22680 66130 01927 87661 11959 
09216 42019 89380 95257 20106 54858 63278 86593 61533 81827 96823 
03019 52035 30185 29689 95773 62259 94138 91249 72177 52834 79131 
51557 48572 42454 15069 59508 29533 11686 17278 55889 07509 83817 
54637 46493 93192 55060 40092 77016 71139 00984 88240 12858 36160 
35637 07660 10471 01819 42955 59619 89467 67837 44944 82553 79774 
72684 71040 47534 64620 80466 84259 06949 12933 13677 02898 91521 
04752 16205 69660 24058 03815 01935 11253 38243 00355 87640 24749 
64732 63914 19927 26042 69922 79678 23547 81636 00934 17216 41219 
92458 63150 30286 18297 45557 06749 83850 54945 88586 92699 56909 
27210 79750 93029 55321 16534 49872 02755 96023 64806 65499 11988 
18347 97753 56636 98074 26542 52786 25518 18417 57467 28909 77772 
79380 00816 47060 01614 52491 92173 21721 47723 50141 44197 35685 
48161 36115 73525 52133 47574 18494 68438 52332 39073 94143 33454 
77624 16862 51898 35694 85562 09921 92221 84272 55025 42568 87671 
79049 46016 53466 80498 86272 32791 78608 57843 83827 96797 66814 
54100 95388 37863 60950 68006 42251 25205 11739 29848 96084 12848 
86269 45604 24196 52850 22210 66118 63067 44278 62203 91949 45047 
12371 37869 60956 36437 19172 87467 76465 75739 62413 89086 58326 
45995 81339 04780 27590 09946 57640 78951 26946 83983 52595 70982 
58226 20522 48940 77267 19478 26848 26014 76990 90264 01363 94437 
45530 50682 03496 25245 17493 99651 43142 98091 90659 25093 72216 
96461 51570 98583 87410 59788 59597 72975 49893 01617 53928 46813 
82686 83868 94277 41559 91855 92524 59539 59431 04997 25246 80845 
98727 36446 95848 65383 67362 22626 09912 46080 51243 88439 04512 
44136 54976 27807 97715 69143 59977 00129 61608 94416 94868 55584 
84063 53422 07222 58284 88648 15845 60285 06016 84273 94522 67467 
67889 52521 38522 54995 46667 27823 98645 65961 16354 88623 05774 
56498 03559 36345 68174 32411 25150 76069 47945 10965 96094 02522 
88797 10893 14566 91368 67228 74894 05601 01503 30861 79286 80920 
87476 09178 24938 58900 97149 09675 98526 13655 49781 89312 97848 
21682 99894 87226 58804 85756 40142 70477 55513 23796 41451 52374 
62343 64542 85844 47952 65867 82105 11413 54735 73952 31134 27166 
10213 59695 36231 44295 24849 37187 11014 57654 03590 27993 44037 
42007 31057 85390 62198 38744 78084 78489 68332 14457 13868 75194 
35064 30218 45319 10484 81005 37061 46806 74919 27819 11979 39952 
06141 96634 28754 44064 37451 23718 19217 99983 91015 91956 18146 
75142 69123 97489 40907 18649 42319 61567 94520 80951 46550 22523 
16038 81930 14209 37621 37855 95663 89377 87083 03906 97920 77346 
72218 25625 99661 50142 15030 68038 44773 45492 02605 41466 59252 
01497 44285 07325 18666 00213 24340 88190 71048 63317 34649 65145 
39057 96268 56100 55081 06658 79699 81635 74736 38405 25714 59102 
89706 41401 10971 20628 04390 39759 51567 71577 00420 33786 99360 
07230 55876 31763 59421 87312 51471 20532 92819 18261 86125 86732 
15791 98414 84882 91644 70609 57527 06957 22091 75671 16722 91098 
16909 15280 17350 67127 48583 22287 18352 09353 96572 51210 83579 
15136 98820 91444 21006 75103 34671 10314 12671 11369 90865 85163 
98315 01970 16515 11685 17143 76576 18351 55650 88490 99898 59982 
38734 55283 31635 50764 79185 35893 22618 54896 32132 93308 98570 
64204 67525 90709 15481 41654 98594 61637 18027 09819 94309 92448 
89575 71282 89059 23233 26097 29971 20844 33573 26548 93823 91193 
25974 63667 30583 60414 28138 83032 03824 90375 89852 43744 17029 
13276 56180 93773 44403 07074 69211 20191 30203 30380 19762 11011 
00449 29321 51608 42444 85963 76698 38952 28684 78312 35526 58213 
14495 76857 26243 34418 93039 68642 62434 10773 22697 80280 73189 
15441 10104 46823 25271 62010 52652 27211 

4.5  La formule de Bailey-Borwein-Plouffe (BBP)

4.5.1  L’historique

En 1995, Plouffe (avec Bailey et Borwein) a découvert la formule de Bailey-Borwein-Plouffe (BBP) qui permet de calculer le n-ième bit de π sans avoir à calculer d’autres bits. Un an plus tard, il publie un nouvel article sur la formule, permettant de déterminer le n-ième chiffre en base 10 de π, mais le temps de calcul, bien que relativement court, n’est pas linéaire.

Il est également le co-auteur de l’Encyclopédie en ligne des suites de nombres entiers.

L’Inverseur de Plouffe est une page web qui contient plus de 200 millions de constantes mathématiques. Un répertoire est accessible et contient plus de 3,93 milliards de constantes à une précision de 64 chiffres décimaux au 21 juillet 2009.

Pour l’anecdote, Simon Plouffe a détenu en 1977 le record Guinness de mémorisation des décimales de π, avec 4 096 décimales. Il en avait mémorisé 4 400, mais en a récité seulement 4 096 parce que 4096=2¹²).

La formule de Plouffe ou formule BBP est

4.5.2  Démonstration de la formule de Plouffe en devoir (Université de Lorraine)

Le but de cet exercice est de montrer la formule de Plouffe

1. Montrer que pour tout a rèel avec 0 < a < 1, pour tous n, p ≥ 1 on a

   ∫

   a 0

   xn−1 1−xp

   =

   +∞ ∑ k=0

   akp+n kp+n

2. En déduire que

   +∞ ∑ k=0

   1 16k

   (

   4 8k+1

   −

   2 8k+4

   −

   1 8k+5

   −

   1 8k+6

   )=16

   ∫

   1 0

   4−2x3−x4−x5 16−x8

   dx

3. Simplifier la fraction 4−2x³−x⁴−x⁵/16−x⁸ (on pourra remarquer que le polynôme 4−2x³−x⁴−x⁵ a une racine simple), puis décomposer la fraction obtenue en éléments simples dans ℚ[X].
4. Conclure

Solution avec l’aide de Xcas On utilise ici Xcas en mode réel pour vérifier ou pour faire des calculs.

1. Pour |x| < 1 et p>1, on a |x^p | < 1. On a donc le développement en série entière :

   xn−1 1−xp

   =xn−1

   +∞ ∑ k=0

   xkp=

   +∞ ∑ k=0

   xkp+n−1

   Avec Xcas, on vérifie et on tape :
   assume(p,integer);sum(x^(k*p),k=0..inf)
   On obtient : x^(-p)/(x^(-p)-1)
   A l’intérieur du rayon de convergence, on peut intégrer terme à terme et comme ∫₀^ax^kp+n−1dx=a^kp+n/kp+n, on obtient :

   ∫

   a 0

   xn−1 1−xp

   dx=

   +∞ ∑ k=0

   akp+n kp+n

   Avec Xcas, on vérifie et on tape :
   int(x^(k*p+n-1),x=0..a)
   On obtient : a^(k*p+n)/(k*p+n)
2. Remarquons tout d’abord que :
   ∫₀¹x^m/16−x⁸dx=∫₀¹x^m/16(1−(x/√2)⁸dx= ∫₀^√2/21/16u^m*√2^m+1/1−u⁸du
   Avec Xcas, on tape pour faire le changement de variable x=u√2 dans l’intégrale :
   subst(’int(x^m/(16-x^8),x,0,1)’,x=sqrt(2)*u)
   On obtient :
   int((sqrt(2)*u)^m*1/(16-(sqrt(2)*u)^8)*sqrt(2),
            u,0,(sqrt(2))/2)
   Puis, on sélectionne sqrt(2)*u)^m*1/(16-(sqrt(2)*u)^8)*sqrt(2), on appelle simplify et on obtient :
   int(-(sqrt(2)*u)^m*sqrt(2)/(16*u^8-16), u,0,(sqrt(2))/2)
   D’aprés ce qui précéde (pour p=8, m=n−1 et a=√2/2) et puisque :
   √2^m+1*a^8k+m+1=1/√2^8*k=1/16^*k, on a :
   

   ∫

   1 0

   xm 16−x8

   dx=

   1 16

   +∞ ∑ k=0

   1 16k*(8k+m+1)

   =

   1 16

   +∞ ∑ k=0

   1 16k*(8k+m+1)

   Donc :
   

   16

   ∫

   1 0

   4−2x3−x4−x5 16−x8

   dx=

   +∞ ∑ k=0

   1 16k

   (

   4 8k+1

   −

   2 8k+4

   −

   2 8k+5

   −

   1 8k+6

   )

3. Factorisation de 4−2x³−x⁴−x⁵
   1 est racine évidente de 4−2x³−x⁴−x⁵.
   On tape :
   factor(4-2x^3-x^4-x^5)
   On obtient :
   -(x-1)*(x^2+2)*(x^2+2*x+2)
   On a donc la factorisation

   x5 + x4 + 2x3 −4 = (x−1)(x2 + 2)(x2 + 2x + 2)

   Factorisation de 16−x⁸
   On tape :
   factor(16-x^8)
   On obtient :
   -(x^2-2)*(x^2+2)*(x^2-2*x+2)*(x^2+2*x+2)
   Simplification de 4−2x³−x⁴−x⁵/16−x⁸
   On tape :
   normal((4-2x^3-x^4-x^5)/(16-x^8))
   On obtient :
   (x-1)/(x^4-2*x^3+4*x-4)
   Décomposition en éléments simples de 4−2x³−x⁴−x⁵/16−x⁸
   On tape :
   partfrac((x-1)/(x^4-2*x^3+4*x-4))
   On obtient :
   x/((x^2-2)*4)+(-x+2)/((x^2-2*x+2)*4)
   Ou bien , on tape directement :
   partfrac((4-2x^3-x^4-x^5)/(16-x^8))
   On obtient :
   x/((x^2-2)*4)+(-x+2)/((x^2-2*x+2)*4)
4. Conclusion
   On a donc :
   16∫₀¹4−2x³−x⁴−x⁵/16−x⁸dx=4∫₀¹x/x²−2dx+4∫₀¹−x+2/x²−2*x+2dx=
   2ln(|1²−2|)−2ln(|0²−2|)−2∫₀¹2x−2/x²−2*x+2dx+4∫₀¹1/x−1)²+1dx=−2ln(2)+2ln(2)+4*atan(0)−4*atan(−1)=π
   On tape :
   normal(16*int((4-2x^3-x^4-x^5)/(16-x^8),x=0..1))
   On obtient :
   pi
   Donc :

   π=

   +∞ ∑ n=0

   1 16k

   (

   4 8k+1

   −

   2 8k+4

   −

   1 8k+5

   −

   1 8k+6

   )

   Avec Xcas, en mode réel, on peut taper directement :
   simplify(sum(1/16^k*(4/(8*k+1)-2/(8*k+4)-1/(8*k+5)-
   1/(8*k+6)),k=0..inf))
   On obtient :
   pi

4.5.3  Le n-ième chiffre après la virgule de π en base 16 avec la formule BBP

La formule de Bailey-Borwein-Plouffe (BBP) permet de calculer le n-ième chiffre après la virgule de π en base 16 sans avoir à calculer d’autres chiffres.
On remarque que si n≥ 1, le n+1-ième chiffre après la virgule de π en base 16 est aussi le 1-ier chiffre après la virgule de 16ⁿπ.
On cherche donc 1-ier chiffre après la virgule de 16ⁿπ en base 16 et on a :

Exemple d’écriture en base 16 d’une approximation de π
Avec Xcas, on tape :
Digits:=25
a:=evalf(pi)
On obtient :
3.1415926535897932384626434
On calcule 16^ka pour k=0..21, on tape :
L:=floor(16^k*a)$(k=0..21)
Puis on convertit le résultat obtenu en base 16 :
LB:=revlist(convert(L[k],base,16))$(k=0..21))
On obtient :
[3],[3,2],[3,2,4],[3,2,4,3],[3,2,4,3,15],[3,2,4,3,15,6],
[3,2,4,3,15,6,10],[3,2,4,3,15,6,10,8],[3,2,4,3,15,6,10,8,8],
[3,2,4,3,15,6,10,8,8,8],[3,2,4,3,15,6,10,8,8,8,5]...
On tape :
L1:=LB[21];
On obtient :
[3,2,4,3,15,6,10,8,8,8,5,10,3,0,8,13,3,1,3,1,9,8]
Le k-ième chiffre après la virgule de π en base 16 est L1[k]
On vérifie, on tape :
evalf(sum(L1[k]/16^k,k=0..10),evalf(pi))
On obtient avec Digits:=21:
3.1415926535897932384626434,3.1415926535897932384626434
On utilise maintenant la formule BBP pour calculer le k-ième chiffre après la virgule de π en base 16.
Posons :

Le calcul des premiers chiffres de S_n(a) permettra d’obtenir ceux de 16ⁿπ, par la relation :
16ⁿπ = 4 S_n(1)−2 S_n(4)−S_n(5)−S_n(6)
Découpons la somme S_n(a) en deux sommes (l’une lorsque k varie entre 0 et n−1 et l’autre lorsque k varie entre n et +∞ :

Calculons les premiers chiffres de A_n(a) et B_n(a).

Les premiers chiffres après la virgule de 4A_n(1)−2A_n(4)−A_n(5)−A_n(6)
Il faut trouver la partie fractionnaire de 16^n−k/8k+a.
La partie entière de 16^n−k/8k+a est tres grande et on va utiliser l’arithmétique modulaire pour la soustraire à 16^n−k/8k+a.
On utilise, pour cela, la division euclidienne de 16^n−k par 8k+a :

c’est à dire 16^n−k=r mod(8k+a).
Donc 16^n−k/8k+a=q+r/8k+a.
Comme q∈ ℤ et r/8k+a<1, la partie fractionnaire de 16^n−k/8k+a est r/8k+a avec 16^n−k=r mod8k+a.
On veut avoir les premiers chiffres après la virgule de r/8k+a si 16^n−k=r mod8k+a.
Il faut donc calculer la partie fractionnaire de la suite SA (le nom de cette suite sur deux lettres fait référence à une variable d’une commande Xcas plus bas) définie par :

où :

Attention
Si frac(SA_n<0) alors la partie fractionnaire de :
4A_n(1)−2A_n(4)−A_n(5)−A_n(6) vaut 1+frac(SA_n) et
Si frac(SA_n>0) alors la partie fractionnaire de :
4A_n(1)−2A_n(4)−A_n(5)−A_n(6) vaut frac(SA_n)

Exemple
Si n=10 et a=1, pour obtenir le 11-ième chiffre apres la virgule de π en base 16, on doit calculer la partie fractionnaire de:
4Ã₁₀(1)−2Ã₁₀(4)−Ã₁₀(5)−Ã₁₀(6) avec
Ã₁₀(a)=∑_k=0⁹(16^10−kmod(8k+a))/8k+a.
On tape pour avoir 4Ã₁₀(1)−2Ã₁₀(4)−Ã₁₀(5)−Ã₁₀(6) :
SA:=evalf(sum(4*powmod(16,10-k,8k+1)/(8k+1)-
2*powmod(16,10-k,8k+4)/(8k+4)-powmod(16,10-k,8k+5)/(8k+5)-
powmod(16,10-k,8k+6)/(8*k+6),k=0..9))
On obtient -2.3654468582027340740994524
Donc la partie fractionnaire de 4A₁₀(1)−2A₁₀(4)−A₁₀(5)−A₁₀(6) est donc :
3+SA soit 3-2.36544685820273407=0.6345531417973

Les premiers chiffres après la virgule de 4B_n(1)−2B_n(4)−B_n(5)−B_n(6)
Posons : b_k=16^n−k/8k+a
Le premier terme de la somme B_n(a) est :
b_n=1/8n+a et donc b_n<1
Calculons b_k/b_k+1 :
b_k/b_k+1=16(8k+8+a)/8k+a=16(1+8/8k+a)>16 Donc pour obtenir B_n(a) avec une précision de p chiffres après la virgule, on va calculer les p+10 premiers termes de la somme pour tenir compte des retenues éventuelles .
Il suffit donc de calculer :
B_n(a)=∑_k=n^n+p+10 16^n−k/8k+a
Exemple (suite)
Si n=10, p=4 et a=1, on a :
B₁₀(1)=∑_k=10^10+4+10 16^10−k/8k+1
On tape pour avoir 4B₁₀(1)−2B₁₀(4)−B₁₀(5)−B₁₀(6) :
SB:= evalf(sum(4*16^(10-k)/(8k+1)-2*16^(10-k)/(8k+4)-
16^(10-k)/(8k+5)-16^(10-k)/(8k+6),k=10..24))
On obtient 0.23002595422921915411944196e-2

Conclusion Au final, pour obtenir les n premiers chiffres de π en base 16, il faut calculer les premiers chiffres après la virgule en base 16 de :
π_n=4S_n(1)−2S_n(4)−S_n(5)−S_N(6)
avec S_n(a)=∑_k=0ⁿ⁻¹ 16^n−k[8k+a]/8k+a+∑_k=n^n+p+10 16^n−k/8k+a
Fin de l’exemple
On tape :
3+SA+SB
On obtient 0.63685340133955811744174205
On tape pour avoir le 11-ième chiffre après la virgule en base 16 de π i.e. le premier chiffre après la virgule en base 16 de 3+SA+SB, on tape :
floor(16*(3+SA+SB)),L1[11]
On obtient 10,10
Puisque p=4, on tape pour avoir les quatre premiers chiffres après la virgule en base 16 de revlist(convert(floor(16^4*(3+SA+SB)),base,16))
On obtient [10,3,0,8]
On vérifie en tapant [L1[k]$(k=11..14)]
On obtient bien :
[10,3,0,8]
Un autre exemple pour avoir de la 13-ième à la 16-ième décimale en base 16 de π
On tape (n=12) :
SA2:=evalf(sum(4*powmod(16,12-k, 8k+1)/(8*k+1)-
2*powmod(16,12-k,8k+4)/(8k+4)- powmod(16,12-k,8k+5)/(8k+5)- powmod(16,12-k,8k+6)/(8*k+6),k=0..11))
On obtient :
-0.11967148118496058569703054e2
On tape (n=12 et p=4) :
SB2:=evalf(sum(4*16^(12-k)/(8k+1)- 2*16^(12-k)/(8k+4)-
16^(12-k)/(8k+5)- 16^(12-k)/(8k+6),k=12..26))
On obtient :
0.16188614229366348745743565e-2
On tape :
floor(16^k*(12+SA2+SB2))$(k=1..4)
On obtient :
0,8,141,2259
On tape :
revlist(convert(0,base,16)),revlist(convert(8,base,16)), revlist(convert(141,base,16)),revlist(convert(2259,base,16))
On obtient :
[8,13],[8,13,3]
Attention il ne faut pas oublier les zéros éventuels qui ne sont pas marqués !!! car on doit renvoyer une liste de longueur 1 pour la 12-ième... et une liste de longueur 4 pour la 15-ième décimale : [0],[0,8],[0,8,13],[0,8,13,3]
Donc la 12-ième décimale est 0 la 14-ième est 8 et les suivantes sont 13 et 3
On vérifie :
L1[k]$(k=13..16)
On obtient bien :
0,8,13,3 Exercice
Trouver la 10000-ième décimale de π en base 16.

On pose n=10⁴−1 et p=4 et on tape :
n:=10000-1
SA3:=evalf(sum((4*powmod(16,n-k, 8k+1)/(8*k+1)-
2*powmod(16,n-k, 8k+4)/(8k+4)- powmod(16,n-k,8k+5)/(8k+5)- powmod(16,n-k,8k+6)/(8*k+6)),k,0,(n-1)))
On obtient :
0.63408883081045966923274262e2
On tape :
SB3:=evalf(sum(4*16^(n-k)/(8k+1)-
2*16^(n-k)/(8k+4)-16^(n-k)/(8k+5)- 16^(n-k)/(8k+6),k=n..n+14))
On obtient :
0.25002812808619387546276337e-8
On tape :
revlist(convert(floor(16^4*(frac(SA3)+SB3)),base,16))
On obtient la 10000-ième, 10001-ième, 10002-ième,10003 ième décimale de π en base 16 :
[6,8,10,12]

4.5.4  Le programme Xcas

À l’aide de la formule BBP, on écrit un programme pour trouver, le n-ième,... le n+p−1-ième chiffre après la virgule, en base 16, de π.
On tape :

pichiffre16(n,p):={
local SA,SB,L,s,k;
si n<0 alors retourne "n doit etre >=0" fsi;
n:=n-1;
SA:=sum((4*powmod(16,n-k, 8k+1)/(8*k+1)-
     2*powmod(16,n-k, 8k+4)/(8k+4)-
     powmod(16,n-k,8k+5)/(8k+5)-
     powmod(16,n-k,8k+6)/(8*k+6)),k,0,(n-1));
SA:=frac(SA);
si SA<0 alors SA:=1+SA fsi;
SB:=sum(4*16^(n-k)/(8k+1)-2*16^(n-k)/(8k+4)
    -16^(n-k)/(8k+5)-16^(n-k)/(8k+6),k=n..n+p+10);
L:=revlist(convert(floor(16^p*frac(SA+SB)),base,16));
s:=size(L);
pour k de 0 jusque p-s-1 faire L:=prepend(L,0); fpour;
retourne L;
}:;

On tape :
pichiffre16(13,5)
On obtient : [0,8,13,3,1]
On tape pour avoir l’ècriture en base 16 d’une approximation de π:
pichiffre16(0,25)
On obtient : [3,2,4,3,15,6,10,8,8,8,5,10,3,0,8,13,3,1,3,1,9,8,10,2,14]] On tape :
pichiffre16(10000,4)
On obtient : [6,8,10,12]

4.5.5  Formule d’Adamchick et Wagon

Avec Xcas, on tape :
simplify(sum((-1)^k/4^k*(2/(4k+1)+2/(4k+2)+1/(4k+3)),k,0,inf))
On obtient : 1/2*pi+atan(1/2)+atan(2)
On tape :
simplify(1/2*pi+atan(1/2)+atan(2))
On obtient : pi

4.5.6  Formule de Plouffe généralisée

Notons S_n=∑_k=0^∞1/16^k(8k+n).
La formule de Plouffe généralisée est :

Avec Xcas, on tape :
simplify(sum(1/16^k*((4+8r)/(8k+1)-8r/(8k+2)-4r/(8k+3)-
(2+8r)/(8k+4)-(1+2r)/(8k+5)-(1+2r)/(8k+6)+r/(8k+7)),k,0,inf)
On obtient : pi
Démonstration
Pour r=0 on retrouve la formule BBP ;
Pour r=−1/4 on retrouve la formule d’Adamchick et Wagon sous une forme plus détaillée.
Posons :
α=1−i=√2e^−iπ/4 et calculons de deux façons l’intégrale :

. Elle est reliée aux S_n car :

Elle calculable par des méthodes élémentaires (en calculant séparément sa partie réelle et sa partie imaginaire ou avec un logarithme complexe), ou avec Xcas en mode complexe :

Ou on tape en mode complexe :
normal(int(1/(1-i-y),y,0,1))
On obtient : i/4*pi+1/2*ln(2)
L’égalité entre ces deux expressions de I équivaut à :

(1)    π=4Im(I)=2S₁+2S₂+S₃−1/2S₅−1/2S₆−1/4S₇,

(2)    ln(2)=2Re(I)= S₁−1/2S₃−1/2S₄−1/4S₅+1/8S₇+1/8S₈

Mais ln(2) peut s’exprimer en fonction des S_n :

(3) ∫₀¹2y/2−y²dy=[−ln(2−y²)]₀¹=ln(2) ∫₀¹y/1−y²/2dy=∑_{k≥ 0}1/(2k+2)2^k
      =S₂+1/2 S₄+1/4 S₆+1/8 S₈

On soustrait (2) et (3) et on a une relation entre les S_n :

(4)     0=2S₁−2S₂−S₃−2S₄−1/2S₅−1/2 S₆+1/4S₇
En multipliant (4) par 1+4r et en ajoutant ce produit à (1), on obtient l’égalité (0).