拉出圆周率的数字.doc

拉出圆周率的数字
David Austin
关键词： 圆周率, 算法
大自然给了我们许多不容易被人类数字系统表示的重要常数，如π和e。例如，π的前50位数字
π=**********…
只提供了这个重要数的一个逼近，因而我们觉得有必要更深入地了解这样的数。我们也可能会问这些重要的常数之间是否有关联。例如，乍一看，π和e似乎在不同的数学领域出现，所以如果发现它们的值可通过简单的方式联系在一起，这将是非常值得注意的。
在这篇文章中，我们将探讨圆的周长与直径的比π怎样以惊人的精度计算出来。事实上，我们将看到如何深入到π的十六进制表示的内部来计算我们希望的任何位数。我们还将介绍相关的新算法，使我们能够探讨各种数值常数是否存在简单的关系。
虽然现实世界的应用不需要准确到50位数字的π，然而计算π的历史是丰富的，其动机主要是出于我们自己的好奇心。对于计算科学家来说，计算问题提出像“数如何有效地相乘”这样的技术性挑战，他们的解决方案有着广泛的应用。
早期的圆周率计算
古希腊数学家阿基米德是仔细研究π的先驱之一。他的技术从一个周长为2π的单位圆开始，然后用内接和外切正多边形的周长来逼近圆周长。
例如，下图所示的两个正六边形的周长将给出π的上、下界。
通过增加多边形的边数，我们可以得到更好的上下界。
如果令an和bn分别表示外切n-边形（标以红色）和内接n-边形（标以蓝色）的周长，则有
bn<2π<an及2π=limn→∞an=limn→∞bn.
当n=6时，容易计算出周长
a6=43√,b6=6.
本质上利用三角中的半角公式，阿基米德确定出当边数加倍后周长如何变化，从而得到循环公式
a2n=2anbnan+bn,b2n=a2nbn−−−−−√.
用这一方法，我们有上下界
22371<π<227.
1600年左右在荷兰工作的德国数学家Ludolph van Ceulen采用了阿基米德的计术算出了π的前35位，并将它的上下界刻在自己的墓碑上。
不久以后，微积分的发展给了计算π的新方法，我们现在介绍当中的一个。记得几何级数具有形式
1+r+r2+r3+….
即，和式中的每一项是前一项乘上一个固定数r。假如r的绝对值小于1，级数收敛并容易赋值。记和为S，
S=1+r+r2+r3+….
再把左右式乘以r
rS=r+r2+r3+r4+….
第一式减去第二式，得到(1−r)S=1，因此
S=1+r+r2+r3+…=11−r.
如果令r=−u2，则得到
11+u2=1−u2+u4−u6+….
两边积分后就有
arctanx===∫x011+u2du∫x0(1−u2+u4−u6+…)dux−x33+x55−x77+…=∑k=0∞(−1)kx2k+12k+1.
如果让x=1，我们得到通常归功于莱布尼茨的π的表达式
π=4arctan1=4(1−13+15−17+…)=4∑k=0∞(−1)k12k+1.
真正计算中，这不是计算π的特别有用的方法，因为它收敛太慢。这里是它的前几项得到的逼近值：
n
前n项的和
1

2

3

4

5

6

⋯
⋯
1000

1001

如你所见，加了1000项后，我们才得到π的前三位数字。
欧拉找到一个更有用的公式：
π4=arctan(12)+arctan(13),
它可用复数相乘简单解释之。记得复数乘积ab的幅角等于a和b的幅角之和。欧拉的公式则由(2+i)(3+i)=5+5i得到。利用上面反正切函数的级数表达式，就有
π=4(1/2−13(1/2)3+15(1/2)5+…) +4(1/3−13(1/3)3+15(1/3)5+…).
相加两个级数的前十项，，精确到前七位。这个公式和涉及反正切的其他类似公式，最终产生π的前几百个数字。
当然，所有这一切都在计算机问世之前实现。最近，Yasumasa Kanada采用了类似的反正切公式，花了一台超级计算机大约600个小时算出π超过一万亿个的位数。
BBP公式
如果我们想计算π的一个特定数字，以上列出的技术要求我们使用高精度的算术来计算它前面的所有数字。20世纪90年代的中期，一个了不起的新的计算公式被David Bailey、Peter Borwein和Simon Plouffe (BBP)发现：
π=∑k=0∞116k(48k+1−28k+4−18k+5−18k+6).