贝利-波尔温-普劳夫公式(BBP公式)提供了一个计算圓周率π的第n位二进制数的spigot算法(spigot algorithm)。这个求和公式是在1995年由西蒙·普勞夫提出的,并以公布这个公式的论文作者大卫·贝利(David H. Bailey)、皮特·波爾溫(Peter Borwein)和普勞夫的名字命名。在论文发表之前,普勞夫已将此公式在他的网站上公布[1][2]。这个公式是:
![\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)
这个公式的发现曾震惊学界。数百年来,求出π的第n位小数而不求出它的前n-1位曾被认为是不可能的。
自从这个发现以来,发现了更多的无理数常数的类似公式,它们都有一个类似的形式:
![\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)
其中α是目标常数,p和q是整系数多项式,b ≥ 2是整数的数制。
这种形式的公式被称为BBP式公式(BBP-type formulas)[3]。由特定的p,q和b可组合出一些著名的常数。但至今尚未找出一种系统的算法来寻找合适的组合,而已知的公式多是通过实验数学得出的。
一个已经得出很多解的BBP式的特例是:
![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)
其中s, b和m是整数,
是一个整数向量。
使用P函数可得到一些解法的省略记法。
已知的BBP式[编辑]
在BBP公式发现之前,就已经有些最简单的此类公式广为所知。他们使用P函数的省略记法如下:
![\begin{align} \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{align}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0217a3eb1cb33c367d6de8127347627528695c67)
![\begin{align} \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{align}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f91cd70cdafbd1369f3a356b4f0945694ef7e390)
普勞夫也对这种形式的反正切函數的級數展開感兴趣(P记法还可以泛化为b不是整数的情形):
![\begin{align} \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{align}](https://wikimedia.org/api/rest_v1/media/math/render/svg/444f2d53c0535d4635eb27b9c50e3fac115a0cb6)
寻找新的等式[编辑]
使用上述P函数,令s = 1, m > 1可得出已知的π的最简单公式。很多已知的公式是令b是2或3的幂,m是2的幂或者其他多因子的值,并令向量A等於零。在计算了一些类似的和之后,这类发现引发了使用计算机搜索类似线性组合的尝试。搜索程序为每个参数,s, b, m分别选择一个定义域,然后求和并检查值,并使用整数关系侦查算法(integer relation algorithm),例如赫拉曼·弗古森(Helaman Ferguson)的PSLQ算法,来找到一个向量A使得这些中间值可以加在一起得到一个著名的常数(A也可能是零)。
π的BBP公式[编辑]
原始的BBP π求和公式是在1995年由Plouffe使用PSLQ算法[4](integer relation algorithm)发现的。它同样可由上述P函数表示:
![\begin{align}\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{align}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6dc5b6efe7403a9af969ac126f3ca4baa62835f0)
这个公式也可以使用下述两个多项式的商来表示:
![\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)
这个等式可以用一个较为简单的方式严格证明。[5]
π的BBP位抽取算法[编辑]
这个公式给出一个抽取π的十六进制位数值的算法。首先我们需要将这个公式化为:
。
在使用此公式实现一个spigot算法之前,还需做一些改动。我们想要找出π在十六进制下的第n位,所以首先我们需要将该求和序列拆为n之前和n之后两部分。对于前述公式的第一项而言,
。
我们再将等式两边乘以16 n,使得小数点恰好落在第n位。
。
由于我们关心的是小数部分,我们观察这个式子的两项,会发现仅有第一项能够出现整数部分,而第二项不会出现整数部分。因为第二项中k > n,第二项中的分子一定小于分母。由此我们需要一个使用一种技巧来从第一项中除去整数。那就是要mod 8k + 1。于是原式第一项的小数部分的公式就是:
。
注意模运算将保证只有小数部分的和会被保留下来。想要快速有效的计算16 n - k mod 8k + 1 ,可使用模幂運算(Modular exponentiation)。
对公式的其余项使用相同的处理办法,并根据原式求出最后的结果。

由于仅有小数部分是准确的,想要抽取特定的小数位需要去掉最终结果的整数部分,并乘16来跳过这一个16进制位(理论上说在精度范围内这一位下面的几个小数位仍然是准确的)。
这个过程类似于長整數乘法(long multiplication),但只需对一些中间列进行求和。由于有些进位没有被计算,而我们只关心最重要的位,计算机通常对很多位(32或者64)同时进行计算。存在一种极低的可能性,就是当极端情况出现,可能会将一个小数(比如1)加到999999999999999上,然后这个错误将会影响最重要的那个位。但在最终计算结果的时候,这种情况如果要发生,那也会是显而易见的。
这个算法被广为应用并有多种程序实现[6][7][8][9]。
使用BBP算法计算π的好处[编辑]
这个算法无需自定义一种包含数千甚至数亿数字的“大数”数据类型。这种算法计算第n位而无需计算前n-1位,因此可以采用简单有效的数据类型。
这种算法对于计算π的第n位(或者第n位的附近几位)是最快最有效的,但使用大数据类型来计算π的前n位的算法仍然比这个算法更快。
D.J.布拉德赫斯特提出了一种BBP算法的泛化形式。[10]这种形式可以用于在接近线性时间和对数空间下求很多其他的常数。例如卡塔兰常数,
,
,阿培裏常数
(其中
是黎曼ζ函數),
,
,
,
,还有很多
和
的不同幂次。这些结果主要是使用多重对数函数(polylogarithm ladder)得到的。
参考资料[编辑]
- ^ Bailey, David H., Borwein, Peter B., and 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).
- ^ 埃里克·韦斯坦因. BBP Formula. MathWorld.
- ^ 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
外部链接[编辑]