Des milliers de décimales de π

Publié le 17/02/17

Le but de cet article est de proposer un algorithme conduisant à une fonction Python de dix lignes calculant rapidement des milliers de décimales du nombre \pi.

1. Une base de numération variable

On admet que : {\pi=2\displaystyle\sum_{k=0}^{+\infty}u_k}, où {u_k=\dfrac{(k!)^2\,2^{k}}{(2k+1)!}}.

On a ici {u_0=1}, et {u_k=\dfrac{k}{2k+1}\,u_{k-1}=\dfrac{1\,\cdot\,2\,\cdot\,3\,\cdots\,k}{3\,\cdot\,5\,\cdot\,7\,\cdots\,(2k+1)}} pour {k\ge1}.

Plus généralement, on considère des développements (supposés convergents) :
\begin{array}{rl}x&=\displaystyle\sum_{k=0}^{+\infty}p_k\,u_k=p_0+\displaystyle\sum_{k=1}^{+\infty}p_k\,\dfrac{1\,\cdot\,2\,\cdot\,3\,\cdots\,k}{3\,\cdot\,5\,\cdot\,7\,\cdots\,(2k+1)}\\\\&=p_0+\dfrac{1}{3}\left(p_1+\dfrac{2}{5}\left(p_2+\dfrac{3}{7}\left(p_3+\dfrac{4}{9}\left(p_4+\cdots\right)\right)\right)\right)\end{array}On notera {x=(p_0,p_1,p_2,\ldots)_{\mathcal{P}}}, donc {\pi=(2,2,2,\ldots)_{\mathcal{P}}}.

On dira que le développement {x=(p_0,p_1,p_2,\ldots)_{\mathcal{P}}} est régulier si :

  • on a {p_k\in\{0,\ldots,2k\}} pour tout {k\ge1},
  • pour tout entier {n\ge1}, il existe {m\ge n} tel que {p_m\lt 2k}.

2. Convergence d’un développement régulier

On se donne ici une développement régulier {x=(p_0,p_1,p_2,\ldots)_{\mathcal{P}}}.

En particulier, il existe {n_0\ge n} tel que {p_{n_0}\lt 2n_0}.

On va voir qu’un tel développement converge dans [0,2[.

Pour tout k\ge1, on a {u_k=\dfrac{k}{2k+1}\,u_{k-1}} donc {ku_k=ku_{k-1}-(k+1)u_k}.

Pour tous entiers {1\le n\le m}, on peut alors écrire :
{\begin{array}{rl}0\le \displaystyle\sum_{k=n}^{m}p_k u_k&\le 2\displaystyle\sum_{k=n}^{m}k u_k\\\\&= 2\displaystyle\sum_{k=n}^{m}(ku_{k-1}-(k+1)u_k)\\\\&=2(nu_{n-1}-(m+1)u_m)\end{array}}

Pour tout {m\ge1}, on a :
{u_m=\dfrac{1\,\cdot\,2\,\cdot\,3\,\cdots\,\cdot\,m}{3\,\cdot\,5\,\cdot\,7\,\cdots\,\cdot\,(2m+1)}\le 2^{-m}}On en déduit {\displaystyle\lim_{m\rightarrow+\infty}(m+1)u_m=0}.
Pour tout {n\ge1}, il en découle : {0\le \displaystyle\sum_{k=n}^{+\infty}p_k\,u_k\le 2n\,u_{n-1}}.

La dernière inégalité est stricte car {p_k\lt 2k} pour au moins un entier {k\ge n}.

Pour tout {n\ge1}, on a obtenu :
{\begin{array}{rl}0\le \displaystyle\sum_{k=n}^{+\infty}p_k\,u_k\lt 2n\,\dfrac{1\,\cdot\,2\,\cdot\,3\,\cdots\,\cdot\,(n-1)}{3\,\cdot\,5\,\cdot\,7\,\cdots\,\cdot\,(2n-1)}=2^{n+1}\dfrac{(n!)^2}{(2n)!}\end{array}}En particulier : {0\le \displaystyle\sum_{k=1}^{+\infty}p_k\,u_k\lt 2}.

Mais remarquons qu’il n’y a donc pas unicité de la représentation d’un réel dans le système {\mathcal{P}}.

Les calculs précédents donnent en effet {\displaystyle\sum_{k=n+1}^{+\infty}k\,u_k=(n+1)u_n}.

Ainsi (en notant 0_k pour le placement de 0 en position k) :
{\begin{array}{l}(0_0,0_1,\ldots,0_{n-1},0_n,n+1,n+2,\ldots)_{\mathcal{P}}\\\\=(0_0,0_1,\ldots,0_{n-1},n+1,0,0,\ldots)_{\mathcal{P}}\end{array}}Par exemple : {1=(1,0,0,\ldots)_{\mathcal{P}}=(0,1,2,3,\ldots)_{\mathcal{P}}}

On encore : {\dfrac23=(0,2,0,0,\ldots)_{\mathcal{P}}=(0,0,2,3,\ldots)_{\mathcal{P}}}.

La partie entière de {(p_0,p_1,p_2,\ldots)} est égale à {p_0} ou à {p_0+1}.

Elle est égale à {1} si {(0,p_1,p_2,p_3\ldots)\ge (0,1,2,3,\ldots)} (ordre lexicographique).

3. Algorithme de réduction à un développement régulier

Soit {x_n=(p_0,p_1,p_2,\ldots,p_{n-1},0,\ldots)_{\mathcal{P}}} un développement régulier fini.

Pour {m\in\mathbb{N}^*}, soit {y=10^m x_n=(10^m p_0,10^m p_1,\ldots,10^mp_{n-1},0,\ldots)_{\mathcal{P}}}.

Ce développement n’est certainement pas régulier.

On va donc proposer un algorithme convertissant {y} en un développement régulier.

Avec les notations de l’énoncé :
{x_{n}=(p_0,p_1,p_2,\ldots,p_{n-1},0,\ldots)_{\mathcal{P}}=\displaystyle\sum_{k=0}^{n-1}p_ku_k}Ainsi {10^mx_n=\displaystyle\sum_{k=0}^{n-1}(10^mp_k)u_k}, à convertir à une forme régulière.

On va se servir de {(2k+1)u_k=k\,u_{k-1}} (valable pour tout {k\ge1}).

  • On effectue la division de {10^m\,p_{n-1}} par {2n-1}.

    Elle s’écrit : {10^m\,p_{n-1}=q_{n-1}(2n-1)+r_{n-1}}, où {0\le r_{n-1}\le 2(n-1)}.
    {\begin{array}{rl}10^mx_{n}&=\displaystyle\sum_{k=0}^{n-3}10^m\,p_k\,u_k+10^m\,p_{n-2}\,u_{n-2}+(q_{n-1}(2n-1)+r_{n-1})\,u_{n-1}\\\\&=\displaystyle\sum_{k=0}^{n-3}10^m\,p_k\,u_k+(10^m\,p_{n-2}+(n-1)q_{n-1})\,u_{n-2}+r_{n-1}\,u_{n-1}\end{array}}La division s’est donc traduite par une retenue {(n-1)q_{n-1}}, reportée sur {10^m\,p_{n-2}}.

  • On divise maintenant {p'_{n-2}=10^m\,p_{n-2}+(n-1)q_{n-1}} par {2n-3}.

    Cela s’écrit {p'_{n-2}=q_{n-2}(2n-3)+r_{n-2}}, où {0\le r_{n-2}\le2(n-2)}:
    {\begin{array}{rl}10^mx_n&=\displaystyle\sum_{k=0}^{n-4}10^m\,p_k\,u_k+(10^m\,p_{n-3}+(n-2)q_{n-2})u_{n-3}\\\\&+r_{n-2}\,u_{n-2}+r_{n-1}\,u_{n-1}\end{array}}

  • On continue ces divisions successives avec retenues.

    A l’avant-dernière étape, on obtient :
    {10^mx_n=10^m\,p_{0}+(10^m\,p_{1}+2q_{2})\,u_1+\displaystyle\sum_{k=2}^{n-1} r_k\,u_k}Il reste à diviser {p'_1=10^m\,p_{1}+2q_{2}} par {3}, en utilisant {3u_1=u_0=1}

    Cette division s’écrit {p'_1=3q_1+r_1}, avec {r_1\in\{0,1,2\}}.

    Finalement, en posant {r_0=10^m\,p_{0}+q_1} :
    {10^mx_n=r_0+r_1\,u_1+\displaystyle\sum_{k=2}^{n-1} r_k\,u_k=\displaystyle\sum_{k=0}^{n-1} r_k\,u_k}

La conversion est terminée, car le développement obtenu est régulier.

Mais dans ce système, on a {[10^m\,x_n]=r_0} ou {[10^m\,x_n]=r_0+1}.

4. Jusqu’où aller pour une bonne approximation de \pi

On sait que {\pi=2\displaystyle\sum_{k=0}^{+\infty}u_k}. Pour {m\ge1}, on pose {\pi_m=2\displaystyle\sum_{k=0}^{m-1}u_k}.

Pour tout {k\ge1}, on a {u_k=\dfrac{k}{2k+1}\,u_{k-1}} donc {u_k\lt \dfrac12 \,u_{k-1}}.

On en déduit, pour tout {m\ge1} et tout {k\gt m} : {u_k\lt\Bigl(\dfrac12\Bigr)^{k-m}u_{m}}.

Ainsi: {\begin{array}{rl}0&\lt \pi-\pi_{m}\lt 2u_m\displaystyle\sum_{k=m}^{+\infty}\Bigl(\dfrac12\Bigr)^{k-m}=4u_m\end{array}}

D’autre part :
{\begin{array}{rl}u_m&=\dfrac{1\,\cdot\,2\,\cdot\,3\cdot\,\ldots\,\cdot m}{3\,\cdot\,5\,\cdot\,7\,\ldots\,\cdot(2m+1)}\\\\&\lt \dfrac{2}{3}\,\dfrac{1\,\cdot\,2\,\cdot\,3\,\cdot\,\ldots\,\cdot m}{2\,\cdot\,4\,\cdot\,6\,\cdot\,\ldots\,\cdot(2m)}=\dfrac{2}{3}\,2^{-m}\end{array}}Montrons maintenant que si {m\ge \dfrac{10}{3}\,n}, alors {0\lt \pi-\pi_m\lt 5\cdot10^{-n}}.

Pour cela il suffit de vérifier que {m\ge \dfrac{10}{3}\,n\Rightarrow 4u_m\lt 5\cdot10^{-n}}.

Mais cela découlera de {m\ge \dfrac{10}{3}\,n\Rightarrow\dfrac{8}{3}\cdot 2^{-m}\lt 5\cdot10^{-n}}

Il reste donc à montrer {8\cdot 2^{-m}\lt 15\cdot 10^{-n}}.

Autrement dire il doit vérifier {m\ln 2>n\ln 10+\ln 8-\ln 15}.

Or on a l’approximation : {\dfrac{\ln 10}{\ln 2}\approx 3.32\lt \dfrac{10}{3}}.

Ainsi : {m\ln 2>\dfrac{10}{3}\,n\ln 2>n\ln 10>n\ln 10+\ln 8-\ln 15}.

On a donc établi : {m\ge \dfrac{10}{3}\,n\Rightarrow 0\lt \pi-\pi_m\lt 5\cdot10^{-n}}, pour {n\ge1}.

5. Mise en oeuvre de l’algorithme

On rappelle la notation : {\pi_m=2\displaystyle\sum_{k=0}^{m-1}u_k}.

Dans la base (\mathbb{P}), \pi_m=(2,2,\ldots,2) de longueur {m}).

On sait maintenant qu’il suffit de choisir {m=[10n/3]} pour que {\pi_m} soit une valeur approchée par défaut de {\pi} à {5\cdot10^{-n}} près. Les deux nombres {\pi} et {\pi_{m}} ont alors les mêmes {m} premiers chiffres (sauf erreur d’arrondi éventuelle sur le tout dernier). Comme nous savons comment former les décimales successives de {\pi_m}, il nous est maintenant possible de trouver celles de {\pi}.

Voici une fonction chiffres_de_pi, d’argument entier n, renvoyant les n premiers chiffres de \pi (donc les n-1 premières décimales) :

Voici le calcul des {1000} premiers chiffres de {\pi} (1/100ième de seconde ur mon Mac mini!) :

‘3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198’

Les temps de calculs sont assez bluffants :