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

Алгоритм Чудновского собой быстрый метод вычисления цифр числа π , основанный на Рамануджана π представляет формулах . Издано братьями Чудновскими в 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 года, и 202 триллиона цифр 28 июня 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 =\lim _{n\to \infty }{\frac {426880{\sqrt {10005}}}{13591409+S(1,n)}}$$

Из первоначального определения $R$, $S(a,b)={\frac {R(a,b)}{Q(a,b)}}$$$\pi =\lim _{n\to \infty }{\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$

См. также

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