前置知识：Binary Splitting

如果不知道这是什么的可以参考这个网页和这篇文章。

这篇文章将讲述另一些级数的 Binary Splitting 算法。

级数 1

\[\frac1\pi=\sum_{k=0}^{+\infty}\frac{(42k+5)}{2^{12k+4}}\binom{2k}{k}^3 \]

定义

\[\begin{aligned} a_k&=\frac{\binom{2k}{k}^3}{2^{12k}},\\ S_a&=\sum_{k=0}^{+\infty}a_k,\\ S_b&=\sum_{k=0}^{+\infty}ka_k \end{aligned} \]

那么 \(\frac1\pi=\frac1{16}(42S_b+S_a)\)。接下来考虑使用 Binary Splitting 计算 \(S_a\) 和 \(S_b\)。

\[\frac{a_k}{a_{k-1}}=\frac{(2k-1)^3}{512k^3} \]

这时就可以定义 \(p_k=(2k-1)^3,q_k=512k^3,r_k=42k+5,\color{red}p_0=1,q_0=1\)）。于是

\[\frac1\pi=\frac1{16}\sum_{k=0}^{+\infty}\frac{r_k\prod_{j=0}^kp_j}{\prod_{j=0}^kq_j} \]

于是端上代码：

def compute_pi(digits):
    def bs(a, b):
        if b - a == 1:
            if a == 0:
                Pab = Qab = 1
            else:
                Pab = (2*a-1)**3
                Qab = 512*a*a*a
            Tab = Pab * (5 + 42*a)
        else:
            m = (a + b) // 2
            Pam, Qam, Tam = bs(a, m)
            Pmb, Qmb, Tmb = bs(m, b)
            Pab = Pam * Pmb
            Qab = Qam * Qmb
            Tab = Qmb * Tam + Pam * Tmb
        return Pab, Qab, Tab
    DIGITS_PER_TERM = math.log10(2) * 6
    N = int(digits/DIGITS_PER_TERM + 1)
    P, Q, T = bs(0, N)
    return 16 * Q / T