Chudnovsky 公式计算圆周率点后百万位 |
您所在的位置：网站首页 › 用matlab求圆周率的程序 › 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 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 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |