算法引入

这篇文章将介绍目前最快的用于计算圆周率的公式之一——Chudnovsky 公式，以及能显著加快其计算速度的 binary splitting 算法。

提示：本文公式较长，使用手机阅读的读者可以尝试横屏阅读（可能需要刷新一下）。

开门见山

\[\frac{1}{\pi{}}=12 \sum_{k=0}^{\infty{}}\frac{(-1)^k \cdot{}(6k)!\cdot{}(13591409+545140134\cdot{}k)}{(3k)!\cdot{}(k!)^3\cdot{}640320^{3k+3/2}} \]

详细证明可参考 [3]。

前置知识

高精度小数的表示和运算

在运算过程中，我们规定所有的小数均为定点小数。假如规定小数点后的位数是 n，那么运算过程中，将所有数放大为原来的 \(10^n\) 倍，存储在 Python 内置的 int 类型中，就可以进行计算了。两个小数的加减法等价于两个 int 的加减法；两个小数的乘法等价于两个 int 相乘再除以 \(10^n\)；两个小数的除法等于被除数乘 \(10^n\) 再与除数相除；小数的开方等价于乘 \(10^n\) 再开方。

朴素算法

这个公式实在是太复杂了，我们来先处理一下。（下面会有很多公式，但是计算过程并不难理解）

\[\frac{1}{\pi{}}=12 \sum_{k=0}^{\infty{}}\frac{(-1)^k \cdot{}(6k)!\cdot{}(13591409+545140134\cdot{}k)}{(3k)!\cdot{}(k!)^3\cdot{}640320^{3k+3/2}} \]

\[=\frac{12}{640320\sqrt{640320}} \sum_{k=0}^{\infty{}}\frac{(-1)^k \cdot{}(6k)!\cdot{}(13591409+545140134\cdot{}k)}{(3k)!\cdot{}(k!)^3\cdot{}640320^{3k}} \]

\[=\frac{1}{426880\sqrt{10005}} \sum_{k=0}^{\infty{}}\frac{(-1)^k \cdot{}(6k)!\cdot{}(13591409+545140134\cdot{}k)}{(3k)!\cdot{}(k!)^3\cdot{}640320^{3k}} \]

\[a_k=\frac{(-1)^k\cdot{}(6k)!}{(3k)!\cdot{}(k!)^3\cdot{640320^3k}} \]

\[b_k=\frac{(-1)^k\cdot{}(6k)!\cdot{}k}{(3k)!\cdot{}(k!)^3\cdot{640320^3k}} \]

\[a=\sum_{k=0}^{\infty{}}a_k \]

\[b=\sum_{k=0}^{\infty{}}b_k \]

那么我们可以发现：

\[\pi{}=\frac{426880\sqrt{10005}}{13591409\cdot{}a+545140134\cdot{}b} \]

我们只需要把 a, b 求出来就可以了。而我们又发现

\[b_k=k\cdot{}a_k \]

\[a_0=1 \]

\[\frac{a_k}{a_{k-1}}=-\frac{(6k-5)(6k-4)(6k-3)(6k-2)(6k-1)6k}{640320^3k^3(3k-2)(3k-1)3k} \]

\[\implies{}a_k=-\frac{24(6k-5)(2k-1)(6k-1)}{640320^3k^3}a_{k-1} \]

思路渐渐地清晰起来了。下面是 Python 的实现：

# 使用 Chudnovsky 公式计算圆周率小数点后 n 位
def chudnovsky_naive(digits):
    import math
    digits *= 2  # 要先乘2，不然计算出来只有 n/2 位
    one = 10 ** digits
    k = 1
    ak = one
    asum = one # 10 ** digits
    bsum = 0

    # ak=0 时说明 ak 太小了，超过了规定的精度，这时停止计算
    while ak:
        ak = 24 * ak * -(6*k-5) * (2*k-1) * (6*k-1) // ((k * 640320) ** 3)
        asum += ak
        bsum += k * ak
        k += 1
    
    denominator = 13591409 * asum + 545140134 * bsum

    # 注意：要用 math.isqrt() 而非 math.sqrt()。
    # math.sqrt() 是将整数转换为浮点数，结果也是浮点数，
    # 10005*one 这么大的整数作为参数传入会报错。
    # math.isqrt(x) 用于计算 int(sqrt(x))，
    # 对于非常大的整数同样生效且计算速度很快。
    numerator = 426880 * one * math.isqrt(10005 * one)

    return numerator // denominator

我们输出看一下

n = int(input("计算圆周率。输入计算位数："))
# 注意计算出来的值是 int(pi * 10**n)，输出的时候要补上小数点
s = str(chudnovsky_naive(n))
print("3." + s[1:])

它给出了正确的结果。在代码开头的地方，我写了一行 digits *= 2，为什么？这是因为 numerator 大概是 \(10^{n+0.5n}\) 这么一个量级的数，denominator 是 \(10^n\) 量级的数，两者相除会得到一个 \(10^{0.5n}\) 量级的数，也就是说实际算出来的位数只有一半。所以我们要将位数乘2，这样算出来的位数就不会少。

我们再测试一下速度

import time
n = int(input("计算圆周率。输入计算位数："))

ts = time.time()
# 测试速度就不输出了，因为二进制转十进制过于耗费时间
x = chudnovsky_naive(n)
te = time.time()
print("耗时：%d 秒" % (te - ts))

我在一台电脑上测试了一下。

处理器：Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz 2.00 GHz

RAM：16.0 GB

Python 3.9.6 64-bit

计算位数	计算时间（秒）
10	0.0
100	0.0
1000	0.00199127197265625
10000	0.0579071044921875
50000	1.4166526794433594
100000	5.523305654525757
200000	21.368163108825684

这个程序虽然在 22s 以内计算出了圆周率小数点后二十万位，但计算点后百万位还是略显吃力。下面介绍 binary splitting 算法。

binary splitting

注意：binary splitting 算法需要配合 \(O(n^2)\) 以下的乘法和除法算法，否则将毫无意义。\(O(n^2)\) 以下的乘法算法包括 , , FFT/NTT；除法算法包括牛顿迭代，分治（两者都需要配合乘法算法使用）

binary splitting 算法（目前我还未看到中文译名）不仅可以加速 Chudnovsky 公式的计算，还可以加快很多级数的运算，只要级数形如

\[S=\sum_{k=0}^{\infty{}}\frac{a_k\cdot{}\prod_{i=0}^{k}p_i}{b_k\cdot{}\prod_{i=0}^{k}q_i} \]

我们先定义一些东西，之后你就会知道这些东西有什么用。定义

\[S(a,b)=\sum_{k=a}^{b-1}\frac{a_k\cdot{}\prod_{i=0}^{k}p_i}{b_k\cdot{}\prod_{i=0}^{k}q_i} \]

\[P(a,b)=\prod_{k=a}^{b-1}p_k \]

\[Q(a,b)=\prod_{k=a}^{b-1}q_k \]

\[B(a,b)=\prod_{k=a}^{b-1}b_k \]

\[T(a,b)=B(a,b)\cdot{}Q(a,b)\cdot{}S(a,b) \]

不难发现，对于任意 m (a <= m < b)，都有

\[S(a,b)=S(a,m)+S(m,b) \]

\[P(a,b)=P(a,m)\cdot{}P(m,b) \]

\[Q(a,b)=Q(a,m)\cdot{}Q(m,b) \]

\[B(a,b)=B(a,m)\cdot{}B(m,b) \]

\[T(a,b)=B(m,b)Q(m,b)T(a,m)+B(a,m)Q(a,m)T(m,b) \]

前四个公式都很好理解；最后一个公式只要把 T(a,b) = B(a,b)Q(a,b)S(a,b) 代入就可以证明了。

定义了这么多式子又推了这么多式子，最后发现：我们只要取 \(m=\frac{1}{2}(a+b)\) 分治去求P, Q, B, T就可以了，当 b = a + 1 时，P, Q, B, T 都可以直接求出。

\[P(a,a+1)=p_a \]

\[Q(a,a+1)=q_a \]

\[B(a,a+1)=b_a \]

\[S(a,a+1)=\frac{a_ap_a}{b_aq_a} \]

\[T(a,a+1)=a_ap_a \]

根据 T(a,b) 和 S(a,b) 的定义式可以得出

\[S(a,b)=\frac{T(a,b)}{B(a,b)Q(a,b)} \]

且 \(S(0,n)\) 就是这个级数的前 n 项和。

我们看一看 Chudnovsky 可不可以变换成这种形式。答案是可以的。

\[\sum_{k=0}^{\infty{}}\frac{(-1)^k \cdot{}(6k)!\cdot{}(13591409+545140134\cdot{}k)}{(3k)!\cdot{}(k!)^3\cdot{}640320^{3k}} \]

\[=\sum_{k=0}^{\infty{}}\frac{(-1)^k \cdot{}(13591409+545140134\cdot{}k)\cdot{}(6k)!/(3k!)}{1\cdot{}\prod_{i=0}^k640320^3i} \]

\[=\sum_{k=0}^{\infty{}}\frac{{(-1)^k \cdot{}(13591409+545140134\cdot{}k)}\cdot{}{\prod_{i=1}^k(6i-5)(2i-1)(6i-1)}}{{1}\cdot{}{\prod_{i=1}^k(640320^3/24)i}} \]

下面这一步

\[\frac{(6k)!}{(3k)!}=\prod_{i=1}^k24(6i-5)(2i-1)(6i-1),k\ne{}0 \]

是这样推出来的：

\[k=0\implies{}\frac{(6k)!}{(3k)!}=1 \]

\[k \ne{}0\implies{}\frac{(6i+6)!}{(3i+3)!}\div{}\frac{(6i)!}{(3i)!}=24(6i-5)(2i-1)(6i-1) \]

还有一个问题：两个乘积式的下界为1，而我们需要一个下界为0的式子，怎么办？强行把它变成下界为0的乘积式！定义 \(p_0=q_0=1\) 即可。也就是说这个乘积式

\[\prod_{i=1}^k24(6i-5)(2i-1)(6i-1) \]

它等于

\[\prod_{i=0}^k 24(6i-5)(2i-1)(6i-1) \space{}\text{if}\space{}i\ne{}0\space{}\text{else}\space{}1 \]

（这里借鉴了 Python 的写法，另一个乘积式同理）

现在证明了 Chudnovsky 公式中的级数可以用 binary splitting 算法求得。我们整理一下公式：

\[a_k=(-1)^k\cdot{}(13591409+545140134\cdot{}k) \]

\[p_0=1,p_k=24(6k-5)(2k-1)(6k-1) \]

\[b_k=1 \]

\[q_0=1,q_k=640320^3k \]

# 使用 Chudnovsky 公式和 binary splitting 算法计算圆周率小数点后 n 位
def chudnovsky_binsplit(digits):
    import math
    digits *= 2

    # 返回 P, Q, B, T
    # 此时 B=1，所以只返回 P, Q, T
    def binsplit(a, b):
        # 直接求的情况
        if b - a == 1:
            # 特殊情况：a = 0
            if a == 0:
                Pab = Qab = 1
            else:
                Pab = (6*a-5) * (2*a-1) * (6*a-1)
                Qab = 640320**3//24 * a**3
            Tab = (13591409 + 545140134 * a) * Pab
            # 对 (-1)^k 这个因子进行处理
            if a & 1:
                Tab = -Tab
            return Pab, Qab, Tab
        else:
            m = (a + b) // 2
            Pam, Qam, Tam = binsplit(a, m)
            Pmb, Qmb, Tmb = binsplit(m, b)
            Pab = Pam * Pmb
            Qab = Qam * Qmb
            Tab = Qmb * Tam + Pam * Tmb
            return Pab, Qab, Tab
    
    # 要计算多少项
    # Chudnovsky 公式每计算一项，正确位数会增加 14 位，
    # 所以要计算 digits//14 + 1 项
    # 想知道原因的读者可以参考 [3] 中 Page 44 的 Theorem 10.13
    terms = digits // 14 + 1
    P, Q, T = binsplit(0, terms)
    return Q * 426880 * math.isqrt(10005 * 10**digits) // T

在同一台电脑测试一下速度

计算位数	计算时间（秒）	朴素算法计算时间（秒）
10	0.0	0.0
100	0.0	0.0
1000	0.0009963512420654297	0.00199127197265625
10000	0.027961015701293945	0.0579071044921875
50000	0.3704197406768799	1.4166526794433594
100000	1.205124855041504	5.523305654525757
200000	4.421649932861328	21.368163108825684

binary splitting 加速的根本原因

binary splitting 算法比朴素算法快的根本原因是：前者更多地计算大数和大数相乘除，而后者则一直在计算大数和小数相乘。前者可以使用 \(O(n^2)\) 以下的乘除法来进行加速，而后者却完全无法从复杂度上优化。

[2] 这篇论文中给出了 binary splitting 算法的时间复杂度：如果进行 n 位数 * n 位数乘法的时间复杂度为 \(M(n)\)，那么 binary splitting 的时间复杂度为 \(O(M(n)(\log n)^2)\)

进一步加速

Python 内置的 int 类型速度还是不够快，想要追求更高的速度就要换用别的高精度整数库。gmpy 就是很好的选择。gmpy 是 gmp 对 Python 的封装。gmp 由 C 语言和汇编写成，是目前最快的开源高精度运算库。

在命令行输入 pip install gmpy2 即可安装 gmpy。接下来我们把代码中的 int 换成 gmpy2.mpz，把 math.isqrt() 换成 gmpy2.isqrt() 就可以了。

# 使用 Chudnovsky 公式和 binary splitting 算法计算圆周率小数点后 n 位
# 使用 gmpy2.mpz 加速
def chudnovsky_binsplit_mpz(digits):
    import gmpy2
    from gmpy2 import mpz
    digits *= 2
    def binsplit(a, b):
        # 直接求的情况
        if b - a == 1:
            # 特殊情况：a = 0
            if a == 0:
                Pab = Qab = 1
            else:
                Pab = (6*a-5) * (2*a-1) * (6*a-1)
                Qab = 640320**3//24 * a**3
            Tab = (13591409 + 545140134 * a) * Pab
            # 对 (-1)^k 这个因子进行处理
            if a & 1:
                Tab = -Tab
            return mpz(Pab), mpz(Qab), mpz(Tab)
        else:
            m = (a + b) // 2
            Pam, Qam, Tam = binsplit(a, m)
            Pmb, Qmb, Tmb = binsplit(m, b)
            Pab = Pam * Pmb
            Qab = Qam * Qmb
            Tab = Qmb * Tam + Pam * Tmb
            return Pab, Qab, Tab
    
    # 要计算多少项
    # Chudnovsky 公式每计算一项，正确位数会增加 14 位，
    # 所以要计算 digits//14 + 1 项
    # 想知道原因的读者可以参考 [3] 中 Page 44 的 Theorem 10.13
    terms = digits // 14 + 1
    P, Q, T = binsplit(0, terms)
    return Q * 426880 * gmpy2.sqrt(10005 * mpz(10)**digits) // T

再次测试速度

计算位数	计算时间（秒）	使用内置int计算时间（秒）
10	0.0	0.0
100	0.0	0.0
1000	0.0033011436462402344	0.0009963512420654297
10000	0.014216423034667969	0.027961015701293945
50000	0.07479691505432129	0.3704197406768799
100000	0.17659258842468262	1.205124855041504
200000	0.35507917404174805	4.421649932861328
500000	1.0489845275878906	N/A
1000000	2.4204533100128174	N/A

速度提升可以说是一次飞跃！而且经过另一次测试，计算小数点后百万位+输出到文件两个操作花费的总时间仅有 3.284137487411499 秒。至此我们也就完成了计算圆周率小数点后一百万位的任务。如果有兴趣的读者还可以尝试使用多线程来加速。我将计算出来的圆周率上传到了 cnblogs，点击这里就可以下载，结果与网络数据对比无误，读者可以用于检测自己的程序是否出错。

参考文献

[1] Nick Craig-Wood. Pi - Chudnovsky. https://www.craig-wood.com/nick/articles/pi-chudnovsky/

[2] Bruno Haible, Thomas Papanikolaou. Fast multiprecision evaluation of series of rational numbers. https://www.ginac.de/CLN/binsplit.pdf

[3] Lorenz Milla. A detailed proof of the Chudnovsky Formula with means of basic complex analysis. https://arxiv.org/pdf/1809.00533.pdf