拉出圆周率的数字

，阿基米德确定出当边数加倍后周长如何变化，从而得到循环公式 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 4.0000 2 2.6667 3 3.4667 4 2.8952 5 3.3397 6 2.9760 ⋯ ⋯ 1000 3.1406 1001 3.1426 如你所见，加了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+…). 相加两个级数的前十项，就得到一个近似3.1415925796，精确到前七位。这个公式和涉及反正切的其他类似公式，最终产生π的前几百个数字。 当然，所有这一切都在计算机问世之前实现。最近，Yasumasa Kanada采用了类似的反正切公式，花了一台超级计算机大约600个小时算出π超过一万亿个的位数。 BBP公式 如果我们想计算π的一个特定数字，以上列出的技术要求我们使用高精度的算术来计算它前面的所有数字。20世纪90年代的中期，一个了不起的新的计算公式被David Bailey、Peter Borwein和Simon Plouffe (BBP)发现： π=∑k=0∞116k(48k+1−28k+4−18k+5−18k+6). 我们稍后描述一下这个公式的发现。首先我们将看到这个公式如何让我们直接计算它的一个十六进制数字，而不必算出前面的所有位数并且不使用高精度算术。 要了解这是如何工作的，让我们考虑计算ln(2)=0.69314718056…的一个类型相似的公式。记得 ln(1+x)=∫x011+udu=∫x0(1−u+u2−u3+…)du=x−x22+x33−x44+…=∑k=1∞(−1)k+1xkk. 故有 ln(2)=−ln(12)=∑k=1∞12kk. 让我们看看ln(2)的二进制数字： ln(2)=.d1d2d3d4…=d12−1+d22−2+d32−3+d42−4+…. 如果我们想计算第n+1个数字dn+1，上式两边乘上2n: 2nln(2)==d12n−1+d22n−2+…+dn+dn+12−1+dn+22−2+…d1d2…dn.dn+1dn+2…. 注意到，我们希望找到的数字dn+1是小数点后面的第一位数字。因此，如果2nln(2)的小数部分小于0.5，则dn+1=0，否则dn+1=1。所以我们只需要计算2nln(2)的小数部分，记为{2nln(2)}。它的计算如下： {2nln(2)}==={2n∑k=1∞12kk}={∑k=1∞2n−kk}{∑k=1n2n−kk}+⎧⎩⎨∑k=n+1∞2n−kk⎫⎭⎬{∑k=1n2n−k(modk)k}+⎧⎩⎨∑k=n+1∞2n−kk⎫⎭⎬. 最后一行中的各项可有效地用著名的二进制幂算法计算。第二个和式的总和小于1/(n+1)，所以如果n很大，它的贡献很小。如果还记得我们只需要确定小数部分是小于或大于0.5，我们只要在第二个和式中计算足够多的项就可解决问题。 回到BBP公式 π=∑k=0∞116k(48k+1−28k+4−18k+5−18k+6), 我们看到，每一项都有的16k能让我们以类似的方式计算π的任一位十六进制数字，而不必计算它之前的数字。 π的十六进制数字从第一万亿零一次开始，实际上是B4466E8D215388C4E014。最近的工作已成功地计算出第一千万亿位(1015)十六进制数字。 PSLQ算法 BBP公式的方法之一是考虑积分 ∫1/2√042√−8x3−42√−8x51−x8dx. 一方面，运用大部分学生所熟悉的部分分式技巧，可以证明这个积分等于π。另一方面，几何级数表明 11−x8=1+x8+x16+x24+…. 这使我们能够将被积函数写成级数形式并逐项积分。把这两个事实放在一起，就给出BBP公式。 但BBP公式是如何被发现的？当然不是通过刚才所描述的积分。相反，BBP公式来自于检测实数集合之间整数关系的所谓的PSLQ算法。 举例来说，如果我们有一组实数x1,x2,…,xn，我们可能会问是否存在整数a1,a2,…,an，其中至少有一个不为零，使得 a1x1+a2x2+…+anxn=0. 对于给定的实数，可能没有一组整数a1,a2,…,an满足上述关系。但如果有，PSLQ算法会找到这样一组整数。 检测整数关系有着悠久的历史。例如，欧几里得在《几何原本》的第十章考虑过这个问题n=2的情形。为了解释欧几里得的算法，我们把两个实数称之为A和B。要求有整数a1和a2使得a1A+a2B=0等价于问是否存在实数C以及整数n和m满足A=nC和B=mC。（在这种情况下，欧几里得说A和B是公度的。） 对正数A和B，欧几里得的著名算法如下运行： 如果A=B，则A和B有整数关系，算法以C=A=B终止。 两个实数的较小者称为B。令A=A−B。 重复上述步骤。 如果该算法终止，那么我们已经发现C，使得A=nC和B=mC。如果我们让 a1=−BC及a2=AC,则a1A+a2B=0, 从而证明一个整数的关系。 算法的两个图示如下给出，其中线段的长度代表了A和B的值。在第一种情况下，我们发现A和B有公度，并得到C。 在下一种情形，A=π,B=e。 经过几个步骤后，该算法还没有终止。当然，这并不意味着如果我们让它运行一段时间它还是不会停止；在实践中，我们还不能确定该算法将无限期地运行下去。然而，我们能注意到，如果算法尚未终止，那么我们可以保证C比A的当前值要小。这反过来又给出整数a1和a2的下限。 自欧几里得的时代起，将欧几里得算法推广到两个实数以上的情形已被欧拉、雅可比、庞加莱等人考虑过。1979年，Helaman Ferguson和Rodney Forcade发现了一个算法，该算法此后被Ferguson、Bailey及Steve Arno精细化成PSLQ算法。 PSLQ算法类似于欧几里得算法，它推广成：如果有一个整数关系，该算法将在发现它时终止。如果算法不终止，就没有整数关系。当然，我们可能永远不知道实践中该算法会不会终止；然而，如果算法运行了一段时间后不终止，我们可以确定关系中整数的上下限。 执行PSLQ算法对自身提出挑战，因为计算机只能使用有限精度的运算。例如，我们为之寻求整数关系的实数在计算机内存中只能用有限多个数字近似。此外，算法中需要的真实运算受制于舍入误差。 正因为如此，PSLQ算法的实施不能够产生明确的整数关系。相反，该算法只建议可能的关系，而整数关系存在的证明则需要与算法独立地提出来。 BBP公式通过应用PSLQ算法找到π和8个常数 Xj=∑k=0∞116k(8k+j),j=1,2,…,8 之间的一个整数关系以及作者称之为“受启发猜测和广泛的搜索”而被发现的。然后该公式通过上述的积分方法证明。 BBP公式使我们能计算π的第n个十六进制数字而不必计算前n−1位。目前我们不知道是否存在类似的公式，使我们以同样的方式来计算任意一位十进制数字。 除了发现BBP公式外，PSLQ算法还找到了其他的一些用途。例如，该算法可用于探讨一个给定的实数是否是代数数（即一个整系数多项式的根）。此外，PSLQ发现了涉及量子场论中费恩曼图赋值所产生的常数的一些恒等式。 顺便说一下，一直从事整数关系探测算法基本数学工作的Helaman Ferguson，也是一个用艺术来表现美和数学理解的雕塑家。 链接: http://www.ams.org/samplings/feature-column/fcarc-pi 作 者: David Austin，Grand Valley State University 翻 译: 丁玖，密执安州立大学博士，南密西西比大学数学教授 校 对: 汤涛，香港浸会大学数学讲座教授