圆周率pi的BBP计算公式之详细证明
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
π 尽管靠着计算机的帮助，人类似乎已经把 的位数计算到了不可逾越的精度。但是以 π 上所有这些算法都有一个共同的特点：为了计算 的第 d 位，将不得不首先计算出 d 位之 π 前的所有位。换句话说，没有任何捷径可以让计算机能够直接计算出 的第 d 位。因为这
个缘故，所有这些已知的算法都表现出需要“吞噬大量内存”的特点，典型地，对内存的需 求量与要计算的总的位数成线性关系。
Rockins Chen <ybc2084@gmail.com> UESTC, Chengdu 4/6/2009
Generated by Foxit PDF Creator © Foxit Software http://www.foxitsoftware.com For evaluation only.
{ } 16 分中引入了取模操作，所以按照公式(5)计算出来的 dπ 可能出现小于 0 的情况。如果计
{ } { } 16 16 算出的 dπ 小于 0，那么需要做一次补偿，即给 dπ 加上一个整数，使得最终结果位
({ }) 16 于 0 到 1 之间。这个用来补偿的整数是 abs
d
Sj
。
(1)
这个公式很容易得到。根据微积分公式和级数展开公式，可以将 tan−1 x 展开成下面的级 数：
∫ ∫ tan−1 x =
x dt 0 1+t2
=
x (1 − t 2 + t 4 − t6 + L)dt
0
(2)
= x − x3 3 + x5 5 − x7 7 + x9 9 −L
代入 x = 1 即可得到 Gregory-Leibniz 公式。Gregory-Leibniz 公式简洁优美，但是可惜的
Generated by Foxit PDF Creator © Foxit Software http://www.foxitsoftware.com For evaluation only.
圆周率 pi 的 BBP 计算公式之详细证明
5
{16 } ∑ ∑ d Sj
=
d k =0
Lindemann 的证明彻底地解决了源自古希腊的“化圆为方”问题，即能否仅用直尺和圆规将 一个圆变为一个正方形的问题。Lindemann 的证明否决了这一猜想。
π 进入 20 世纪 50 年代，凭借着计算机技术的发展， 先是被计算到了上千位，随后又 π 被计算到了上百万位。一些新的技术被引进到 的计算中，比如在 1965 年的时候人们发现
1 0
2
x1−1 1− x8
dx
−
2
×
24
2
1 0
2
x4−1 1− x8
dx
−
25
2
1 0
2
x5−1 1− x8
dx
−
26
2
1 0
x 2 6−1 1 − x8 dx
∫=
1 0
24
2 − 8x3 − 4 2x4 − 8x5
1− x8
dx
在上式中代入 y = 2x ，得
Rockins Chen <ybc2084@gmail.com> UESTC, Chengdu 4/6/2009
(3)
π 1873 年，Shanks 采用 Machin 的方案将 计算到了小数点后 707 位，不过可惜的是，
Rockins Chen <ybc2084@gmail.com> UESTC, Chengdu 4/6/2009
Generated by Foxit PDF Creator © Foxit Software http://www.foxitsoftware.com For evaluation only.
π 是，这个公式中的级数收敛得非常慢，以至于仅需要计算 的小数点后两位有效数字时，
就需要计算级数中的前面数百项。
随后，一个和 Gregory-Leibniz 公式类似但是收敛速度快得多的公式被提了出来，这个 公式要归功于 Machin：
π 4 = 4 tan−1(1 5) − tan−1(1 239)
π π2
log(2)
log
2
(
2)
log(9 10)
提出了一种计算某些超越数（包括 ， ，
，
，
等）小数点后
第 d 位数字的算法。这些算法具有如下特点：
l 算法计算出的小数点后第 d 位的数字不是以 10 为基数（radix）的，而是以 2 为 基数或者以 16 为基数，也就是实际计算出的是以二进制表示或者十六进制表示
∫ ∫ ∑ ∑ ∑ ( ) 1 0
2
x p−1 1 − x8
dx
=
1 0
2 ∞ x p−1+8k dx = ∞
1 x p+8k
k =0
k=0 p + 8k
1 0
2
=
1 2p 2
∞
k =0 16k
1 8k + p
于是有
∑ ( ) ∫ ∞
k =0 16k
1 8k + p
= 2p 2 1 0
x 2 p−1 1 − x8 dx
∞。对第一部分，之所以在分子中加入模 8k + j 操作，是因为公式(7)中只关注商的小数部 分，因此通过取模操作简化运算。第一部分包含总共 d 项的连加和，其中的每一项都是一个
不大于 8k + j 的整数（16d−k mod8k + j ）除以 8k + j 的商，所以用通常的计算机浮点算术
运算即可完成每一项的除法操作，并求连加和。
Generated by Foxit PDF Creator © Foxit Software http://www.foxitsoftware.com For evaluation only.
圆周率 pi 的 BBP 计算公式之详细证明
1
圆周率 pi 的 BBP 计算公式之详细证明
Rockins Chen <ybc2084@gmail.com> UESTC, Chengdu 4/6/2009
16d −k 8k + j
+
∞ k =d
+1
16d −k 8k + j
∑ ∑ =
k
d =0
16d
−
k
mod 8k +
8k j
+
j
+
k
∞ =d
+1
16 8k
d −k
+j
(7)
{ } 16 公式(7)把
d
Sj
分成两部分来计算，第一部分计算 k 从 0 到 d，第二部分计算从 d 到
圆周率 pi 的 BBP 计算公式之详细证明
3
的超越数的第 d 位后的数字。具体以何种基数表示取决于要计算的超越数的性
π 质，如计算 时就是以 16 为基数
l 算法不需要多精度浮点算术运算的支持，在支持 IEEE 浮点运算的通用计算机上 即可进行计算
l 算法只需要非常少的内存
l 算法的时间复杂度与要计算的位 d 成线性关系，这使得算法能够在普通的工作站
傅立叶变换（FFT）可以用来完成高精度的乘法运算，而运算速度却比传统的计算方案快得 多。
尽管有了这样一些显著的进步，在 20 世纪 70 年代之前所有计算机都是采用古典公式来
π 计算 的，通常是 Machin 公式的变种。到了 1976 年，Salamin 和 Richard Brent 独立发现了 π π 一个新的计算 的算法。这个算法收敛到 的速度比任何古典公式都要快得多，该算法只 π 需二十五次迭代就足以将 计算到超过四千五百万十进制位的精度。
π 因此，当一个能够直接计算出 的第 d 位的全新公式被发现的时候，带给人们的绝不
是一个小小的惊喜。这个公式，被称为 Bailey-Borwein-Plouffe 公式，简称为 BBP 公式，相 应的算法，被称为 BBP 算法。
2 BBP 公式的证明
BBP 公式是在 1996 年发现的[1]，该算法的正式论文发表于 1997 年[3]，在这篇文章中，
，换句话说就是 3.1408...
<π
<
3.1428... 。公元
π 5 世纪，中国古代数学家祖冲之采用阿基米德方法的一种变体将 求到了小数点后 7 位。
π 到了 17 世纪，微积分被牛顿和莱布尼兹发现之后，大量计算 的公式被发现了。著名
的 Gregory-Leibniz 公式就是其中一个：
π 4 = 1 −1 3 +1 5 −1 7 +1 9 −1 11+ L
π log(2)
上通过几个小时的计算求出
或
的十亿个二进制位。在文献[3]中，作者展
π 示了将 计算到十六进制小数点后一百亿位的结果
π 计算 的 BBP 公式是：
∞
∑ π =
1( 4
−
2
−
1
−
1
)
16 k=0 k 8k + 1 8k + 4 8k + 5 8k + 6
(4)
此公式是计算机推导出的。但要证明此公式并不难。 【证明】首先注意到对任何 p < 8 ，有下式成立：
圆周率 pi 的 BBP 计算公式之详细证明
2
后来发现 Shanks 的计算结果中第 527 位以后的数字是错误的。但这仍然是一个很大的进步。
π π 在不断追求计算更高精度 的过程中，一些和 有关的重要问题被解决了。17 世纪晚 π π 期，Lambert 和 Legendre 证明了 是无理数。1882 年，Lindemann 证明了 是超越数。
−
2 ln
y2 − 2y + 2
1 0
+4 tan−1
(
y
− 1)
1 0
= 2ln (−1) − 2ln (−2) + 2ln 2 − 4 tan−1 (−1)
= 2ln (i2 ) − 2ln (2) − 2ln (i2 ) + 2ln 2 + π
=π
证毕。■
3 BBP 公式的算法实现
π 利用公式(4)，可以得到一个计算十六进制 的算法，算法能够在不计算小数点后前 d π 位的情况下，直接计算出小数点后从 d 位开始的一串数字。尽管该算法不是目前求 的效
Generated by Foxit PDF Creator © Foxit Software http://www.foxitsoftware.com For evaluation only.
圆周率 pi 的 BBP 计算公式之详细证明
4
∫1 16 y −16
0
y4
−
2y3
+
4y
−
dy 4
( )( ) ∫1
=
16 y −16
dy
0 y2 − 2 y2 − 2y + 2
∫ ∫ =
1 0
4y y2 −
dy 2
−
1 4y −8
0
y2
−
2
y
+
dy 2
∫ ∫ ∫ =
2
1 0
2y
y2
−
dy 2
−2
1 0
y
2
2y −2 −2y +
dy 2
+
4
1
0 (y
1
− 1)2
dy +1
( ) ( ) = 2ln
y2 − 2
1 0
π 16 位之后的数字，先将小数点移动到 d 位之后，也就等价于对 乘以 d ，假定 d > 0 ，那
么上式可以写成：
{ } { } { } { } { } { } 16dπ
=
4 16dS1
−
2
16
d
S4
−
wk.baidu.com16
d
S5
− 16dS6
(5)
其中
∑ S j
=
∞
1
k=0 16k (8k
+
j)
(6)
注意到
Rockins Chen <ybc2084@gmail.com> UESTC, Chengdu 4/6/2009
对于公式(7)中的第二部分，根据计算机浮点算术运算单元的精度，通常计算开始的几 项就够了，因为通过观察可知，第二部分中的每一项都是一个分数除以一个整数的商，随着 k 的增大，计算出来的每一项的值将迅速逼近计算机浮点算术运算单元所能表示的最小精 度。
{ } 16 简单分析可知，公式(5)中 dπ 的值是不可能小于 0 的，但是由于在公式(7)的第一部
摘要：本文给出了用于计算圆周率
pi
的
BBP
[1,2,3]
公式的详细证明。主要是对文献 的翻译和
解释。
关键字：圆周率 pi BBP 公式
1 求解 pi 的历史
π 第一个尝试严格计算 的人是公元前 250 年的阿基米德，他采用内接多边形和外切多
π 边形逼近的办法确定出
310 < π
的范围是 71
< 31 7
率最高的算法，但该算法最大的优势在于：它能够直接计算出从小数点后某一位开始的一串
π 数字，而不依赖于 d 位之前的计算结果，这使它可以独立地验证 中某一位之后的几位数
字。
π π 对于十六进制表示的 ，将小数点向右移动一位等价于对 乘了 16。这里为了表述方 π 便，先引入一个算子{•} ，表示一个实数的小数部分。因此，为了计算十六进制 的第 d
因此公式(4)(4)等号右边可以写成
∞
∑16 k =0
1
k
4 8k +
1
−
2 8k +
4
−
1 8k +
5
−
1 8k +
6
=
∞
4
∑k 16 k =0
1
(8k
+ 1)
−
∞
2
∑k 16 k =0
1
(8k
+
4)
−
∞
∑k 16 k =0
1
(8k
+
5)
−
∞
∑k 16 k =0
1
(8k
+
6)
∫ ∫ ∫ ∫ = 4 × 21 2
接下来面临的一个问题是：公式(7)的第一个连续和中，每一项的分子16d−k mod8k + j 应
当如何快速地求解。通常，对于指数计算，采用的方法是将指数因式分解，然后将指数二元 展开（binary expansion）。比如，可以这样求解 317 ：
个缘故，所有这些已知的算法都表现出需要“吞噬大量内存”的特点，典型地，对内存的需 求量与要计算的总的位数成线性关系。
Rockins Chen <ybc2084@gmail.com> UESTC, Chengdu 4/6/2009
Generated by Foxit PDF Creator © Foxit Software http://www.foxitsoftware.com For evaluation only.
{ } 16 分中引入了取模操作，所以按照公式(5)计算出来的 dπ 可能出现小于 0 的情况。如果计
{ } { } 16 16 算出的 dπ 小于 0，那么需要做一次补偿，即给 dπ 加上一个整数，使得最终结果位
({ }) 16 于 0 到 1 之间。这个用来补偿的整数是 abs
d
Sj
。
(1)
这个公式很容易得到。根据微积分公式和级数展开公式，可以将 tan−1 x 展开成下面的级 数：
∫ ∫ tan−1 x =
x dt 0 1+t2
=
x (1 − t 2 + t 4 − t6 + L)dt
0
(2)
= x − x3 3 + x5 5 − x7 7 + x9 9 −L
代入 x = 1 即可得到 Gregory-Leibniz 公式。Gregory-Leibniz 公式简洁优美，但是可惜的
Generated by Foxit PDF Creator © Foxit Software http://www.foxitsoftware.com For evaluation only.
圆周率 pi 的 BBP 计算公式之详细证明
5
{16 } ∑ ∑ d Sj
=
d k =0
Lindemann 的证明彻底地解决了源自古希腊的“化圆为方”问题，即能否仅用直尺和圆规将 一个圆变为一个正方形的问题。Lindemann 的证明否决了这一猜想。
π 进入 20 世纪 50 年代，凭借着计算机技术的发展， 先是被计算到了上千位，随后又 π 被计算到了上百万位。一些新的技术被引进到 的计算中，比如在 1965 年的时候人们发现
1 0
2
x1−1 1− x8
dx
−
2
×
24
2
1 0
2
x4−1 1− x8
dx
−
25
2
1 0
2
x5−1 1− x8
dx
−
26
2
1 0
x 2 6−1 1 − x8 dx
∫=
1 0
24
2 − 8x3 − 4 2x4 − 8x5
1− x8
dx
在上式中代入 y = 2x ，得
Rockins Chen <ybc2084@gmail.com> UESTC, Chengdu 4/6/2009
(3)
π 1873 年，Shanks 采用 Machin 的方案将 计算到了小数点后 707 位，不过可惜的是，
Rockins Chen <ybc2084@gmail.com> UESTC, Chengdu 4/6/2009
Generated by Foxit PDF Creator © Foxit Software http://www.foxitsoftware.com For evaluation only.
π 是，这个公式中的级数收敛得非常慢，以至于仅需要计算 的小数点后两位有效数字时，
就需要计算级数中的前面数百项。
随后，一个和 Gregory-Leibniz 公式类似但是收敛速度快得多的公式被提了出来，这个 公式要归功于 Machin：
π 4 = 4 tan−1(1 5) − tan−1(1 239)
π π2
log(2)
log
2
(
2)
log(9 10)
提出了一种计算某些超越数（包括 ， ，
，
，
等）小数点后
第 d 位数字的算法。这些算法具有如下特点：
l 算法计算出的小数点后第 d 位的数字不是以 10 为基数（radix）的，而是以 2 为 基数或者以 16 为基数，也就是实际计算出的是以二进制表示或者十六进制表示
∫ ∫ ∑ ∑ ∑ ( ) 1 0
2
x p−1 1 − x8
dx
=
1 0
2 ∞ x p−1+8k dx = ∞
1 x p+8k
k =0
k=0 p + 8k
1 0
2
=
1 2p 2
∞
k =0 16k
1 8k + p
于是有
∑ ( ) ∫ ∞
k =0 16k
1 8k + p
= 2p 2 1 0
x 2 p−1 1 − x8 dx
∞。对第一部分，之所以在分子中加入模 8k + j 操作，是因为公式(7)中只关注商的小数部 分，因此通过取模操作简化运算。第一部分包含总共 d 项的连加和，其中的每一项都是一个
不大于 8k + j 的整数（16d−k mod8k + j ）除以 8k + j 的商，所以用通常的计算机浮点算术
运算即可完成每一项的除法操作，并求连加和。
Generated by Foxit PDF Creator © Foxit Software http://www.foxitsoftware.com For evaluation only.
圆周率 pi 的 BBP 计算公式之详细证明
1
圆周率 pi 的 BBP 计算公式之详细证明
Rockins Chen <ybc2084@gmail.com> UESTC, Chengdu 4/6/2009
16d −k 8k + j
+
∞ k =d
+1
16d −k 8k + j
∑ ∑ =
k
d =0
16d
−
k
mod 8k +
8k j
+
j
+
k
∞ =d
+1
16 8k
d −k
+j
(7)
{ } 16 公式(7)把
d
Sj
分成两部分来计算，第一部分计算 k 从 0 到 d，第二部分计算从 d 到
圆周率 pi 的 BBP 计算公式之详细证明
3
的超越数的第 d 位后的数字。具体以何种基数表示取决于要计算的超越数的性
π 质，如计算 时就是以 16 为基数
l 算法不需要多精度浮点算术运算的支持，在支持 IEEE 浮点运算的通用计算机上 即可进行计算
l 算法只需要非常少的内存
l 算法的时间复杂度与要计算的位 d 成线性关系，这使得算法能够在普通的工作站
傅立叶变换（FFT）可以用来完成高精度的乘法运算，而运算速度却比传统的计算方案快得 多。
尽管有了这样一些显著的进步，在 20 世纪 70 年代之前所有计算机都是采用古典公式来
π 计算 的，通常是 Machin 公式的变种。到了 1976 年，Salamin 和 Richard Brent 独立发现了 π π 一个新的计算 的算法。这个算法收敛到 的速度比任何古典公式都要快得多，该算法只 π 需二十五次迭代就足以将 计算到超过四千五百万十进制位的精度。
π 因此，当一个能够直接计算出 的第 d 位的全新公式被发现的时候，带给人们的绝不
是一个小小的惊喜。这个公式，被称为 Bailey-Borwein-Plouffe 公式，简称为 BBP 公式，相 应的算法，被称为 BBP 算法。
2 BBP 公式的证明
BBP 公式是在 1996 年发现的[1]，该算法的正式论文发表于 1997 年[3]，在这篇文章中，
，换句话说就是 3.1408...
<π
<
3.1428... 。公元
π 5 世纪，中国古代数学家祖冲之采用阿基米德方法的一种变体将 求到了小数点后 7 位。
π 到了 17 世纪，微积分被牛顿和莱布尼兹发现之后，大量计算 的公式被发现了。著名
的 Gregory-Leibniz 公式就是其中一个：
π 4 = 1 −1 3 +1 5 −1 7 +1 9 −1 11+ L
π log(2)
上通过几个小时的计算求出
或
的十亿个二进制位。在文献[3]中，作者展
π 示了将 计算到十六进制小数点后一百亿位的结果
π 计算 的 BBP 公式是：
∞
∑ π =
1( 4
−
2
−
1
−
1
)
16 k=0 k 8k + 1 8k + 4 8k + 5 8k + 6
(4)
此公式是计算机推导出的。但要证明此公式并不难。 【证明】首先注意到对任何 p < 8 ，有下式成立：
圆周率 pi 的 BBP 计算公式之详细证明
2
后来发现 Shanks 的计算结果中第 527 位以后的数字是错误的。但这仍然是一个很大的进步。
π π 在不断追求计算更高精度 的过程中，一些和 有关的重要问题被解决了。17 世纪晚 π π 期，Lambert 和 Legendre 证明了 是无理数。1882 年，Lindemann 证明了 是超越数。
−
2 ln
y2 − 2y + 2
1 0
+4 tan−1
(
y
− 1)
1 0
= 2ln (−1) − 2ln (−2) + 2ln 2 − 4 tan−1 (−1)
= 2ln (i2 ) − 2ln (2) − 2ln (i2 ) + 2ln 2 + π
=π
证毕。■
3 BBP 公式的算法实现
π 利用公式(4)，可以得到一个计算十六进制 的算法，算法能够在不计算小数点后前 d π 位的情况下，直接计算出小数点后从 d 位开始的一串数字。尽管该算法不是目前求 的效
Generated by Foxit PDF Creator © Foxit Software http://www.foxitsoftware.com For evaluation only.
圆周率 pi 的 BBP 计算公式之详细证明
4
∫1 16 y −16
0
y4
−
2y3
+
4y
−
dy 4
( )( ) ∫1
=
16 y −16
dy
0 y2 − 2 y2 − 2y + 2
∫ ∫ =
1 0
4y y2 −
dy 2
−
1 4y −8
0
y2
−
2
y
+
dy 2
∫ ∫ ∫ =
2
1 0
2y
y2
−
dy 2
−2
1 0
y
2
2y −2 −2y +
dy 2
+
4
1
0 (y
1
− 1)2
dy +1
( ) ( ) = 2ln
y2 − 2
1 0
π 16 位之后的数字，先将小数点移动到 d 位之后，也就等价于对 乘以 d ，假定 d > 0 ，那
么上式可以写成：
{ } { } { } { } { } { } 16dπ
=
4 16dS1
−
2
16
d
S4
−
wk.baidu.com16
d
S5
− 16dS6
(5)
其中
∑ S j
=
∞
1
k=0 16k (8k
+
j)
(6)
注意到
Rockins Chen <ybc2084@gmail.com> UESTC, Chengdu 4/6/2009
对于公式(7)中的第二部分，根据计算机浮点算术运算单元的精度，通常计算开始的几 项就够了，因为通过观察可知，第二部分中的每一项都是一个分数除以一个整数的商，随着 k 的增大，计算出来的每一项的值将迅速逼近计算机浮点算术运算单元所能表示的最小精 度。
{ } 16 简单分析可知，公式(5)中 dπ 的值是不可能小于 0 的，但是由于在公式(7)的第一部
摘要：本文给出了用于计算圆周率
pi
的
BBP
[1,2,3]
公式的详细证明。主要是对文献 的翻译和
解释。
关键字：圆周率 pi BBP 公式
1 求解 pi 的历史
π 第一个尝试严格计算 的人是公元前 250 年的阿基米德，他采用内接多边形和外切多
π 边形逼近的办法确定出
310 < π
的范围是 71
< 31 7
率最高的算法，但该算法最大的优势在于：它能够直接计算出从小数点后某一位开始的一串
π 数字，而不依赖于 d 位之前的计算结果，这使它可以独立地验证 中某一位之后的几位数
字。
π π 对于十六进制表示的 ，将小数点向右移动一位等价于对 乘了 16。这里为了表述方 π 便，先引入一个算子{•} ，表示一个实数的小数部分。因此，为了计算十六进制 的第 d
因此公式(4)(4)等号右边可以写成
∞
∑16 k =0
1
k
4 8k +
1
−
2 8k +
4
−
1 8k +
5
−
1 8k +
6
=
∞
4
∑k 16 k =0
1
(8k
+ 1)
−
∞
2
∑k 16 k =0
1
(8k
+
4)
−
∞
∑k 16 k =0
1
(8k
+
5)
−
∞
∑k 16 k =0
1
(8k
+
6)
∫ ∫ ∫ ∫ = 4 × 21 2
接下来面临的一个问题是：公式(7)的第一个连续和中，每一项的分子16d−k mod8k + j 应
当如何快速地求解。通常，对于指数计算，采用的方法是将指数因式分解，然后将指数二元 展开（binary expansion）。比如，可以这样求解 317 ：