Algorithme de Chudnovski

Pour les articles homonymes, voir 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, ^[1] 31 400 milliards de chiffres en septembre 2018. –Janvier 2019, ^[2] 50 000 milliards de chiffres le 29 janvier 2020, ^[3] 62 800 milliards de chiffres le 14 août 2021^[4], 100 000 milliards de chiffres le 21 mars 2022^[5], et 105 000 milliards de chiffres en mars 14, 2024^[6].

Algorithme

L'algorithme est basé sur l'opposé du nombre de Heegner ${\displaystyle d=-163}$, la fonction j ${\displaystyle 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:

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

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 ${\displaystyle O\left(n(\log n)^{3}\right)}$^[7].

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

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

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

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

${\displaystyle {\frac {f(n)}{f(n-1)}}}$

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

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

${\displaystyle f(0)=1}$

${\displaystyle f}$

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

${\displaystyle f}$

${\displaystyle n=0}$

${\displaystyle f}$

${\displaystyle n=1}$

$${\displaystyle {\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 )}}$$

${\displaystyle P(a,b)=\prod _{j=a}^{b-1}{-(6j-1)(2j-1)(6j-5)}}$

${\displaystyle Q(a,b)=\prod _{j=a}^{b-1}{10939058860032000j^{3}}}$

$${\displaystyle {\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 )}}$$

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

${\displaystyle R(a,b)=Q(a,b)\cdot S(a,b)}$

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

${\displaystyle S(1,\infty )}$

${\displaystyle S(1,n)}$

${\displaystyle n}$

${\displaystyle \infty }$

${\displaystyle \pi }$

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

${\displaystyle R}$

${\displaystyle S(a,b)={\frac {R(a,b)}{Q(a,b)}}}$

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

Calcul récursif des fonctions

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

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

Construction de la récursion

Si ${\displaystyle b=a+1}$

• ${\displaystyle P(a,a+1)=-(6a-1)(2a-1)(6a-5)}$
• ${\displaystyle Q(a,a+1)=10939058860032000a^{3}}$
• ${\displaystyle S(a,a+1)={\frac {P(a,a+1)}{Q(a,a+1)}}\cdot (545140134a+13591409)}$
• ${\displaystyle 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

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

${\displaystyle 640320^{3}/24=10939058860032000}$

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

${\displaystyle 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).

• Portail de l'analyse