אלגוריתם צ'ודנובסקי

אלגוריתם צ'ודנובסקי הוא אלגוריתם מהיר לחישוב הספרות של פאי (מספר). האלגוריתם פורסם בשנת 1989 על ידי האחים דויד וגריגורי וולפוביץ' צ'ודנובסקי. האלגוריתם שימש לשבירת השיא העולמי לחישוב הספרות של פאי מספר פעמים: בשיאים שחישבו את 2.7 טריליון (כלומר 12^10) הספרות של פאי ב-2009, את 10 טריליון הספרות ב-2011, בנובמבר 2016 חישבו בעזרת האלגוריתם 22.4 טריליון ספרות של פאי ובספטמבר 2018 עד ינואר 2019 חושבו באמצעותו 31.4 טריליון ספרות].

סיבוכיות האלגוריתם כדי לחשב $n$ ספרות היא $O(n\log ^{3}n)$.

האלגוריתם מבוסס על הטור האינסופי:

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

אשר ניתן לפשטו לצורה ${\frac {(640320)^{3/2}}{12\pi }}={\frac {426880{\sqrt {10005}}}{\pi }}=\sum _{k=0}^{\infty }{\frac {(6k)!(545140134k+13591409)}{(3k)!(k!)^{3}\left(-262537412640768000\right)^{k}}}\!$ יישנם שלושה ביטויים מסוגים שונים, שהם סדרות בקצבים שונים וקבוע: בקצב פולינומי M_k, בקצב ליניארי L_k, בקצב אקספוננציאלי X_k, והקבוע C, שעל בסיסם ניתנת הנוסחה שהיא:

$\pi =C\left(\sum _{k=0}^{\infty }{\frac {M_{k}\cdot L_{k}}{X_{k}}}\right)^{-1}$,

כאשר: $C=426880{\sqrt {10005}}$
M_k+1 = M_k * (12k+2) * (12k+6) * (12k+10) / (k+1)^3, M₀ = 1 [M_k = (6k)! / ((3k)! * (k!)^3)]
L_k+1 = L_k + 545,140,134, L₀ = 13,591,409 [L_k = 13591409 + 545140134*k]
X_k+1 = X_k * -262,537,412,640,768,000, X₀ = 1 [X_k = (-640320)^3k = (-262537412640768000)^k]

ניתן אף לשפר את הנוסחה על ידי הצבה של: K_k+1 = K_k + 12, K₀ = 6
M_k+1 = M_k * (K_k^3 - 16K_k) / (k+1)^3 ,M₀ = 1

דוגמת קוד בפייתון של האלגוריתם

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(f"1 = {chudnovsky(2)}")  # 3.141592653589793238462643384

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