圆周率pi的bbp计算公式之详细证明.pdf

圆周率 pi 的 BBP 计算公式之详细证明 Rockins Chen UESTC, Chengdu 4/6/2009 1圆周率圆周率 pi 的的 BBP 计算公式之详细证明计算公式之详细证明 Rockins Chen UESTC, Chengdu 4/6/2009 摘要：本文给出了用于计算圆周率 pi 的 BBP 公式的详细证明主要是对文献[1,2,3]的翻译和 解释 关键字：圆周率 pi BBP 公式 1 求解求解 pi 的历史的历史 第一个尝试严格计算π的人是公元前 250 年的阿基米德，他采用内接多边形和外切多边形逼近的办法确定出π的范围是10133717π UESTC, Chengdu 4/6/2009 2后来发现 Shanks 的计算结果中第 527 位以后的数字是错误的 但这仍然是一个很大的进步 在不断追求计算更高精度π的过程中，一些和π有关的重要问题被解决了17 世纪晚期，Lambert 和 Legendre 证明了π是无理数1882 年，Lindemann 证明了π是超越数Lindemann 的证明彻底地解决了源自古希腊的“化圆为方”问题，即能否仅用直尺和圆规将 一个圆变为一个正方形的问题。

Lindemann 的证明否决了这一猜想 进入 20 世纪 50 年代，凭借着计算机技术的发展，π先是被计算到了上千位，随后又被计算到了上百万位一些新的技术被引进到π的计算中，比如在 1965 年的时候人们发现傅立叶变换（FFT）可以用来完成高精度的乘法运算，而运算速度却比传统的计算方案快得 多 尽管有了这样一些显著的进步， 在 20 世纪 70 年代之前所有计算机都是采用古典公式来计算π的，通常是 Machin 公式的变种到了 1976 年，Salamin 和 Richard Brent 独立发现了一个新的计算π的算法这个算法收敛到π的速度比任何古典公式都要快得多，该算法只需二十五次迭代就足以将π计算到超过四千五百万十进制位的精度 尽管靠着计算机的帮助，人类似乎已经把π的位数计算到了不可逾越的精度但是以上所有这些算法都有一个共同的特点：为了计算π的第 d 位，将不得不首先计算出 d 位之前的所有位换句话说，没有任何捷径可以让计算机能够直接计算出π的第 d 位因为这个缘故，所有这些已知的算法都表现出需要“吞噬大量内存”的特点，典型地，对内存的需 求量与要计算的总的位数成线性关系 因此，当一个能够直接计算出π的第 d 位的全新公式被发现的时候，带给人们的绝不是一个小小的惊喜。

这个公式，被称为 Bailey-Borwein-Plouffe 公式，简称为 BBP 公式，相 应的算法，被称为 BBP 算法 2 BBP 公式的证明公式的证明 BBP 公式是在 1996 年发现的[1]，该算法的正式论文发表于 1997 年[3]，在这篇文章中，提出了一种计算某些超越数（包括π，2π，log(2)，2(2)log，log(9 10)等）小数点后第 d 位数字的算法这些算法具有如下特点： l 算法计算出的小数点后第 d 位的数字不是以 10 为基数（radix）的，而是以 2 为 基数或者以 16 为基数，也就是实际计算出的是以二进制表示或者十六进制表示Generated by Foxit PDF Creator © Foxit Software For evaluation only.圆周率 pi 的 BBP 计算公式之详细证明 Rockins Chen UESTC, Chengdu 4/6/2009 3的超越数的第 d 位后的数字具体以何种基数表示取决于要计算的超越数的性质，如计算π时就是以 16 为基数 l 算法不需要多精度浮点算术运算的支持，在支持 IEEE 浮点运算的通用计算机上 即可进行计算 l 算法只需要非常少的内存 l 算法的时间复杂度与要计算的位 d 成线性关系， 这使得算法能够在普通的工作站上通过几个小时的计算求出log(2)或π的十亿个二进制位。

在文献[3]中，作者展示了将π计算到十六进制小数点后一百亿位的结果 计算π的 BBP 公式是： 014211()8184858616k kkkkkπ∞==−−−++++∑ (4) 此公式是计算机推导出的但要证明此公式并不难 【证明证明】首先注意到对任何8p UESTC, Chengdu 4/6/2009 4()()() ()()()()()() ( )( )( )14301220112200111222000212111 0001221616 244 1616 222448 222 2221224222112ln22ln224tan12ln12ln22ln24tan12ln2ln 22ln2lydyyyy ydyyyyyydydyyyy yydydydyyyyyyyyyii−−− −+− −=−−+−=−−−+ −=−+−−+−+=−−−++−=−−−+−−=−−+∫∫∫∫∫∫∫n2ππ+=证毕■ 3 BBP 公式的算法实现公式的算法实现 利用公式(4)，可以得到一个计算十六进制π的算法，算法能够在不计算小数点后前 d 位的情况下，直接计算出小数点后从 d 位开始的一串数字尽管该算法不是目前求π的效 率最高的算法， 但该算法最大的优势在于： 它能够直接计算出从小数点后某一位开始的一串 数字，而不依赖于 d 位之前的计算结果，这使它可以独立地验证π中某一位之后的几位数 字。

对于十六进制表示的π，将小数点向右移动一位等价于对π乘了 16这里为了表述方便，先引入一个算子{ }•，表示一个实数的小数部分因此，为了计算十六进制π的第 d位之后的数字，先将小数点移动到 d 位之后，也就等价于对π乘以16d，假定0d >，那么上式可以写成： {}{}{} {} {}{}1456421616161616dddddSSSSπ=−−− (5) 其中 01 16 (8)jk kSkj∞==+∑(6) 注意到 Generated by Foxit PDF Creator © Foxit Software For evaluation only.圆周率 pi 的 BBP 计算公式之详细证明 Rockins Chen UESTC, Chengdu 4/6/2009 5{} 01011616 8816mod816 8816d kd kdd j kk ddkd kdkk dSkjkjkj kjkj−−∞== +−−∞== +=++++=+++∑∑∑∑(7) 公式(7)把{}16d jS分成两部分来计算，第一部分计算 k 从 0 到 d，第二部分计算从 d 到∞。

对第一部分，之所以在分子中加入模8kj+操作，是因为公式(7)中只关注商的小数部分，因此通过取模操作简化运算第一部分包含总共 d 项的连加和，其中的每一项都是一个不大于8kj+的整数（16mod8d kkj−+）除以8kj+的商，所以用通常的计算机浮点算术运算即可完成每一项的除法操作，并求连加和 对于公式(7)中的第二部分，根据计算机浮点算术运算单元的精度，通常计算开始的几 项就够了，因为通过观察可知，第二部分中的每一项都是一个分数除以一个整数的商，随着 k 的增大，计算出来的每一项的值将迅速逼近计算机浮点算术运算单元所能表示的最小精 度 简单分析可知， 公式(5)中{}16dπ的值是不可能小于 0 的， 但是由于在公式(7)的第一部分中引入了取模操作，所以按照公式(5)计算出来的{}16dπ可能出现小于 0 的情况如果计算出的{}16dπ小于 0，那么需要做一次补偿，即给{}16dπ加上一个整数，使得最终结果位于 0 到 1 之间这个用来补偿的整数是{}()16d jabsS 接下来面临的一个问题是： 公式(7)的第一个连续和中， 每一项的分子16mod8d kkj−+应当如何快速地求解。

通常，对于指数计算，采用的方法是将指数因式分解，然后将指数二元展开（binary expansion） 比如，可以这样求解173： ()()222172333=•利用此方法，可以仅通过 5 次乘法运算计算出173的结果，而如果采用常规方法的话， 需要 16 次乘法操作 在公式(7)中，指数运算之后紧接着就是一次取模运算通过简单修改一下求指数的算 法可以将指数运算与取模运算合并到一个算法中以下就是指数取模的算法过程： 要计算modnrbc=（其中 r，b，n，c 为正整数） ，首先令 t 为 2 的幂，即2xt =，且满 足122xxn+≤ UESTC, Chengdu 4/6/2009 6if 1t ≥ then 2modrrc←; goto A; endif 这里取模运算符（mod）定义为：[]mod:xyxx y y=−，其中[ ]•算子表示取整注意到上述指数取模算法完全是在不超过2c的正整数上执行，因此算法能够在不产生舍入误差 的情况下正确执行 以上就是关于求解π的 BBP 算法的完整说明 作为示例，这里以 d=1,000,000 为例根据上面的算法可以求得： {}16dπ=0.423429797567540358599812674109⋯⋯ 但这还不是最终结果， 因为 BBP 算法求得的是第 d 位之后的十六进制数字 （即 d+1, d+2, d+3, …位上的十六进制数字） 。

上述结果还需要转换为十六进制表示，转换过程很简单：不 断将计算结果乘上 16，然后将相乘之后的结果中的整数部分取出，剩下的小数部分再做为计算结果，重复上述转换过程比如，对 d=1,000,000 时{}16dπ的值，按此转换过程可得如下序列：6，12，6，5，14，5，2，12，11，4，5，9，3，5，0，0，5，0，12，4，11， 11， 1， ⋯， 表示成十六进制符号就是： 6C65E52CB459350050E4BB1⋯， 这就是π从 1,000,001 位开始的 24 位十六进制数字 由于 BBP 算法采用常规浮点算术运算单元进行计算（而不是采用多精度浮点运算） ，因 此计算出的 d 位之后的十六进制数字中有多少位正确依赖于执行除法与求和运算的浮点算 术单元的精度根据作者 David H. Bailey 在文献[2]中给出的结论，在采用 IEEE 64 位浮点算 术的计算机上， 对于合理范围内的 d， 典型地能够产生 9 位甚至更多十六进制位； 在采用 128 位浮点算术的计算机上，典型地可以产生 24 位或者更多十六进制位能够产生正确的结果的d的范围大致是这样： 对于采用IEEE 64位浮点算术的计算机， BBP算法在d小于71.18 10× 时产生正确的结果； 对于采用 128 位浮点算术的计算机， BBP 算法在 d 小于1410时产生正确 的结果。

另外，由于 BBP 算法所具有的一个独特优点：能够独立计算某一位之后的数字因此 这也提供了一种验证结果是否正确的方法比如，如果已经得到了 d 位之后 n 位上的数字， 那么，可以再计算 d-1 位之后 n 位上的数字，以及 d+1 位之后 n 位上的数字，由于三次计算 的结果中有一些位是重叠的， 因此通过比较这些重叠的位上的数字是否相同就可以判断计算 结果是否正确 参考文献参考文献 [1] David H. Bailey, Jonathan M. Borwein, Peter B. Borwein, et al. The Quest for Pi. Mathematical Intelligencer. 1997, vol. 19 no. 1:50-57 [2] David H. Bailey. The BBP Algorithm for Pi (manuscrip。