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