Chudnovsky 公式计算圆周率点后百万位

算法引入

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

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

Update on 2022/04/24 修改了一些细节。
Update on 2022/04/30 修改了一些细节。

$$\require{color}$$

开门见山

\[\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)\) 以下的乘法算法包括 Karatsuba, Toom-Cook, 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)\cdot{}Q(a,b)\cdot{}S(a,b)\) 代入就可以证明了。

至此，我们可以得到一个重要的结论：我们只要取 \(m=(a+b)/2\) 分治去求 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(k!\cdot 640320^k)^3} \]

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

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

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

下面这一步

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

是这样推出来的：

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

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

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

\[\therefore\frac{(6k)!}{(3k)!}=\prod_{i=1}^k24(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