Algorithme de Chudnovski

L' algorithme de Chudnovsky est une méthode rapide pour calculer les chiffres de π, basée sur les formules de π de Ramanujan. Il a été publié par les frères Chudnovsky en 1988.

Il a été utilisé pour le calcul du record mondial de 2 700 milliards de chiffres de π en décembre 2009, 10 000 milliards de chiffres en octobre 2011, 22 400 milliards de chiffres en novembre 2016, 31 400 milliards de chiffres en septembre 2018. – janvier 2019, 50 000 milliards de chiffres le 29 janvier 2020, 62 800 milliards de chiffres le 14 août 2021, 100 000 milliards de chiffres le 21 mars 2022, et 105 000 milliards de chiffres en mars 14, 2024.

Algorithme

L'algorithme est basé sur l'opposé du nombre de Heegner $d=-163$, la fonction j $j\left({\tfrac {1+i{\sqrt {163}}}{2}}\right)=-640320^{3}$, et sur la série hypergéométrique généralisée à convergence rapide ci-dessous:

$${\frac {1}{\pi }}=12\sum _{k=0}^{\infty }{\frac {(-1)^{k}(6k)!(545140134k+13591409)}{(3k)!(k!)^{3}(640320)^{3k+3/2}}}$$

Une preuve détaillée de cette formule peut être trouvée ici :

Cette identité est similaire à certaines formules de Ramanujan impliquant π, et est un exemple de série Ramanujan-Sato.

La complexité temporelle de l'algorithme est en $O\left(n(\log n)^{3}\right)$.

Optimisations

La technique d'optimisation utilisée pour les calculs du record du monde est appelée division binaire.

Fractionnement binaire

Un facteur de ${\textstyle 1/{640320^{3/2}}}$ peut être retiré de la somme et simplifié en

$${\frac {1}{\pi }}={\frac {1}{426880{\sqrt {10005}}}}\sum _{k=0}^{\infty }{\frac {(-1)^{k}(6k)!(545140134k+13591409)}{(3k)!(k!)^{3}(640320)^{3k}}}$$

Soit $f(n)={\frac {(-1)^{n}(6n)!}{(3n)!(n!)^{3}(640320)^{3n}}}$, et en substituant dans la somme.

$${\frac {1}{\pi }}={\frac {1}{426880{\sqrt {10005}}}}\sum _{k=0}^{\infty }{f(k)\cdot (545140134k+13591409)}$$

${\frac {f(n)}{f(n-1)}}$ peut être simplifié en ${\frac {-(6n-1)(2n-1)(6n-5)}{10939058860032000n^{3}}}$, donc

$$f(n)=f(n-1)\cdot {\frac {-(6n-1)(2n-1)(6n-5)}{10939058860032000n^{3}}}$$

$f(0)=1$ par définition de $f$, donc

$$f(n)=\prod _{j=1}^{n}{\frac {-(6j-1)(2j-1)(6j-5)}{10939058860032000j^{3}}}$$

Cette définition de $f$ n'est pas défini pour $n=0$, on calcule donc le premier terme de la somme et l'ajoute à la nouvelle définition de $f$ en partant de $n=1$

$${\frac {1}{\pi }}={\frac {1}{426880{\sqrt {10005}}}}{\Bigg (}13591409+\sum _{k=1}^{\infty }{{\Bigg (}\prod _{j=1}^{k}{\frac {-(6j-1)(2j-1)(6j-5)}{10939058860032000j^{3}}}{\Bigg )}\cdot (545140134k+13591409)}{\Bigg )}$$

Soit $P(a,b)=\prod _{j=a}^{b-1}{-(6j-1)(2j-1)(6j-5)}$ et $Q(a,b)=\prod _{j=a}^{b-1}{10939058860032000j^{3}}$, donc

$${\frac {1}{\pi }}={\frac {1}{426880{\sqrt {10005}}}}{\Bigg (}13591409+\sum _{k=1}^{\infty }{{\frac {P(1,k+1)}{Q(1,k+1)}}\cdot (545140134k+13591409)}{\Bigg )}$$

Laisser $S(a,b)=\sum _{k=a}^{b-1}{{\frac {P(a,k+1)}{Q(a,k+1)}}\cdot (545140134k+13591409)}$ et $R(a,b)=Q(a,b)\cdot S(a,b)$

$$\pi ={\frac {426880{\sqrt {10005}}}{13591409+S(1,\infty )}}$$

$S(1,\infty )$ est une limite qui ne peut être qu'approché, on calcule à la place $S(1,n)$ et lorsque $n$ se rapproches de $\infty$, l'approximation de $\pi$ l'approximation s'améliore.

$$\pi \approx {\frac {426880{\sqrt {10005}}}{13591409+S(1,n)}}$$

Par définition originale de $R$, $S(a,b)={\frac {R(a,b)}{Q(a,b)}}$

$$\pi \approx {\frac {426880{\sqrt {10005}}\cdot Q(1,n)}{13591409Q(1,n)+R(1,n)}}$$

Calcul récursif des fonctions

• $P(a,b)=P(a,m)\cdot P(m,b)$

• $Q(a,b)=Q(a,m)\cdot Q(m,b)$
• $S(a,b)=S(a,m)+{\frac {P(a,m)}{Q(a,m)}}S(m,b)$
• $R(a,b)=Q(m,b)R(a,m)+P(a,m)R(m,b)$

Construction de la récursion

Si $b=a+1$

• $P(a,a+1)=-(6a-1)(2a-1)(6a-5)$
• $Q(a,a+1)=10939058860032000a^{3}$
• $S(a,a+1)={\frac {P(a,a+1)}{Q(a,a+1)}}\cdot (545140134a+13591409)$
• $R(a,a+1)=P(a,a+1)\cdot (545140134a+13591409)$

Code Python

import decimal

def binary_split(a, b):
    if b == a + 1:
        Pab = -(6*a - 5)*(2*a - 1)*(6*a - 1)
        Qab = 10939058860032000 * a**3
        Rab = Pab * (545140134*a + 13591409)
    else:
        m = (a + b) // 2
        Pam, Qam, Ram = binary_split(a, m)
        Pmb, Qmb, Rmb = binary_split(m, b)
        
        Pab = Pam * Pmb
        Qab = Qam * Qmb
        Rab = Qmb * Ram + Pam * Rmb
    return Pab, Qab, Rab

def chudnovsky(n):
    """Chudnovsky algorithm."""
    P1n, Q1n, R1n = binary_split(1, n)
    return (426880 * decimal.Decimal(10005).sqrt() * Q1n) / (13591409*Q1n + R1n)

print(chudnovsky(2))  # 3.141592653589793238462643384

decimal.getcontext().prec = 100
for n in range(2,10):
    print(f"{n=} {chudnovsky(n)}")  # 3.14159265358979323846264338...

Remarques

$e^{\pi {\sqrt {163}}}\approx 640320^{3}+743.99999999999925\dots$

$640320^{3}/24=10939058860032000$

$545140134=163\cdot 127\cdot 19\cdot 11\cdot 7\cdot 3^{2}\cdot 2$

$13591409=13\cdot 1045493$

Voir aussi

• Série Ramanujan-Sato
• Formule Bailey – Borwein – Plouffe
• L'algorithme de Borwein
• Approximations de π

Références

• (en) Cet article est partiellement ou en totalité issu de l’article de Wikipédia en anglais intitulé « Chudnovsky algorithm » (voir la liste des auteurs).