高精度 π值计算的若千问题

高精度 π值计算的若千问题
李学武
(天津师范大学计算机科学系)
(根据天津师大学报(自然科学版) 第16卷第4期1996年12月同名论文删改)

摘要 本文阐述了高精度π值计算的重要意义，并介绍了作者在π值计算、π值理论研究及算法优化等方面的工作进展。
关键词 圆周率 π 高精度 算法的优化 计算公式

0 引 言
1995年11月，我们用PASAL语言编制程序，在一台486DX/l00微机上计算π(圆周率)值，结果如下:
计算的小数位数: 1000位 10000位 50000位
运行时间 8 s 10 min 6 h 22 min
所得结果与文献[1]所刊登的数据完全一致，早在1989年，美国的G.V.库德诺夫斯基与D.V.库德诺夫斯基两兄弟分别在巨型计算机CRAY-2与IBM3090上，已经将π值计算到小数点后10亿多位[2]。针对这一现实，本文将着重阐述我们的工作所具有的多方面的意义，最后通过两个命题给出作者在π值算法优化方面的部分成果。

1 π值计算的现状与意义
π的本原含义是圆周率，即圆的周长与直径之比。在大多数科学计算中，所用的π值一般不超过10位小数，如果用一个球面把目前人类所探测到的宇宙空间包起来，再计算该球面大圆的周长，为了使误差不超过一个氢原子的直径，只需使用π的39位小数[3]，迄今我们还不知道还有什么实际问题的计算，需要比这更精确的π值。
然而，从古至今，人们始终没有中断对更精确的%值的渴求。有关对π的研究与计算的历史，将另文阐述，这里我们只通过下表列举儿个重要的事件。有关数据大都依据原始资料进行了订正。

年代 作 者 计算位数 备 注
公元429
--500之间 祖冲之(中国) 7位 这一成果领先世界约1000年
1706 J.Machin（英） 100位 [4]
1949 G.W.Reitweisner(美) 2037位 ENIAC机，运行70h[4]
1961 D.Shaks & J.W.Wrench（美） 10万位 IBM7090机，运行8h43min[1]
1973 Guilloud&Bouyer(法) 100万位 CDC7600机[4]
1987 金田康正（日） 134217000位(1.3 亿位) NEC SX-2机[4]
（英译名为Kanada）
1989 1989 G.V.Chudnovsky & 1,011,196,691位（10亿位） IBM3090机，CRAY-2机[2]
D.V.Chudnovsky(美)

尽管π值目前已计算出10亿多位，但对π和π的计算的研究依然没有终止，库德诺夫斯基在完成他们的计算后曾指出，关于π，还有很多工作可做。事实上，1989年以后，在一些同际上颇有影响的学术刊物上，仍不时出现有关π的论文。据此，我们认为，高精度π值的计算，至少具有以下几方

面的意义:
1.1检验计算机硬件和软件的性能
例如，我们在486DX/100微机上计算5万位π值，大约需要6h，运行中，任何一个微小的机器故障，都将导致错误的结果。因而，这对于计算机的运行速度，以及运行的可靠性，稳定性，都是一个很好的检验?为此，我们把有关计算程序与参考数据制成应用软件，并请某计算机销售公司协助惟广，该软件尤其适合在购买新机器时使用。与此同时，我们还把PASCAL源程序分别在Turbo Pascal V6.0与Turbo Pascal V7.0环境下编译，结果发现后者比前者大约节省10%的运行时间。用这一方法可对各种高级语言的编译系统进行检验。1986年，将π值计算到2900万位的贝利(Bailey)指出，他做这一计算的直接动因，就是为了检验超级计算机CLAY-2的硬件环境、FORTRAN编译系统与操作系统的可靠性与完善性[5]。许多文献都指出，高精度π值计算可作为检验计算机性能的一个标准计算。
1.2 关干π的小数位的随机分布性质的研究
这一点在理论上与应用上都具有重要价值。1771年，兰伯特(Lambert)证明了π是无理数，1882年，林德曼(Lindemann)进一步证明了π是超越数(即它不是任何一个有理系数代数方程的根)，但关于π的各位小数的随机分布的性质，目前还主要依赖于统计分析，计算的位数越多，统计分析的结果就越可靠。随机数在多种应用与研究中都占有重要的地位，目前许多重要的研究课题(如利用蒙特卡洛方法研究中子迁移)都需要分布很均匀的随机数(其均匀分布性可作这样直观的理解:在一个相当长的由0-9这10个数字组成的序列中，每个数字 (如'2')出现的频率大约是十分之一，每个相连两位数字(如'37')出现的频率大约是百分之一，…等等)。为此我们设计了一个用蒙特卡洛方法作定积分计算的试验，分别采用QBASlC，Turbo C，Turbo Pascal等语言系统提供的随机数和π的小数位为随机数进行计算，初步结果表明，用π的小数位计算的效果最好，最后的报告将在实验全部完成后给出。文献[2]也曾指出，从理论上讲，由--个序列产生的随机游动，应满足由这一序列的布朗运动所确定的一种多重对数定律，他们采用多种常用的伪随机数序列进行模拟，都不能通过这种检验，惟独用π的小数位产生的随机游动得到了很好的吻合?
在六十年代以前，人们使用随机数，上要依赖于布丰表等人为的一些随机数表，有了计算机之后，多种高级语言都提供了用某种算法产生的随机数(即伪随机数)，它只占用很少空间，但有两个致命的弱点:周期不够长;分布不够均匀。这对于一般的应用影响不大，但对某些重要的科学计算，将会产生不能容忍的误差。因此，在必要时，

用π或由π生成的其它数作为随机数，似乎更为理想。在外存中存放几十万位或数百万位π值，耍占用一定的存储空间，但与目前正方兴未艾的多媒体技术所需空间相比，这只不过是沧海一粟。
1.3 促进相关学科发展
π值的各种现代算法与高斯的算术儿何平均数列，以及印度数学家拉玛努贾(Ramanujan，1887-1920)的模方程理论有着密切的联系，这要涉及到解析数论、特殊函数论特别是椭圆函数理论等较深奥的数学分支。在拉玛努贾去世后三十多年，人们发现了他的--个笔记本，里面有数百个公式，大都没有证明。令人惊讶的是，有些公式在他去世若干年后，才陆续被人们发现。而大多数却不为人所知，甚至至今仍没有得到证明，而其中不少公式都可在数学上、物理上得到应用。有关这方面的研究已有很多成果，其中不少都与π有关[6][7][8]。
在现代算法申，普遍认为效率最高的是一种四阶收敛算法(详见本文第2部分)，其初值取
a0=6-4√2，最后ak收敛于1/π。根据椭圆函数的理论，这--类迭代算法的极限都是初值的函数，因而，--些文献指出，π要算到多少位，初值中的√2就要算到多少位[5]，我们的计算表明，直接用1/π的前若干位作初值，确实不收敛于1/π，但也不一定需要把√2计算得与π一样精确。如果这一点在理论上得到证实，不仅可在很大程度上减少运算量(前几步可只用较少的位数迭代)，而且必将深化人们对迭代方法的误差传播理论的认识。
高精度π值的计算，要涉及多方面的计算机算法。例如，大多数计算都要使用高精度乘法，这就要设法利用快速Fourier变换或快速数论变换，它可将一个乘法的位复杂性由0(n2)降至O(nlog2)，此外，还要涉及各种并行算法，以便充分利用并行 (向量)计算机所提供的功能。例如，贝利在计算π值时设计的并行算法，比串行算法快了20倍[5]。这些工作的意义都远远超出了对π值本身的计算。中科院计算中心以孙家昶教授为首的一个研究小组，于1991年1月，在一个由5片Transputer组成的并行系统上，将n值计算到小数点后100万位，用了102 h。该计算的直接目的，是为了研究"任意精度计算"，这是纳人国家 “七五”规划某项目中的一个子课题[9]。(据最新消息 (国际电子报，1995.11.2) ，该组在4台SGIPowerChallenge可移植环境PVM上，在1995年，用18 h并行计算π值到2016万位)。
1.4 π值可用于密码编制
π作为超越数，特别是利用它的均匀分布的性质，可在密码研制上发挥独特的作用，它可看作是一个取之不尽的“码源”，源文中每个字符使用密码的个数可与它出现的频率成正比，从而使用传统的统计分析方法几乎无法破译。唯一的

缺陷是密文的长度可能是源文的几倍，但按照现代的通信技术，这已不成为障碍，我系90级李刚同学的毕业论文，在这方面作了很有意义的尝试。
(待续)
2 π值计算公式

π值的计算公式(或方法)可找到上百种甚至更多一些，仅此一项就可写成一本书。但不是每个公式都可用于高精度计算 (待续)