Algorithme de Chudnovsky

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)}}}$ peut être simplifié en ${\displaystyle {\frac {-(6n-1)(2n-1)(6n-5)}{10939058860032000n^{3}}}}$, donc$${\displaystyle f(n)=f(n-1)\cdot {\frac {-(6n-1)(2n-1)(6n-5)}{10939058860032000n^{3}}}}$$${\displaystyle f(0)=1}$ par définition de ${\displaystyle f}$, donc$${\displaystyle f(n)=\prod _{j=1}^{n}{\frac {-(6j-1)(2j-1)(6j-5)}{10939058860032000j^{3}}}}$$Cette définition de ${\displaystyle f}$ n'est pas défini pour ${\displaystyle n=0}$, on calcule donc le premier terme de la somme et l'ajoute à la nouvelle définition de ${\displaystyle f}$ en partant de ${\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 )}}$$Soit ${\displaystyle P(a,b)=\prod _{j=a}^{b-1}{-(6j-1)(2j-1)(6j-5)}}$ et ${\displaystyle Q(a,b)=\prod _{j=a}^{b-1}{10939058860032000j^{3}}}$, donc$${\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 )}}$$Soit ${\displaystyle S(a,b)=\sum _{k=a}^{b-1}{{\frac {P(a,k+1)}{Q(a,k+1)}}\cdot (545140134k+13591409)}}$ et ${\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 )}$ est une limite qui ne peut être qu'approché, on calcule à la place ${\displaystyle S(1,n)}$ et lorsque ${\displaystyle n}$ se rapproches de ${\displaystyle \infty }$, l'approximation de ${\displaystyle \pi }$ l'approximation s'améliore.$${\displaystyle \pi \approx {\frac {426880{\sqrt {10005}}}{13591409+S(1,n)}}}$$Par définition originale de ${\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)}}}$$

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...