【并行计算】六种方法计算圆周率

并行计算实验之计算圆周率。从传统的蒙特卡洛方法计算圆周率开始,依次给出了计算圆周率的几种算法,并从并行计算的角度对这些算法的性能和特点进行了研究。

计算π的历史

圆周率π可能是科学界内外最广为人知的自然常数了。早在公元5世纪时,南朝宋数学家祖冲之用割圆法将圆周率计算到小数点后7位数字。大约同一时间,印度的数学家也将圆周率计算到了小数点后5位。历史上首个π精确无穷级数公式(即莱布尼茨公式)直到约1000年后才由印度数学家发现。20世纪以来,随着计算机技术的快速发展,π的精度也在极速提高。截至2019年,π的十进制精度已高达1013位。虽然几乎所有的科学研究对π的精度要求都不会超过几百位,当前依然有许多科学家和爱好者为了打破记录、测试超级计算机的计算能力和高精度乘法算法等原因,不断地向π的更高位数发起挑战。

圆周率有许多优美的数学性质,由此产生的计算圆周率的公式和算法也有成千上万。借《并行计算》实验课的机会,我对其中比较有趣和有意义的一些算法进行了一波研究,实现了代码并进行了一些测试。

方法0:蒲丰投针与蒙特卡洛方法

在介绍蒙特卡洛方法之前,先讲一个概率学上的故事:

在1777年,法国数学家蒲丰(也有译作布丰)设计了一个投针试验,实验方式如下:
1) 取一张白纸,在上面画上许多条间距为a的平行线。
2) 取一根长度为l(l≤a) 的针,随机地向画有平行直线的纸上掷n次,观察针与直线相交的次数,记为m。
3)计算针与直线相交的概率.

蒲丰经过计算,得出这个概率是: 。由此蒲丰惊奇地发现:有利的扔出与不利的扔出两者次数的比,是一个包含π的表示式.如果针的长度等于a/2,那么有利扔出的概率为1/π.扔的次数越多,由此能求出越为精确的π的值。

一些数学家在这个问题上进行了实验,结果如下:

在概率学上,蒲丰投针实验是第一个用几何形式表达概率问题的例子,他首次使用随机实验处理确定性数学问题,为概率论的发展起到了一定的推动作用。

那么在计算机中,我们也可以借鉴蒲丰投针的思想,利用随机事件发生的频率来预估出圆周率的值。当然,由于计算机能够很方便地生成随机数,没有必要使用投针的方式,只需随机生成点就可以计算出圆周率。以下描述一种简单的算法:

首先在(0,1)中生成两个随机数x,y,令(x,y)为坐标系中的一个点,那么它将均匀地分布在(0,0),(1,1)的正方形内。然后统计其中到原点的距离小于1的点的数量。这些点分布在以原点为圆心、半径为1的四分之一个圆上,则其占点的总数的比值的期望为π/4。因此进行多次随机,然后将比值乘以4就可以得到近似的圆周率值。

算法的执行过程如图所示:


这个算法有一个明显的好处:它天生支持并行或分布式计算,因为每次投掷的试验都是独立的,可以在很多台设备上分别进行随机试验,之后统计结果即可。

方法1:面积积分

上面的方法已经可以在1秒内计算出圆周率的前8位小数,但是如果想要结果更加精确,难度很大。因为这个算法完全依靠随机次数的增加来逼近理论上的概率,而大数定律的收敛很慢。那么既然上面的算法是通过四分之一个单位圆的面积来得到π的值,我们可不可以直接计算圆的面积呢?

由这个思路引出了计算圆周率的第二个方法:利用面积积分。能够得到π的函数很多,因此这个思路得到的算法也不唯一。这里选用了反正切函数y=arctan(x),利用公式 ,可以得到下图中曲线与坐标轴之间围成的面积为π/4。


在计算机中,计算图形的面积一般可以用划分+近似为梯形面积+求和的方式来计算。具体算法为:将(0,1)的区间等分为n个部分,如第一部分为(0,1/n]、第二部分为(1/n,2/n]等,然后分别计算x坐标在这个区间范围内对应的图形的面积。由于此时对应的图形很“细”,可以近似地看做一个梯形或一个矩形,则可以套用梯形或矩形的面积计算公式求得结果。最后将这些矩形的面积相加,就得出了不规则图形的面积。可以看出,这个算法的精确程度取决于一开始的划分数量,划分得越细,计算的精确度就越高。

同时这个算法也支持并行计算:在程序的开始对梯形的划分进行预分配,然后各个进程分别计算自己被分配到的梯形的面积,最后在进行求和即可。

并行程序的代码如下:


在单机4进程的并行环境下运行程序的结果如下图所示。可以看到,经过1e9次分割,误差2.49e-14,耗时0.589267秒,计算速度约为27位/秒。


方法2:幂级数

除了利用将函数积分转化为分段求和的方式来逼近圆周率之外,数学上还有不少可以用来计算π的公式,比如利用幂级数对arctan(x)进行展开可以得到:


对这个式子进行简单的循环求和,也可以逐渐逼近圆周率的值。

同时这种方法也可以实现并行计算,首先由主节点确定要计算的项数,将计算任务分配到各个节点,然后上多个节点可以分别计算式子的不同部分,然后将结果相加即可。

以下是这种方法的代码:



在单机4进程的并行环境下运行程序的结果如下图所示。可以看到,经过2e9次分割,误差4.98e-10,耗时0.877899秒,计算速度约为14位/秒。


方法3:改进的幂级数

以上两种方法本质上都是通过逼近一个极限为π的数列来计算圆周率的。上面的两个数列每一项衰减得比较慢,要精确到10-N大致需要计算2*10N项。

对于幂级数而言,当x越接近于0时,收敛越快。比如在上面的例子中,我们选取的x=1,离0有相当的距离。如果对arctan(x)在x=1/5处展开,可以得到级数:


以上级数由于含有1/5和1/239的指数运算,收敛速度很快。事实上,当n=4时,就有1/9*59<10-6;而右边收敛的更快。

采用和方法2相似的并行计算方式,我们就可以很快地计算圆周率了。从这种方法开始,由于计算出的圆周率的精度远超double或long double的精度,为了测试这种方案的效率上限,我使用了boost库中的multiprecision::cpp_bin_float这个高精度库,用来进行高精度数字的计算。该高精度库使用了FFT加速的大数乘法等运算,效率较高。

由于高精度的变量类型不在Mpi的基本变量类型范围内,因此改用了Mpi_Gather函数对其进行收集后手动叠加。此外,由于高精度场景下加法、乘法、除法等运算的复杂度不再是O(1),导致求和式中的每一项的计算所需的时间复杂度也不相同。因此并行时采用了均摊的思想,每个程序计算第(i+nk)项,其中i为程序的rank,n为线程数量,k为正整数。

本方法的代码如下:



在单机4进程的并行环境下运行程序的结果如下图所示。可以看到,经过1600次迭代,误差1.17666e-1999,耗时1.6秒,计算速度约为1250位/秒。


进行了更多数据的测试,结果如下:

720次循环,误差为6.28261e-999(精确到1000位),花费时间0.503472秒。
1440次循环,误差为1.17666e-1999(精确到2000位),花费时间1.447734秒。
2160次循环,误差为5.76950e-2999(精确到3000位),花费时间为4.64656秒。
3500次循环,误差为6.62508e-5000(精确到5000位),花费时间为24.713461秒。

通过测试结果可以看出,随着循环次数的增多,循环次数(级数的累加次数)与圆周率计算的位数呈线性增长,增长率约为0.72(循环/位)。但可能是由于高精度运算比较耗时,如对高精度的指数、除法等运算的复杂度较高,可以看到随着位数的增加,运行时间的增长非常迅速,可能为O(n3)。因此使用这个方法来计算更多位数的圆周率比较困难。

方法4:BBP公式

以上这些圆周率计算方法都是基于简单的数学结论,已经能够获得不错的计算结果,但还有很大的改进空间。随着人们对计算圆周率更高位数的不断追求,其算法也在不断地被发明和改进。这个链接给出了一个圆周率位数的世界纪录年表,可以在这里看到各位前辈在这个问题下做出过的努力:Chronology of computation of π。下面研究并介绍了其中的几个典型的“超级圆周率算法”。

BBP公式,全称Bailey–Borwein–Plouffe formula,发现于1995年,以三位发表者的名字命名。


这个公式的特殊之处就在于公式中的(1/16)^k。式子的前半部分是一个单调递减、极限为0的式子,这意味着对于任意正整数k,求和式的第k+1项起每项都(远)小于(1/16)^k。这个性质使得该算法有一个非常神奇的特点:它可以跳过圆周率前面的位数直接计算目标位数,比如直接计算圆周率的第100000位,而不用耗费内存来存储和计算前99999位数字。

若想要计算圆周率的第n位(此处指十六进制下的第n位,后面会考虑如何计算十进制),对该公式进行如下变形:


然后计算右边的小数部分。其中右半部分可以写成:


其中mod指求余运算。这样,前半部分可以使用快速幂算法加速运算,复杂度为O(n*log(n));后半部分收敛速度很快,可以通过计算循环的前100项来代替。所以整个过程的复杂度为O(n2log(n))。

该算法的源论文地址为:math.CA/9803067

由于本算法是基于分别计算圆周率的各位数码,因此很适合改造为并行计算。在实验中我首先实现了BBP算法,然后将其改造成了相应的并行计算程序。

在细节处理方面,其中值得注意的是对运算结果进行按照16进制取模的过程中,为了保证其浮点数存储没有误差(因为分段运算,一点点误差都会导致后面的数值白算),先后尝试了两种方法:

double solve(int id){
	double res, series(int, int);
	res = 4.*series(1, id) - 2.*series(4, id) - series(5, id) - series(6, id);
	res = res - floor(res);
	///方法一
	//auto temp = (*((long long*)(&res)))&0xFFFFFFFFFFF00000L;
	//return *(double*)&temp;
	//方法二
	return floor(res*4294967296.l)/4294967296.0;
}

第一种方法是对double进行截断,即根据C++中浮点数的存储规则,对其超出精度的部分进行截取,相当于在16进制下对16^-9取模;另一种方法是利用floor函数进行截断。经过试验,第二种方法比较稳定,在计算到10000位时都没有遇到问题,而第一种方法会产生误差。

以下是BBP程序的完整代码:


程序的运行截图如下:


测试结果:
计算到1000位,所需时间0.812535秒;
计算到2000位,所需时间1.839296秒;
计算到5000位,所需时间5.365835秒;
计算到10000位,所需时间14.599410秒。

可以看出,用这种方法计算圆周率的效果较好,速度增长率接近线性(慢于线性的原因猜测是内存开销变大和高精度操作导致的缓慢),而且很短时间内计算出了前10000位的数值。

此外,这种方法的一个显著优点是它可以从圆周率的中间开始计算,如直接算出圆周率的第100000位而不需要知道前面的数字,既能很方便地转化为并行计算,又能提供一种计算圆周率较高位数的捷径。

阅读参考资料得知,这种算法可以被用于验证圆周率的超高位数计算的结果是否正确,方法是利用其直接计算某一位的性质,对计算结果进行随机抽测。

方法5: Chudnovsky 算法

Chudnovsky 算法在1988年被乔纳斯基兄弟发明,是至今圆周率计算位数世界纪录保持者使用的算法。Google Cloud在2018年九月-2019年一月的五个月里,用该算法计算出了圆周率的前31.4兆位小数。

该算法的详细证明可参见这篇论文:JSTOR 40391165。它基于这个收敛速度极快的广义超几何级数:


对式子进行化简,使其便于计算机计算,可以得到:


因此算法可以写作一个迭代式:


其中


由于算法基于迭代,难以并行化。因此简单地对该算法实现了非并行的代码如下:


程序的运行截图如下:


该算法在三分钟左右计算出了圆周率的前25000位结果。由于这是在单核单进程的运行环境下进行的运算,达到这样的效率,比起BBP算法又是进步了不少。

事实上,虽然该算法从逻辑上难以并行化,但查阅资料得知,大整数乘法的FFT运算本身其实是可以并行化的,甚至可以放在GPU、TPU等设备中进行加速运算,达到更高的计算速度。如果以后有机会在这些设备上编程时,也许我会想办法在上面写一个Chudnovsky 算法,测测看它们的性能。

结论

经过这些实验可以发现,同样一个计算圆周率的问题,不同的数学原理得出了不同的算法,其结果和速度也会天差地别。
同时,随着计算位数的提升,不仅是公式的收敛速度、算法的时间复杂度在限制着计算的速度,高精度运算的时间消耗、内存消耗也开始出现瓶颈。查看世界纪录计算数亿上兆位圆周率的计算报告可知,在更高位数上计算时,RAM、Cache、磁盘、RAID、cpu指令集等因素也对结果有很大影响,都被优化到了极致。

此外,通过这次实验我也体会到了并行计算的强大之处:从蒙特卡洛算法到积分分割,再到级数计算和BBP算法计算数码,各种不同的算法通过各种方式被并行化,多个节点分别计算部分任务,然后再Gather到主节点上,达到了数倍的性能提升。

算法的优化是永无止境的。对圆周率这个超越数的探索过程,也是对计算机计算能力的探索过程。从公式、算法到并行计算、硬件优化,也许这就是圆周率计算的魅力所在吧!

参考资料

1. 《并行计算实验指导书——圆周率π的计算》 – 北邮软件学院 卢本捷

2. 那些打破圆周率小数位计算的记录是怎么判断计算得正不正确的? - 知乎 酱紫君

3. Pi Formulas, Algorithms and Computations - Fabrice Bellard

4. 圆周率 - Wikipedia

5. BBP Formula - Wikipedia