查看“︁贝利-波尔温-普劳夫公式”︁的源代码

'''贝利-波尔温-普劳夫公式'''（'''BBP公式'''）提供了一个计算'''[[圓周率|圓周率π]]'''的第''n''位[[二进制]]数的{{le|spigot算法|spigot algorithm}}（{{lang|en|spigot algorithm}}）。这个[[总和|求和]]公式是在1995年由[[西蒙·普勞夫]]提出的，并以公布这个公式的论文作者[[大卫·贝利]]、{{le|彼得·博温|Peter Borwein}}和普勞夫的名字命名。在论文发表之前，普勞夫已将此公式在他的网站上公布<ref>{{cite journal
 | author = [[大卫·贝利|Bailey, David H.]], {{le|皮特·波爾溫|Peter Borwein|Borwein, Peter B.}}, and [[西蒙·普勞夫|Plouffe, Simon]]
 | title = On the Rapid Computation of Various Polylogarithmic Constants
 | journal = Mathematics of Computation
 | volume = 66 | issue = 218 | pages = 903–913
<!-- | url = http://crd.lbl.gov/~dhbailey/dhbpapers/digits.pdf   PDF is free at doi link, and is better formatted. -->
 | doi = 10.1090/S0025-5718-97-00856-9
 |date=April 1997}}</ref><ref>{{Cite web |url=http://groups.google.com/group/sci.math.symbolic/msg/5b7e62ce42ae0cb6 |title=Controversy surrounding who among the three actually invented this algorithm |accessdate=2010-04-25 |archive-date=2010-04-17 |archive-url=https://web.archive.org/web/20100417015611/http://groups.google.com/group/sci.math.symbolic/msg/5b7e62ce42ae0cb6 |dead-url=no }}</ref>。这个公式是： 
: <math> \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] </math>

这个公式的发现曾震惊学界。数百年来，求出π的第''n''位小数而不求出它的前''n''-1位曾被认为是不可能的。

自从这个发现以来，发现了更多的无理数常数的类似公式，它们都有一个类似的形式：
: <math>\alpha = \sum_{k = 0}^{\infty}\left[ \frac{1}{b^k} \frac{p(k)}{q(k)} \right]</math>
其中α是目标常数，''p''和''q''是整系数[[多项式]]，''b'' ≥ 2是整数的[[记数系统|数制]]。

这种形式的公式被称为'''BBP式公式'''（BBP-type formulas）<ref>{{cite mathworld| title=BBP Formula| urlname=BBPFormula}}</ref>。由特定的''p'',''q''和''b''可组合出一些著名的常数。但至今尚未找出一种系统的算法来寻找合适的组合，而已知的公式多是通过{{le|实验数学|Experimental mathematics}}得出的。

== 特例 ==
一个已经得出很多解的BBP式的特例是：
: <math> 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]</math>
其中''s'', ''b''和''m''是整数，<math>A =(a_1, a_2, \dots , a_m)</math>是一个整数向量。
使用P函数可得到一些解法的省略记法。

=== 已知的BBP式 ===
在BBP公式发现之前，就已经有些最简单的此类公式广为所知。他们使用P函数的省略记法如下：
: <math>\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}</math>
: <math>\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}</math>

普勞夫也对这种形式的[[反正切|反正切函數]]的級數展開感兴趣（P记法还可以泛化为''b''不是整数的情形）:
: <math>\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}</math>

=== 寻找新的等式 ===
使用上述P函数，令''s&nbsp;=&nbsp;1'', ''m&nbsp;>&nbsp;1''可得出已知的π的最简单公式。很多已知的公式是令b是2或3的幂，m是2的幂或者其他多因子的值，并令向量A等於零。在计算了一些类似的和之后，这类发现引发了使用计算机搜索类似线性组合的尝试。搜索程序为每个参数，s, b, m分别选择一个定义域，然后求和并检查值，并使用{{le|整数关系侦查算法|Integer relation algorithm}}（{{lang|en|integer relation algorithm}}），例如{{le|赫拉曼·弗古森|Helaman Ferguson}}（{{lang|en|Helaman Ferguson}}）的PSLQ算法，来找到一个向量''A''使得这些中间值可以加在一起得到一个著名的常数（''A''也可能是零）。

=== π的BBP公式 ===
原始的BBP π求和公式是在1995年由Plouffe使用{{le|整数关系侦查算法|Integer relation algorithm|PSLQ算法}}<ref>{{Cite web |url=http://www.ams.org/journals/mcom/1999-68-225/S0025-5718-99-00995-3/home.html |title=ANALYSIS OF PSLQ |accessdate=2011-11-21 |archive-date=2016-03-04 |archive-url=https://web.archive.org/web/20160304205518/http://www.ams.org/journals/mcom/1999-68-225/S0025-5718-99-00995-3/home.html |dead-url=no }}</ref>（{{lang|en|integer relation algorithm}}）发现的。它同样可由上述''P''函数表示：
: <math>\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}</math>
这个公式也可以使用下述两个多项式的商来表示：
: <math>\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]</math>

这个等式可以用一个较为简单的方式严格证明。<ref name=QuestForPi>{{Cite web |url=http://crd.lbl.gov/~dhbailey/dhbpapers/pi-quest.pdf |title=The Quest for Pi |accessdate=2010-04-25 |archive-date=2011-09-27 |archive-url=https://web.archive.org/web/20110927044106/http://crd.lbl.gov/~dhbailey/dhbpapers/pi-quest.pdf |dead-url=yes }}</ref>

==== π的BBP位抽取算法 ====
这个公式给出一个抽取π的[[十六进制]]位数值的算法。首先我们需要将这个公式化为：
: <math>\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)}</math>。

在使用此公式实现一个spigot算法之前，还需做一些改动。我们想要找出π在[[十六进制]]下的第''n''位，所以首先我们需要将该求和序列拆为n之前和n之后两部分。对于前述公式的第一项而言，
: <math>\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)}</math>。

我们再将等式两边乘以16<sup>&nbsp;''n''</sup>，使得小数点恰好落在第''n''位。
: <math>\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}</math>。

由于我们关心的是小数部分，我们观察这个式子的两项，会发现仅有第一项能够出现整数部分，而第二项不会出现整数部分。因为第二项中''k'' > ''n'',第二项中的分子一定小于分母。由此我们需要一个使用一种技巧来从第一项中除去整数。那就是要mod&nbsp;8''k''&nbsp;+&nbsp;1。于是原式第一项的小数部分的公式就是：
: <math>\sum_{k = 0}^{n} \frac{16^{n-k} \mod 8k+1}{8k+1} + \sum_{k = n + 1}^{\infty} \frac{16^{n-k}}{8k+1}</math>。

注意[[模]]运算将保证只有小数部分的和会被保留下来。想要快速有效的计算16<sup>&nbsp;''n'' - ''k''</sup>&nbsp;mod&nbsp;8''k''&nbsp;+&nbsp;1 ,可使用[[模幂運算]]（{{lang|en|Modular exponentiation}}）。

对公式的其余项使用相同的处理办法，并根据原式求出最后的结果。
: <math>4 \Sigma_1 - 2 \Sigma_2 - \Sigma_3 - \Sigma_4. \,\!</math>

由于仅有小数部分是准确的，想要抽取特定的小数位需要去掉最终结果的整数部分，并乘16来跳过这一个16进制位（理论上说在精度范围内这一位下面的几个小数位仍然是准确的）。

这个过程类似于長整數乘法（{{lang|en|long multiplication}}），但只需对一些中间列进行求和。由于有些[[进位]]没有被计算，而我们只关心最重要的位，计算机通常对很多位（32或者64）同时进行计算。存在一种极低的可能性，就是当极端情况出现，可能会将一个小数（比如1）加到999999999999999上，然后这个错误将会影响最重要的那个位。但在最终计算结果的时候，这种情况如果要发生，那也会是显而易见的。

这个算法被广为应用并有多种程序实现<ref>{{Cite web |url=http://www.nedprod.com/programs/portable/Pi/index.html |title=A C++ implementation of the BBP algorithm for π（portable, SSE2 and OpenMP versions） |accessdate=2010-04-25 |archive-date=2010-11-27 |archive-url=https://web.archive.org/web/20101127074508/http://nedprod.com/programs/portable/Pi/index.html |dead-url=no }}</ref><ref>{{Cite web |url=http://en.literateprograms.org/Pi_with_the_BBP_formula |title=(Python){{!}} A Python implementation of the BBP algorithm for π |accessdate=2010-04-25 |archive-date=2016-03-04 |archive-url=https://web.archive.org/web/20160304083253/http://en.literateprograms.org/Pi_with_the_BBP_formula |dead-url=no }}</ref><ref>{{Cite web |url=http://www.314159265358979323846264338327950288419716939937510.net/ |title=A Ruby implementation of the BBP algorithm for π |access-date=2018-04-24 |archive-date=2008-06-08 |archive-url=https://web.archive.org/web/20080608125003/http://www.314159265358979323846264338327950288419716939937510.net/ |dead-url=no }}</ref><ref>{{Cite web |url=https://issues.apache.org/jira/browse/HADOOP-5052 |title=Computing π on a distributed cluster of computers |accessdate=2010-04-25 |archive-date=2010-02-05 |archive-url=https://web.archive.org/web/20100205020542/http://issues.apache.org/jira/browse/HADOOP-5052 |dead-url=no }}</ref>。



==== 使用BBP算法计算π的好处 ====
这个算法无需自定义一种包含数千甚至数亿数字的“大数”数据类型。这种算法计算第''n''位而'''无需'''计算前''n''-1位，因此可以采用简单有效的数据类型。

这种算法对于计算π的第''n''位（或者第''n''位的附近几位）是最快最有效的，但使用大数据类型来计算π的前''n''位的算法仍然比这个算法更快。

== 推广 ==
D.J.布拉德赫斯特提出了一种BBP算法的泛化形式。<ref>D.J. Broadhurst, "[http://arxiv.org/abs/math.CA/9803067 Polylogarithmic ladders, hypergeometric series and the ten millionth digits of ζ(3) and ζ(5)] {{Wayback|url=http://arxiv.org/abs/math.CA/9803067 |date=20190713081616 }}", (1998) ''arXiv'' math.CA/9803067</ref>这种形式可以用于在接近线性时间和对数空间下求很多其他的常数。例如[[卡塔兰常数]]，<math>\pi^3</math>，<math>\log^32</math>，[[阿培裏常数]]<math>\zeta(3)</math>（其中<math>\zeta(x)</math>是[[黎曼ζ函數]]），<math>\pi^4</math>，<math>\log^42</math>，<math>\log^52</math>，<math>\zeta (5)</math>，还有很多<math>\pi</math>和<math>\log2</math>的不同幂次。这些结果主要是使用[[多重对数函数]]（{{lang|en|polylogarithm ladder}}）得到的。

== 参考资料 ==
{{reflist}}

== 外部链接 ==
* [http://rjlipton.wordpress.com/2009/03/15/cooks-class-contains-pi/ Cook’s Class Contains Pi]{{Wayback|url=http://rjlipton.wordpress.com/2009/03/15/cooks-class-contains-pi/ |date=20100610000539 }}

[[Category:圆周率算法]]