圆周率之随笔（3）

 

今天是圆周率节，赶工怒发博文！

在古代，为了计算高精度的圆周率，数学家们往往要用纸和笔辛苦计算几个月甚至几年。而今天我们点点鼠标就可以悠闲的分分秒把圆周率计算到几百万甚至上亿位。这一章我们来看看，在现代计算机上，如何在最短时间里计算出最长的圆周率。

不少人想到计算圆周率可能会联想到SuperPi这个邪恶小程序。这个程序是日本人金田康正在1995年写的，可以在家用电脑Windows系统上把圆周率计算到2**32位。国内很多的装机商攒机商，都用这个程序来测试新装机的超频性能。这其实是很不科学的。首先这个程序是1995年开发的，距离今天已经有快20年，很多新CPU的特性都没有利用上。第二，这个程序是用单线程开发的，对于现代CPU多核的特性也没法支持。第三，单一的用Superpi做超频测试，也容易让硬件厂商进行有针对性的作弊。现在其实有比较新的Hyper PI 可以做为Superpi的替代品。

Superpi 和 Hyper PI 远远不是现在在家用电脑上计算圆周率的主流和先锋。他们采用的是高斯-勒让德算法。其算法也早已过时，有更好的可以替代。

由于圆周率的计算算法实在是五花八门，如果都总结下来的话，作者可能可以骗几十章博客更新，读者可能也懒得再看直接关窗口。那么这里只说一个目前为止最先进和最好的算法。这就是楚德诺夫斯基(Chudnovsky)算法。

楚德诺夫斯基是生于乌克兰的美国数学家，1992年纽约人杂志把楚德诺夫斯基评为世界上最厉害的数学家之一。1987年他发明的计算圆周率的楚德诺夫斯基算法，一直到今天仍然保持着圆周率的最佳计算记录。

该算法是如何来的，请看专业论文，这里直接给出公式：

看上去超级复杂是不是？其实不然，通过简化，代码只有短短的几行。该算法的特色就是，整个计算过程只用最后做一次除法和一次平方根运算。其余全都是加法和乘法！

a = lambda i: (1,1,13591409) if i==0 else ((6*i-5)*(2*i-1)*(6*i-1),i**3*10939058860032000,        (6*i-5)*(2*i-1)*(6*i-1)*(13591409+545140134*i)*(-1 if i&1 else 1))b = lambda x,y: (x[0]*y[0], x[1]*y[1], y[1]*x[2]+x[0]*y[2])bs = lambda i,j: a(i) if j-i==1 else b(bs(i,(i+j)/2),bs((i+j)/2,j))sqrtC = lambda digits: isqrt(1823176476672000*10**(2*digits))def isqrt(n):    x,d = 2**(n.bit_length()-1),1    while d:        d = (n/x-x)/2        x += d    return xdef pi(digits):    N = int(digits/14.1816474627254776555+1)    P, Q, T = bs(0, N)    return (Q*sqrtC(digits)) / T

这段代码在Python 2.X上可执行，运行print pi(100)看看？圆周率的前100位就出来咯。

对于这段代码稍作一些解释：上文看似复杂的公式，其实被拆成了两个函数：a 和 b。整个圆周率计算被简化成了一个mapreduce操作：reduce(b, map(a, range(N))，也就是 N 次 a 运算和 N-1 次 b 运算。为了继续优化，我们设计了 bs 函数，做的是二分切割，将每次计算区域切割两份，然后分头计算，最后做归并。这个计算量其实是一样的，还是 N 次 a 运算和 N-1 次的 b 运算。那么二分切割为什么会运行的更快呢？其主要原因是进行二分切割可以将高精度的大数乘法尽量拖后，执行的层数越低，数字的长度才成指数倍提高。

楚德诺夫斯基算法号称是目前计算圆周率最快的算法，然而这段代码在笔者的电脑上运行计算 pi(10000) 居然要 60 秒之久。问题出在哪呢？经过简单的 benchmark 我们发现，其主要时间消耗在了 isqrt 函数里。由于 python 默认没有高精度整数的平方根操作，所以这个用牛顿拉普森法模拟的平方根操作很慢。

为了快速求平方根，下面我们祭出高精度运算的大杀器 gmp。同时采用python可以用的 gmpy 来加速平方根操作。经过修改的代码如下：

from gmpy2 import *a = lambda i: (1,1,13591409) if i==0 else (mpz((6*i-5)*(2*i-1)*(6*i-1)),mpz(i)**3*10939058860032000,        mpz((6*i-5)*(2*i-1)*(6*i-1))*(13591409+545140134*i)*(-1 if i&1 else 1))b = lambda x,y: (x[0]*y[0], x[1]*y[1], y[1]*x[2]+x[0]*y[2])bs = lambda i,j: a(i) if j-i==1 else b(bs(i,(i+j)/2),bs((i+j)/2,j))sqrtC = lambda digits: isqrt(1823176476672000*mpz(10)**(2*digits))def pi(digits):    N = int(digits/14.1816474627254776555+1)    P, Q, T = bs(0, N)    return (Q*sqrtC(digits)) / T

经过修改的代码也更短更简洁了，只有短短的10行。让我们测试一下速度如何？(需要先安装 gmpy2 的 python 包才可执行)

位数时间10,0000.04s100,0000.11s1,000,0001.07s10,000,00017.76s

20秒内，在笔记本上轻松将圆周率计算到了 1000 万位！1987年在巨型超级计算机上才创造的圆周率世界纪录，在今天的任何一台个人电脑上，都可以轻松碾压了。

然而我们还不满足。这个代码依然没有把楚德诺夫斯基算法的优势发挥出来。楚德诺夫斯基算法最大的优势是可以很好的并行化。在多核时代，单线程运行的代码实在没法体现出CPU的优势出来。

然而由于python的GIL的设计问题，多线程的代码无法很好的利用多核。因此我们用多进程的方法实现了楚德诺夫斯基算法的并行化。下面是在4核CPU的笔记本的运行环境下比较不同进程数量的程序性能。计算位数依然是1000万位：

计算核心(进程数)时间单核版本17.76s118.58s212.66s310.39s49.88s59.79s610.17s

可以看到，通过并行化代码，利用CPU的多核，我们将性能提升了1倍！在计算核心继续增加超过CPU本身的4核之后，性能不升反降，这也是预料之中的事。

有人会问，为什么用4个计算核心提升的性能不是4倍？这里其实牵扯到并行化算法中比较复杂的问题。4核心最理想的性能当然是4倍的提升，但是由于算法本身限制，以及任务细分之后之间的依赖关系，计算核心之间的数据交换，导致最理想情况很难做到。算法中的任务划分也非常有学问，分的太粗会对多内核利用率不高，分的太细又会大大增加进程/线程间的通讯和同步的消耗。这里就不做更深入探讨了。

在这里我们不得不提到 gmp 官方给的一个 C 语言实现的楚德诺夫斯基算法的圆周率计算程序。这个程序对算法细节进行了大量的优化和微调，代码长度接近800行。其中 pgmp 是利用多线程做的并行楚德诺夫斯基算法。有兴趣的同学可以阅读一下代码了解更多优化细节。

单核C: gmp-chudnovsky 并行C: pgmp-chudnovsky 作者的山寨并行Python版: pi-chudnovsky

下面我们用作者自己的无优化山寨 Python 版本和 gmp 官方给的 C 语言版本的圆周率计算程序进行一下性能测试比较。测试环境是一台8核的Amazon c1.xlarge虚拟机。测试位数是圆周率的前1亿位(100,000,000)。

计算核心(进程数)时间(C版本)时间(Python版本)单核版本242.3s397.0s1309.4s419.4s2192.3s226.9s3159.9s168.1s4137.2s135.1s5137.5s116.7s6130.8s105.3s7129.7s102.4s8126.4s99.8s

可见尽管作者的 python 实现的山寨程序尽管算法所用的数据结构优化不够，因此在计算核心较少情况下表现不佳，不过在并行化方面还是做的很好，充分利用到了多核CPU的每一个计算核心。100多行的代码用 python 操作符重载的特性来简化了原本应该很复杂的并行算法的编写。

当然，这些还远远没有到极限。关于圆周率计算的竞赛自古到今一直在进行，在互联网时代越发的激烈。计算竞争的平台也从之前的大型机小型机转到了个人计算机。这里介绍一个著名的圆周率项目叫y-cruncher。

y-cruncher是目前的圆周率计算世界纪录保持者，经作者测试，在刚才的同样测试环境下，计算圆周率1亿位只需要 36.9 秒。y-cruncher 的基本圆周率算法依然还是楚德诺夫斯基算法，但是其代码长度已经到达了17万行！由于优化的足够好，其主要的计算上限还是在内存的大小，y-cruncher的实现主要精力也是在底层内存交换、CPU底层特性的利用以及基础高精度算法的优化上。具体细节大家还是去看他的主页吧。

最后，这篇断断续续写了很久的博客总算完成了。由于篇幅关系，关于博大精深的圆周率以及圆周率文化只是浅尝辄止，希望大家看完以后也能和我一样从此喜欢上深邃的圆周率。