贝利-波尔温-普劳夫公式 （英語：Bailey–Borwein–Plouffe formula ）提供了一个计算圓周率π n 位二进制 数的spigot算法 spigot algorithm 求和 公式是在1995年由西蒙·普勞夫 提出的，并以公布这个公式的论文作者大衛·H·貝利 彼得·博温 [ 1] [ 2]
π
=
∑
k
=
0
∞
[
1
16
k
(
4
8
k
+
1
−
2
8
k
+
4
−
1
8
k
+
5
−
1
8
k
+
6
)
]
{\displaystyle \pi =\sum _{k=0}^{\infty }\left[{\frac {1}{16^{k}}}\left({\frac {4}{8k+1}}-{\frac {2}{8k+4}}-{\frac {1}{8k+5}}-{\frac {1}{8k+6}}\right)\right]}
这个公式的发现曾震惊学界。数百年来，求出π的第n 位小数而不求出它的前n -1位曾被认为是不可能的。
自从这个发现以来，发现了更多的无理数常数的类似公式，它们都有一个类似的形式：
α
=
∑
k
=
0
∞
[
1
b
k
p
(
k
)
q
(
k
)
]
{\displaystyle \alpha =\sum _{k=0}^{\infty }\left[{\frac {1}{b^{k}}}{\frac {p(k)}{q(k)}}\right]}
其中α是目标常数，p 和q 是整系数多项式 ，b ≥ 2是整数的数制 。
这种形式的公式被称为BBP式公式 （BBP-type formulas）[ 3] p ,q 和b 可组合出一些著名的常数。但至今尚未找出一种系统的算法来寻找合适的组合，而已知的公式多是通过实验数学 得出的。
编辑
一个已经得出很多解的BBP式的特例是：
P
(
s
,
b
,
m
,
A
)
=
∑
k
=
0
∞
[
1
b
k
∑
j
=
1
m
a
j
(
m
k
+
j
)
s
]
{\displaystyle P(s,b,m,A)=\sum _{k=0}^{\infty }\left[{\frac {1}{b^{k}}}\sum _{j=1}^{m}{\frac {a_{j}}{(mk+j)^{s}}}\right]}
其中s , b 和m 是整数，
A
=
(
a
1
,
a
2
,
…
,
a
m
)
{\displaystyle A=(a_{1},a_{2},\dots ,a_{m})}
编辑
在BBP公式发现之前，就已经有些最简单的此类公式广为所知。他们使用P函数的省略记法如下：
ln
9
10
=
−
1
10
−
1
200
−
1
3
000
−
1
40
000
−
1
500
000
−
⋯
=
−
∑
k
=
1
∞
1
10
k
⋅
k
=
−
1
10
∑
k
=
0
∞
[
1
10
k
(
1
k
+
1
)
]
=
−
1
10
P
(
1
,
10
,
1
,
(
1
)
)
{\displaystyle {\begin{aligned}\ln {\frac {9}{10}}&=-{\frac {1}{10}}-{\frac {1}{200}}-{\frac {1}{3\ 000}}-{\frac {1}{40\ 000}}-{\frac {1}{500\ 000}}-\cdots \\&=-\sum _{k=1}^{\infty }{\frac {1}{10^{k}\cdot k}}=-{\frac {1}{10}}\sum _{k=0}^{\infty }\left[{\frac {1}{10^{k}}}\left({\frac {1}{k+1}}\right)\right]\\&=-{\frac {1}{10}}P\left(1,10,1,(1)\right)\end{aligned}}}
ln
2
=
1
2
+
1
2
⋅
2
2
+
1
3
⋅
2
3
+
1
4
⋅
2
4
+
1
5
⋅
2
5
+
⋯
=
∑
k
=
1
∞
1
2
k
⋅
k
=
1
2
∑
k
=
0
∞
[
1
2
k
(
1
k
+
1
)
]
=
1
2
P
(
1
,
2
,
1
,
(
1
)
)
{\displaystyle {\begin{aligned}\ln 2&={\frac {1}{2}}+{\frac {1}{2\cdot 2^{2}}}+{\frac {1}{3\cdot 2^{3}}}+{\frac {1}{4\cdot 2^{4}}}+{\frac {1}{5\cdot 2^{5}}}+\cdots \\&=\sum _{k=1}^{\infty }{\frac {1}{2^{k}\cdot k}}={\frac {1}{2}}\sum _{k=0}^{\infty }\left[{\frac {1}{2^{k}}}\left({\frac {1}{k+1}}\right)\right]\\&={\frac {1}{2}}P\left(1,2,1,(1)\right)\end{aligned}}}
普勞夫也对这种形式的反正切函數 的級數展開感兴趣（P记法还可以泛化为b 不是整数的情形）:
arctan
1
b
=
1
b
−
1
b
3
3
+
1
b
5
5
−
1
b
7
7
+
1
b
9
9
+
⋯
=
∑
k
=
1
∞
[
1
b
k
sin
k
π
2
k
]
=
1
b
∑
k
=
0
∞
[
1
b
4
k
(
1
4
k
+
1
+
−
1
4
k
+
3
)
]
=
1
b
P
(
1
,
b
4
,
4
,
(
1
,
0
,
−
1
,
0
)
)
{\displaystyle {\begin{aligned}\arctan {\frac {1}{b}}&={\frac {1}{b}}-{\frac {1}{b^{3}3}}+{\frac {1}{b^{5}5}}-{\frac {1}{b^{7}7}}+{\frac {1}{b^{9}9}}+\cdots \\&=\sum _{k=1}^{\infty }\left[{\frac {1}{b^{k}}}{\frac {\sin {\frac {k\pi }{2}}}{k}}\right]={\frac {1}{b}}\sum _{k=0}^{\infty }\left[{\frac {1}{b^{4k}}}\left({\frac {1}{4k+1}}+{\frac {-1}{4k+3}}\right)\right]\\&={\frac {1}{b}}P\left(1,b^{4},4,(1,0,-1,0)\right)\end{aligned}}}
编辑
使用上述P函数，令s = 1 , m > 1 可得出已知的π的最简单公式。很多已知的公式是令b是2或3的幂，m是2的幂或者其他多因子的值，并令向量A等於零。在计算了一些类似的和之后，这类发现引发了使用计算机搜索类似线性组合的尝试。搜索程序为每个参数，s, b, m分别选择一个定义域，然后求和并检查值，并使用整数关系侦查算法 integer relation algorithm 赫拉曼·弗古森 Helaman Ferguson A 使得这些中间值可以加在一起得到一个著名的常数（A 也可能是零）。
编辑
原始的BBP π求和公式是在1995年由Plouffe使用PSLQ算法 [ 4] integer relation algorithm P 函数表示：
π
=
∑
k
=
0
∞
[
1
16
k
(
4
8
k
+
1
−
2
8
k
+
4
−
1
8
k
+
5
−
1
8
k
+
6
)
]
=
P
(
1
,
16
,
8
,
(
4
,
0
,
0
,
−
2
,
−
1
,
−
1
,
0
,
0
)
)
{\displaystyle {\begin{aligned}\pi &=\sum _{k=0}^{\infty }\left[{\frac {1}{16^{k}}}\left({\frac {4}{8k+1}}-{\frac {2}{8k+4}}-{\frac {1}{8k+5}}-{\frac {1}{8k+6}}\right)\right]\\&=P\left(1,16,8,(4,0,0,-2,-1,-1,0,0)\right)\end{aligned}}}
这个公式也可以使用下述两个多项式的商来表示：
π
=
∑
k
=
0
∞
[
1
16
k
(
120
k
2
+
151
k
+
47
512
k
4
+
1024
k
3
+
712
k
2
+
194
k
+
15
)
]
{\displaystyle \pi =\sum _{k=0}^{\infty }\left[{\frac {1}{16^{k}}}\left({\frac {120k^{2}+151k+47}{512k^{4}+1024k^{3}+712k^{2}+194k+15}}\right)\right]}
这个等式可以用一个较为简单的方式严格证明。[ 5]
编辑
这个公式给出一个抽取π的十六进制 位数值的算法。首先我们需要将这个公式化为：
π
=
4
∑
k
=
0
∞
1
(
16
k
)
(
8
k
+
1
)
−
2
∑
k
=
0
∞
1
(
16
k
)
(
8
k
+
4
)
−
∑
k
=
0
∞
1
(
16
k
)
(
8
k
+
5
)
−
∑
k
=
0
∞
1
(
16
k
)
(
8
k
+
6
)
{\displaystyle \pi =4\sum _{k=0}^{\infty }{\frac {1}{(16^{k})(8k+1)}}-2\sum _{k=0}^{\infty }{\frac {1}{(16^{k})(8k+4)}}-\sum _{k=0}^{\infty }{\frac {1}{(16^{k})(8k+5)}}-\sum _{k=0}^{\infty }{\frac {1}{(16^{k})(8k+6)}}}
在使用此公式实现一个spigot算法之前，还需做一些改动。我们想要找出π在十六进制 下的第n 位，所以首先我们需要将该求和序列拆为n之前和n之后两部分。对于前述公式的第一项而言，
∑
k
=
0
∞
1
(
16
k
)
(
8
k
+
1
)
=
∑
k
=
0
n
1
(
16
k
)
(
8
k
+
1
)
+
∑
k
=
n
+
1
∞
1
(
16
k
)
(
8
k
+
1
)
{\displaystyle \sum _{k=0}^{\infty }{\frac {1}{(16^{k})(8k+1)}}=\sum _{k=0}^{n}{\frac {1}{(16^{k})(8k+1)}}+\sum _{k=n+1}^{\infty }{\frac {1}{(16^{k})(8k+1)}}}
我们再将等式两边乘以16 n ，使得小数点恰好落在第n 位。
∑
k
=
0
∞
16
n
−
k
8
k
+
1
=
∑
k
=
0
n
16
n
−
k
8
k
+
1
+
∑
k
=
n
+
1
∞
16
n
−
k
8
k
+
1
{\displaystyle \sum _{k=0}^{\infty }{\frac {16^{n-k}}{8k+1}}=\sum _{k=0}^{n}{\frac {16^{n-k}}{8k+1}}+\sum _{k=n+1}^{\infty }{\frac {16^{n-k}}{8k+1}}}
由于我们关心的是小数部分，我们观察这个式子的两项，会发现仅有第一项能够出现整数部分，而第二项不会出现整数部分。因为第二项中k > n ,第二项中的分子一定小于分母。由此我们需要一个使用一种技巧来从第一项中除去整数。那就是要mod 8k + 1。于是原式第一项的小数部分的公式就是：
∑
k
=
0
n
16
n
−
k
mod
8
k
+
1
8
k
+
1
+
∑
k
=
n
+
1
∞
16
n
−
k
8
k
+
1
{\displaystyle \sum _{k=0}^{n}{\frac {16^{n-k}\mod 8k+1}{8k+1}}+\sum _{k=n+1}^{\infty }{\frac {16^{n-k}}{8k+1}}}
注意模 运算将保证只有小数部分的和会被保留下来。想要快速有效的计算16 n - k mod 8k + 1 ,可使用模幂運算 （Modular exponentiation
对公式的其余项使用相同的处理办法，并根据原式求出最后的结果。
4
Σ
1
−
2
Σ
2
−
Σ
3
−
Σ
4
.
{\displaystyle 4\Sigma _{1}-2\Sigma _{2}-\Sigma _{3}-\Sigma _{4}.\,\!}
由于仅有小数部分是准确的，想要抽取特定的小数位需要去掉最终结果的整数部分，并乘16来跳过这一个16进制位（理论上说在精度范围内这一位下面的几个小数位仍然是准确的）。
这个过程类似于長整數乘法（long multiplication 进位 没有被计算，而我们只关心最重要的位，计算机通常对很多位（32或者64）同时进行计算。存在一种极低的可能性，就是当极端情况出现，可能会将一个小数（比如1）加到999999999999999上，然后这个错误将会影响最重要的那个位。但在最终计算结果的时候，这种情况如果要发生，那也会是显而易见的。
这个算法被广为应用并有多种程序实现[ 6] [ 7] [ 8] [ 9]
编辑
这个算法无需自定义一种包含数千甚至数亿数字的“大数”数据类型。这种算法计算第n 位而无需 计算前n -1位，因此可以采用简单有效的数据类型。
这种算法对于计算π的第n 位（或者第n 位的附近几位）是最快最有效的，但使用大数据类型来计算π的前n 位的算法仍然比这个算法更快。
编辑
编辑
^ Bailey, David H. , Borwein, Peter B. Plouffe, Simon . On the Rapid Computation of Various Polylogarithmic Constants. Mathematics of Computation. April 1997, 66 (218): 903–913. doi:10.1090/S0025-5718-97-00856-9 ^ Controversy surrounding who among the three actually invented this algorithm . [2010-04-25 ] . （原始内容存档 于2010-04-17）. ^ Weisstein, Eric W. (编). BBP Formula . at MathWorld Wolfram Research, Inc. （英语） . ^ ANALYSIS OF PSLQ . [2011-11-21 ] . （原始内容存档 于2016-03-04）. ^ The Quest for Pi (PDF) . [2010-04-25 ] . （原始内容 (PDF) 存档于2011-09-27）. ^ A C++ implementation of the BBP algorithm for π（portable, SSE2 and OpenMP versions） . [2010-04-25 ] . （原始内容存档 于2010-11-27）. ^ (Python)| A Python implementation of the BBP algorithm for π . [2010-04-25 ] . （原始内容存档 于2016-03-04）. ^ A Ruby implementation of the BBP algorithm for π . [2018-04-24 ] . （原始内容存档 于2008-06-08）. ^ Computing π on a distributed cluster of computers . [2010-04-25 ] . （原始内容存档 于2010-02-05）. ^ D.J. Broadhurst, "Polylogarithmic ladders, hypergeometric series and the ten millionth digits of ζ(3) and ζ(5) （页面存档备份 ，存于互联网档案馆 ）", (1998) arXiv math.CA/9803067
编辑
![{\displaystyle \pi =\sum _{k=0}^{\infty }\left[{\frac {1}{16^{k}}}\left({\frac {4}{8k+1}}-{\frac {2}{8k+4}}-{\frac {1}{8k+5}}-{\frac {1}{8k+6}}\right)\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/af6bc360851499dd2ab2a90bee03fbe2040089d5)
![{\displaystyle \alpha =\sum _{k=0}^{\infty }\left[{\frac {1}{b^{k}}}{\frac {p(k)}{q(k)}}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fa85e25344327bcebd7461413dc9e28dc1a060ce)
![{\displaystyle P(s,b,m,A)=\sum _{k=0}^{\infty }\left[{\frac {1}{b^{k}}}\sum _{j=1}^{m}{\frac {a_{j}}{(mk+j)^{s}}}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/79b71ecf9ae82b1b067d955f13f3c46dc89bb0ed)
![{\displaystyle {\begin{aligned}\ln {\frac {9}{10}}&=-{\frac {1}{10}}-{\frac {1}{200}}-{\frac {1}{3\ 000}}-{\frac {1}{40\ 000}}-{\frac {1}{500\ 000}}-\cdots \\&=-\sum _{k=1}^{\infty }{\frac {1}{10^{k}\cdot k}}=-{\frac {1}{10}}\sum _{k=0}^{\infty }\left[{\frac {1}{10^{k}}}\left({\frac {1}{k+1}}\right)\right]\\&=-{\frac {1}{10}}P\left(1,10,1,(1)\right)\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0217a3eb1cb33c367d6de8127347627528695c67)
![{\displaystyle {\begin{aligned}\ln 2&={\frac {1}{2}}+{\frac {1}{2\cdot 2^{2}}}+{\frac {1}{3\cdot 2^{3}}}+{\frac {1}{4\cdot 2^{4}}}+{\frac {1}{5\cdot 2^{5}}}+\cdots \\&=\sum _{k=1}^{\infty }{\frac {1}{2^{k}\cdot k}}={\frac {1}{2}}\sum _{k=0}^{\infty }\left[{\frac {1}{2^{k}}}\left({\frac {1}{k+1}}\right)\right]\\&={\frac {1}{2}}P\left(1,2,1,(1)\right)\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f91cd70cdafbd1369f3a356b4f0945694ef7e390)
![{\displaystyle {\begin{aligned}\arctan {\frac {1}{b}}&={\frac {1}{b}}-{\frac {1}{b^{3}3}}+{\frac {1}{b^{5}5}}-{\frac {1}{b^{7}7}}+{\frac {1}{b^{9}9}}+\cdots \\&=\sum _{k=1}^{\infty }\left[{\frac {1}{b^{k}}}{\frac {\sin {\frac {k\pi }{2}}}{k}}\right]={\frac {1}{b}}\sum _{k=0}^{\infty }\left[{\frac {1}{b^{4k}}}\left({\frac {1}{4k+1}}+{\frac {-1}{4k+3}}\right)\right]\\&={\frac {1}{b}}P\left(1,b^{4},4,(1,0,-1,0)\right)\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/444f2d53c0535d4635eb27b9c50e3fac115a0cb6)
![{\displaystyle {\begin{aligned}\pi &=\sum _{k=0}^{\infty }\left[{\frac {1}{16^{k}}}\left({\frac {4}{8k+1}}-{\frac {2}{8k+4}}-{\frac {1}{8k+5}}-{\frac {1}{8k+6}}\right)\right]\\&=P\left(1,16,8,(4,0,0,-2,-1,-1,0,0)\right)\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6dc5b6efe7403a9af969ac126f3ca4baa62835f0)
![{\displaystyle \pi =\sum _{k=0}^{\infty }\left[{\frac {1}{16^{k}}}\left({\frac {120k^{2}+151k+47}{512k^{4}+1024k^{3}+712k^{2}+194k+15}}\right)\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ba24829fa8716fb94e690df022bb7a9e3e48ed23)
