貝利-波爾溫-普勞夫公式

貝利-波爾溫-普勞夫公式（BBP公式）提供了一個計算圓周率π的第n位二進制數的spigot算法（英語：spigot algorithm）（spigot algorithm）。這個求和公式是在1995年由西蒙·普勞夫提出的，並以公佈這個公式的論文作者大衛·貝利、彼得·博溫（英語：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]

這個公式的發現曾震驚學界。數百年來，求出π的第n位小數而不求出它的前n-1位曾被認為是不可能的。

自從這個發現以來，發現了更多的無理數常數的類似公式，它們都有一個類似的形式：

\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可組合出一些著名的常數。但至今尚未找出一種系統的算法來尋找合適的組合，而已知的公式多是通過實驗數學（英語：Experimental mathematics）得出的。

特例

一個已經得出很多解的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]

其中s, b和m是整數， $A=(a_{1},a_{2},\dots ,a_{m})$ 是一個整數向量。使用P函數可得到一些解法的省略記法。

已知的BBP式

在BBP公式發現之前，就已經有些最簡單的此類公式廣為所知。他們使用P函數的省略記法如下：

{\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}}

{\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不是整數的情形）:

{\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）（integer relation algorithm），例如赫拉曼·弗古森（英語：Helaman Ferguson）（Helaman Ferguson）的PSLQ算法，來找到一個向量A使得這些中間值可以加在一起得到一個著名的常數（A也可能是零）。

π的BBP公式

原始的BBP π求和公式是在1995年由Plouffe使用PSLQ算法（英語：Integer relation algorithm）^[4]（integer relation algorithm）發現的。它同樣可由上述P函數表示：

{\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}}

這個公式也可以使用下述兩個多項式的商來表示：

\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]

π的BBP位抽取算法

這個公式給出一個抽取π的十六進制位數值的算法。首先我們需要將這個公式化為：

\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之後兩部分。對於前述公式的第一項而言，

\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位。

\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。於是原式第一項的小數部分的公式就是：

\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）（Modular exponentiation）。

對公式的其餘項使用相同的處理辦法，並根據原式求出最後的結果。

4\Sigma _{1}-2\Sigma _{2}-\Sigma _{3}-\Sigma _{4}.\,\!

由於僅有小數部分是準確的，想要抽取特定的小數位需要去掉最終結果的整數部分，並乘16來跳過這一個16進制位（理論上說在精度範圍內這一位下面的幾個小數位仍然是準確的）。

這個過程類似於長整數乘法（long multiplication），但只需對一些中間列進行求和。由於有些進位沒有被計算，而我們只關心最重要的位，計算機通常對很多位（32或者64）同時進行計算。存在一種極低的可能性，就是當極端情況出現，可能會將一個小數（比如1）加到999999999999999上，然後這個錯誤將會影響最重要的那個位。但在最終計算結果的時候，這種情況如果要發生，那也會是顯而易見的。

這個算法被廣為應用並有多種程序實現^[6]^[7]^[8]^[9]。

使用BBP算法計算π的好處

這個算法無需自定義一種包含數千甚至數億數字的「大數」數據類型。這種算法計算第n位而無需計算前n-1位，因此可以採用簡單有效的數據類型。

這種算法對於計算π的第n位（或者第n位的附近幾位）是最快最有效的，但使用大數據類型來計算π的前n位的算法仍然比這個算法更快。

推廣

D.J.布拉德赫斯特提出了一種BBP算法的泛化形式。^[10]這種形式可以用於在接近線性時間和對數空間下求很多其他的常數。例如卡塔蘭常數， $\pi ^{3}$ ， $\log ^{3}2$ ，阿培裏常數 $\zeta (3)$ （其中 $\zeta (x)$ 是黎曼ζ函數）， $\pi ^{4}$ ， $\log ^{4}2$ ， $\log ^{5}2$ ， $\zeta (5)$ ，還有很多 $\pi$ 和 $\log 2$ 的不同冪次。這些結果主要是使用多重對數函數（polylogarithm ladder）得到的。