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

貝利-波爾溫-普勞夫公式（Bailey–Borwein–Plouffe formula）提供了一個計算圓周率π的第n位二進位數的spigot演算法（spigot algorithm）。這個求和公式是在1995年由西蒙·普勞夫提出的，並以公布這個公式的論文作者大衛·H·貝利 (數學家)、彼得·博溫和普勞夫的名字命名。在論文發表之前，普勞夫已將此公式在他的網站上公布。這個公式是：

$$\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）。由特定的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]$$

其中s, b和m是整數，$A =(a_1, a_2, \dots , a_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}$$

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

普勞夫也對這種形式的反正切函數的級數展開感興趣（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}$$

### 尋找新的等式

使用上述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使用整數關係偵查演算法（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}$$

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

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

這個等式可以用一個較為簡單的方式嚴格證明。

#### π的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，使得小數點恰好落在第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）。

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

$$4 \Sigma_1 - 2 \Sigma_2 - \Sigma_3 - \Sigma_4. \,\!$$

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

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

這個演算法被廣為應用並有多種程序實現。

#### 使用BBP演算法計算π的好處

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

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

## 推廣

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

## 外部連結

- Cook’s Class Contains Pi