前置知识：Binary Splitting

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

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

\(\pi\)

\[\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

\(\zeta(3)\)

这里是公式链接。我们使用第三个公式。该级数每项提供 \(\log_{10}(717445350000)\) 位精度。

\[\zeta(3) = \frac{1}{48}\sum_{n=1}^{\infty}\frac{(-1)^{(n-1)}\,P(n)}{R(n)\cdot\binom{5n}{n}\binom{5n}{2n}\binom{9n}{4n}\binom{10n}{5n}\binom{12n}{6n}} \]

其中

\[\begin{matrix} P(n)= 1565994397644288\,n^{11} - 6719460725627136\,n^{10} + 12632254526031264\,n^9\\ - 13684352515879536\,n^8 + 9451223531851808\,n^7 - 4348596587040104\,n^6\\ + 1352700034136826\,n^5 - 282805786014979\,n^4 + 38721705264979\,n^3\\ - 3292502315430\,n^2 + 156286859400\,n - 3143448000\\ R(n)=n^5\,(2n-1)^3\,(3n-1)\,(3n-2)\,(4n-1)\,(4n-3)\,(6n-1)\,(6n-5) \end{matrix} \]

于是根据上述法则，假设 \(a_n=\frac{1}{R(n)\cdot\binom{5n}{n}\binom{5n}{2n}\binom{9n}{4n}\binom{10n}{5n}\binom{12n}{6n}}\)，则

\[\begin{aligned} \frac{a_n}{a_{n-1}}&=\frac{R(n-1)\cdot\binom{5n-5}{n-1}\binom{5n-5}{2n-2}\binom{9n-9}{4n-4}\binom{10n-10}{5n-5}\binom{12n-12}{6n-6}}{R(n)\cdot\binom{5n}{n}\binom{5n}{2n}\binom{9n}{4n}\binom{10n}{5n}\binom{12n}{6n}}\\ &=\frac{(n-1)^5(2n-3)^3(3n-5)(3n-4)(4n-7)(4n-5)(6n-11)(6n-7)}{270(9n-8)(9n-7)(9n-5)(9n-4)(9n-2)(9n-1)(10n-9)(10n-7)(10n-3)(10n-1)(12n-11)(12n-7)(12n-5)(12n-1)} \end{aligned} \]

这时就可以定义

\[\begin{matrix} p_n=(n-1)^5(2n-3)^3(3n-5)(3n-4)(4n-7)(4n-5)(6n-11)(6n-7),\\ q_n=270(9n-8)(9n-7)(9n-5)(9n-4)(9n-2)(9n-1)\\ (10n-9)(10n-7)(10n-3)(10n-1)(12n-11)(12n-7)(12n-5)(12n-1),\\ r_n=P(n),\color{red}p_1=1,q_1=44,008,272,000 \end{matrix} \]

def compute_apery(digits):
    def polp(n):
        return (-1)**(n-1) * (1565994397644288*n**11 - 6719460725627136*n**10 + 12632254526031264*n**9\
        - 13684352515879536*n**8 + 9451223531851808*n**7 - 4348596587040104*n**6 \
        + 1352700034136826*n**5 - 282805786014979*n**4 + 38721705264979*n**3\
        - 3292502315430*n**2 + 156286859400*n - 3143448000)
    def bs(a, b):
        if b - a == 1:
            if a == 1:
                Pab = 1
                Qab = 44008272000
            else:
                n = a
                Pab = (n-1)**5*(2*n-3)**3*(3*n-5)*(3*n-4)*(4*n-7)*(4*n-5)*(6*n-11)*(6*n-7)
                Qab = 270*(9*n-8)*(9*n-7)*(9*n-5)*(9*n-4)*(9*n-2)*(9*n-1)\
                     *(10*n-9)*(10*n-7)*(10*n-3)*(10*n-1)*(12*n-11)*(12*n-7)*(12*n-5)*(12*n-1)
            Tab = Pab * polp(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(717445350000)
    N = int(digits/DIGITS_PER_TERM + 3)
    P, Q, T = bs(1, N)
    return T / Q / 48

\(\ln 2\)

ln(2) = 144*arccoth(251) + 54*arccoth(449) - 38*arccoth(4801) + 62*arccoth(8749)
ln(3) = 228*arccoth(251) + 86*arccoth(449) - 60*arccoth(4801) + 98*arccoth(8749)
ln(5) = 334*arccoth(251) + 126*arccoth(449) - 88*arccoth(4801) + 144*arccoth(8749)
ln(7) = 404*arccoth(251) + 152*arccoth(449) - 106*arccoth(4801) + 174*arccoth(8749)

这样可以通过计算这四个 arccoth 值来计算四个 \(\ln\) 值。

根据这个网页的写法，可以得到如下代码：

def arccoth(x, digits=1000):
    def bs(a, b):
        if b - a == 1:
            Pab = 2*a+3
            Qab = (2*a+3)*x^2
            Tab = 1
        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 = 2 * math.log10(x)
    N = int(digits/DIGITS_PER_TERM + 3)
    P, Q, T = bs(0, N)
    return (1 + T / Q) / x