圆周率pi的BBP计算公式之详细证明
Generated by Foxit PDF Creator Foxit Softwarehttp://www.foxitsoftware.com For evaluation only.
圆周率
pi 的BBP 计算公式之详细证明
1
圆周率pi 的BBP 计算公式之详细证明
Rockins Chen UESTC, Chengdu 4/6/2009
摘要：本文给出了用于计算圆周率pi 的BBP 公式的详细证明。主要是对文献[1,2,3]的翻译和
解释。
关键字：圆周率pi BBP 公式
1 求解pi 的历史
第一个尝试严格计算边形逼近的办法确定出
π的人是公元前250年的阿基米德，他采用内接多边形和外切多
π
π的范围是
5世纪，中国古代数学家祖冲之采用阿基米德方法的一种变体将
π求到了小数点后7位。 π的公式被发现了。著名
(1)
到了17世纪，微积分被牛顿和莱布尼兹发现之后，大量计算的Gregory-Leibniz 公式就是其中一个：
=1−+−+−+L
这个公式很容易得到。根据微积分公式和级数展开公式，可以将tan −1x 展开成下面的级
数：
x dt
tan x =∫=∫(1−t 2+t 4−t 6+L ) dt 00
3579
=x −x +x −x +x −L −1
x
(2)
代入x =1即可得到Gregory-Leibniz 公式。Gregory-Leibniz 公式简洁优美，但是可惜的是，这个公式中的级数收敛得非常慢，以至于仅需要计算就需要计算级数中的前面数百项。
随后，一个和Gregory-Leibniz 公式类似但是收敛速度快得多的公式被提了出来，这个公式要归功于Machin ：
π的小数点后两位有效数字时，
π=4tan −1−tan −1
1873年，Shanks 采用Machin 的方案将
(3)
π计算到了小数点后707位，不过可惜的是，
Generated by Foxit PDF Creator Foxit Softwarehttp://www.foxitsoftware.com For evaluation only.
圆周率
pi 的BBP 计算公式之详细证明
2
后来发现Shanks 的计算结果中第527位以后的数字是错误的。但这仍然是一个很大的进步。
在不断追求计算更高精度
π的过程中，一些和π有关的重要问题被解决了。17世纪晚π是无理数。1882年，Lindemann 证明了π是超越数。
期，Lambert 和Legendre 证明了
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)，等）小数点后
提出了一种计算某些超越数（包括π，π，log(2)，
2
第d 位数字的算法。这些算法具有如下特点：
l
算法计算出的小数点后第d 位的数字不是以10为基数（radix ）的，而是以2为基数或者以16为基数，也就是实际计算出的是以二进制表示或者十六进制表示
Generated by Foxit PDF Creator Foxit Softwarehttp://www.foxitsoftware.com For evaluation only.
圆周率
pi 的BBP 计算公式之详细证明
3
的超越数的第d 位后的数字。具体以何种基数表示取决于要计算的超越数的性质，如计算l l l
π时就是以16为基数
算法不需要多精度浮点算术运算的支持，在支持IEEE 浮点运算的通用计算机上即可进行计算
算法只需要非常少的内存
算法的时间复杂度与要计算的位d 成线性关系，这使得算法能够在普通的工作站上通过几个小时的计算求出log(2)或示了将
π的十亿个二进制位。在文献[3]中，作者展
π计算到十六进制小数点后一百亿位的结果
1
4211
−−− 8k +18k +48k +58k +6
计算
π的BBP 公式是：
π=∑
k =0∞
16
(4)
此公式是计算机推导出的。但要证明此公式并不难。 【证明
】首先注意到对任何p
∫
=∫0∞
x
k =0
p −1+8k
dx =∑
k =0
∞
1
x p +=
11
k =0
∞
于是有
1p 2=∫k =0168k +p ∞
因此公式(4)(4)等号右边可以写成
14211−−− k =0
16=4k =0∞
∞
1
168k +112
2∑
k =0
∞
1
168k +4−∑
k =0
∞
1
16
8k +5−∑
k =0
∞
1
16
8k +64252
−2×2∫−2∫=4×2∫=∫62
−2∫
在上式中代入y =，得
Generated by Foxit PDF Creator Foxit Softwarehttp://www.foxitsoftware.com For evaluation only.
圆周率
pi 的BBP 计算公式之详细证明
4
16y 16
∫0y −2y +4y −4116y −16=∫0
y −2y −2y +21
14y 4y 8
=∫−∫0y −20y −2y +212y 112y 21=2∫−2∫+4∫ 0y −20y −2y +20
y −1+11
2−1=2ln (y 2−2)10−2ln (y −2y +2)0+4tan (y −1)0
=2ln (−1)−2ln (−2)+2ln 2−4tan −1(−1)=2ln (i 2)−2ln (2)−2ln (i 2)+2l n 2+π=π
证毕。■
3 BBP公式的算法实现
利用公式(4)，可以得到一个计算十六进制的算法，算法能够在不计算小数点后前d 位的情况下，直接计算出小数点后从d 位开始的一串数字。尽管该算法不是目前求的效率最高的算法，但该算法最大的优势在于：它能够直接计算出从小数点后某一位开始的一串数字，而不依赖于d 位之前的计算结果，这使它可以独立地验证中某一位之后的几位数字。
π
π
π
π，将小数点向右移动一位等价于对π乘了16。这里为了表述方
便，先引入一个算子{•}，表示一个实数的小数部分。因此，为了计算十六进制π的第d
对于十六进制表示的
位之后的数字，先将小数点移动到d 位之后，也就等价于对么上式可以写成：
π乘以16
6
d
，假定d >0，那
{16π}={4{16S }−2{16S }−{16S }−{16S }}
d
d
d
d
d
1
4
5
(5)
其中
S j =∑
注意到
1
16(8k +j ) k =0
∞
(6)
Generated by Foxit PDF Creator Foxit Softwarehttp://www.foxitsoftware.com For evaluation only.
圆周率
pi 的BBP 计算公式之详细证明
5
{16S }
d
j
d −k ∞ 16d −k d 16 = +∑ k =08k +j k =d +18k +j
∞ 16d −k d 16d −k mod8k j = +∑8k +j 8k +j k =0k =d +1
(7)
公式(7)把
{16S }分成两部分来计算，第一部分计算k 从0到d ，第二部分计算从d 到
d
j
∞。对第一部分，之所以在分子中加入模8k +j 操作，是因为公式(7)中只关注商的小数部分，因此通过取模操作简化运算。第一部分包含总共d 项的连加和，其中的每一项都是一个不大于8k +j 的整数（16d −k mod8k +j ）除以8k +j 的商，所以用通常的计算机浮点算术运算即可完成每一项的除法操作，并求连加和。
对于公式(7)中的第二部分，根据计算机浮点算术运算单元的精度，通常计算开始的几
项就够了，因为通过观察可知，第二部分中的每一项都是一个分数除以一个整数的商，随着k 的增大，计算出来的每一项的值将迅速逼近计算机浮点算术运算单元所能表示的最小精度。
{16π}的值是不可能小于0的，但是由于在公式(7)的第一部
分中引入了取模操作，所以按照公式(5)计算出来的{16π}可能出现小于0的情况。如果计算出的{16π}小于0，那么需要做一次补偿，即给{16π}加上一个整数，使得最终结果位于0到1之间。这个用来补偿的整数是 abs ({16S }) 。
简单分析可知，公式(5)中
d
d
d
d
d
j
接下来面临的一个问题是：公式(7)的第一个连续和中，每一项的分子16d −k mod8k +j 应当如何快速地求解。通常，对于指数计算，采用的方法是将指数因式分解，然后将指数二元展开（binary expansion）。比如，可以这样求解317：
2
317= (32)
()
2
•3
2
利用此方法，可以仅通过5次乘法运算计算出317的结果，而如果采用常规方法的话，需要16次乘法操作。
在公式(7)中，指数运算之后紧接着就是一次取模运算。通过简单修改一下求指数的算法可以将指数运算与取模运算合并到一个算法中。以下就是指数取模的算法过程： 要计算r =b n mod c （其中r ，b ，n ，c 为正整数），首先令t 为2的幂，即t =2x ，且满足2x ≤n
A: if n ≥t then r ←br mod c ; n ←n −t ; endif
t ←
Generated by Foxit PDF Creator Foxit Softwarehttp://www.foxitsoftware.com For evaluation only.
圆周率
pi 的BBP 计算公式之详细证明
6
if t ≥1 then r ←r 2mod c ; goto A; endif
这里取模运算符（mod ）定义为：x mod y :=x −[]y ，其中[•]算子表示取整。注意到上述指数取模算法完全是在不超过c 2的正整数上执行，因此算法能够在不产生舍入误差的情况下正确执行。
以上就是关于求解
π的BBP 算法的完整说明。
作为示例，这里以d =1,000,000为例。根据上面的算法可以求得：
{16π}=0.[***********][1**********]9…
d
但这还不是最终结果，因为BBP 算法求得的是第d 位之后的十六进制数字（即d +1, d +2, d +3, …位上的十六进制数字）。上述结果还需要转换为十六进制表示，转换过程很简单：不断将计算结果乘上16，然后将相乘之后的结果中的整数部分取出，剩下的小数部分再做为计算结果，重复上述转换过程。比如，对d =1,000,000时
{16π}的值，按此转换过程可得
d
如下序列：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 小于1.18×107时产生正确的结果；对于采用128位浮点算术的计算机，BBP 算法在d 小于1014时产生正确的结果。
另外，由于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 (manuscript). available at http://crd.lbl.gov/~dhbailey/dhbpapers/bbp-alg.pdf
[3] David Bailey, Peter Borwein and Simon Plouffe. On The Rapid Computation Of Various Polylogarithmic Constants. Mathematics of Computation. 1997, 66(218):903-913