贝利-波尔温-普劳夫公式

**貝利 - 波爾溫 - 普勞夫公式**（**BBP 公式**）提供了一個計算**圓周率 π** 的第 _n_ 位二進位數的 spigot 算法（spigot algorithm）。這個求和公式是在 1995 年由西蒙・普勞夫提出的，並以公布這個公式的論文作者大衛・貝利（David H. Bailey）、皮特・波爾溫（Peter Borwein）和普勞夫的名字命名。在論文發表之前，普勞夫已將此公式在他的網站上公布。這個公式是：

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

### π 的 BBP 公式

原始的 BBP π 求和公式是在 1995 年由 Plouffe 使用 PSLQ 算法（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]$

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

#### π 的 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 8_k_+ 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 8_k_+ 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 ^{3}2$ ，阿培裏常數 $\zeta (3)$ （其中 $\zeta (x)$ 是黎曼 ζ 函數）， $\pi ^{4}$ ， $\log ^{4}2$ ， $\log ^{5}2$ ， $\zeta (5)$ ，還有很多 $\pi $ 和 $\log 2$ 的不同冪次。這些結果主要是使用多重對數函數（polylogarithm ladder）得到的。

## 參考資料

## 外部連結

* Cook’s Class Contains Pi（頁面存檔備份，存於網際網路檔案館）