алгоритм Чудновского

Алгоритм Чудновского собой быстрый метод вычисления цифр числа π , основанный на Рамануджана π представляет формулах . Издано братьями Чудновскими в 1988 году. его использовали для вычисления числа π с точностью до миллиарда десятичных знаков.

Он использовался при расчете мирового рекорда 2,7 триллионов цифр числа π в декабре 2009 года. 10 триллионов цифр в октябре 2011 года, 22,4 триллиона цифр в ноябре 2016 года, 31,4 триллиона цифр в сентябре 2018 г. – январе 2019 г., 50 триллионов цифр 29 января 2020 года, 62,8 триллиона цифр на 14 августа 2021 года, 100 триллионов цифр 21 марта 2022 года, и 105 триллионов цифр 14 марта 2024 года.

Алгоритм

Алгоритм основан на отрицательном числе Хегнера. $d=-163$, j -функция $j\left({\tfrac {1+i{\sqrt {163}}}{2}}\right)=-640320^{3}$и на следующих быстро сходящихся обобщенных гипергеометрических рядах :

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

Подробное доказательство этой формулы можно найти здесь:

Это тождество похоже на некоторые формулы Рамануджана , включающие π , и является примером серии Рамануджана-Сато .

Временная сложность алгоритма составляет $O\left(n(\log n)^{3}\right)$.

Оптимизации

Техника оптимизации, используемая для вычислений мировых рекордов, называется бинарным расщеплением .

Бинарное расщепление

Фактор ${\textstyle 1/{640320^{3/2}}}$ можно вынести из суммы и упростить до

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

Позволять $f(n)={\frac {(-1)^{n}(6n)!}{(3n)!(n!)^{3}(640320)^{3n}}}$и подставим это в сумму.

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

${\frac {f(n)}{f(n-1)}}$ можно упростить до ${\frac {-(6n-1)(2n-1)(6n-5)}{10939058860032000n^{3}}}$, так

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

$f(0)=1$ из первоначального определения $f$, так

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

Это определение $f$ не определено для $n=0$, поэтому вычислите первый член суммы и используйте новое определение $f$

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

Позволять $P(a,b)=\prod _{j=a}^{b-1}{-(6j-1)(2j-1)(6j-5)}$ и $Q(a,b)=\prod _{j=a}^{b-1}{10939058860032000j^{3}}$, так

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

Позволять $S(a,b)=\sum _{k=a}^{b-1}{{\frac {P(a,k+1)}{Q(a,k+1)}}\cdot (545140134k+13591409)}$ и $R(a,b)=Q(a,b)\cdot S(a,b)$

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

$S(1,\infty )$ никогда не может быть вычислено, поэтому вместо этого вычисляйте $S(1,n)$ и как $n$ подходы $\infty$, $\pi$ аппроксимация улучшится.

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

Из первоначального определения $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)}}$$

Рекурсивное вычисление функций

Рассмотрим значение $m$ такой, что $a<m<b$

• $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)$

Базовый вариант рекурсии

Учитывать $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)$

Код 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(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...

Примечания

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

См. также

• Серия Рамануджана – Сато
• Формула Бейли – Борвейна – Плуффа
• Алгоритм Борвейна
• Приближения π