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 Modèle:Nombre de chiffres de π en décembre 2009, Modèle:Nombre de chiffres en octobre 2011, Modèle:Nombre de chiffres en novembre 2016, Modèle:Nombre de chiffres en septembre 2018. – janvier 2019, Modèle:Nombre de chiffres le 29 janvier 2020, Modèle:Nombre de chiffres le 14 août 2021, Modèle:Nombre de chiffres le 21 mars 2022, et Modèle:Nombre de chiffres le 14 mars 2024.

Algorithme

L'algorithme est basé sur l'opposé du nombre de Heegner $d=-163$, la fonction j $j\left(\frac{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}^{\mathrm{\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(n(\log n{)}^3)$.

Optimisations

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

Fractionnement binaire

Un facteur de $1/640320^{3/2}$ peut être retiré de la somme et simplifié en$$\frac{1}{\pi }=\frac{1}{426880\sqrt{10005}}\sum _{k=0}^{\mathrm{\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}^{\mathrm{\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}}(13591409+\sum _{k=1}^{\mathrm{\infty }}\left(\prod _{j=1}^k\frac{-(6j-1)(2j-1)(6j-5)}{10939058860032000j^3}\right)\cdot (545140134k+13591409))$$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}}(13591409+\sum _{k=1}^{\mathrm{\infty }}\frac{P(1,k+1)}{Q(1,k+1)}\cdot (545140134k+13591409))$$Soit $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,\mathrm{\infty })}$$$S(1,\mathrm{\infty })$ est une limite qui ne peut être qu'approché, on calcule à la place $S(1,n)$ et lorsque $n$ se rapproches de $\mathrm{\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 π