円周率.jp (http://xn--w6q13e505b.jp/formula/bbp.html)

このページでは BBP の公式と呼ばれる公式やそれに類した公式を紹介する。 BBP の公式は 16 進数 (2 進数) 表記での $n$ 桁数目を $O(n)$ で求められる方式として知られる。

$\pi$ 全体を計算するという観点で見る場合、 単純な実装で考えると $O(n^2)$ で求めることができるが、 Binary Splitting 法や DRM 法を適用すれば $O(n(\log n)^3)$ で計算することもできる。

BBPの公式

\[ \pi = \sum_{n=0}^{\infty} \frac{1}{16^n}\left(\frac{4}{8n+1} - \frac{2}{8n+4} - \frac{1}{8n+5} - \frac{1}{8n+6}\right) \]

これが1995年9月19日[JB02]に発見された "BBP の公式" と呼ばれる公式である。 発見者 3 名（David Bailey, Peter Borwein, Simon Plouffe） の頭文字を取って名付けられた。[FB02] この式が正しい事の証明は高校数学程度の知識で十分可能である。 別ページにその流れを記した。

Bellard の公式[FT10]

\[ \pi = \frac{1}{64} \sum_{n=0}^{\infty} \frac{(-1)^n}{1024^n} \left( -\frac{32}{4n+1} - \frac{1}{4n+3} + \frac{256}{10n+1} -\frac{64}{10n+3} - \frac{4}{10n+5} - \frac{4}{10n+7} + \frac{1}{10n+9} \right) \]

現在知られている BBP 系統の公式類の中でも最も効率よく計算することができる。

その他の類似公式[FB02]

\[ \pi = \sum_{n=0}^{\infty} \frac{(-1)^n}{4^n} \left( \frac{2}{4n+1} + \frac{2}{4n+2} + \frac{1}{4n+3} \right) \quad \text{(by Adamchik and Wagon)} \]

BBP 公式の発見により、16 進数（2 進数） で特定の桁を算出するという計算記録も生み出された。 以下にその記録を示す。 なお、桁数は 16 進数で数えた小数点以下の桁数であり、 「結果」枠は該当桁以降の数字を示す。 つまり記録が 100 桁だとすると、100 桁目、101 桁目、… を表示している。 （ただし 1997 年の Bellard の記録については最初の '1' の部分だけ 2 進数であり、 上位 3bit 分は不明であったが、 後の計算で 16 進数は 7 ということが分かった。）

参考資料は [JB02] [FW01] [FW02] [FW03] [FW05] [JW08]

年	桁位置	結果	達成者	公式	備考
1995/9	10,000,000,000	9 21C73C683	Bailey, P. Borwein, Plouffe	BBP の公式
1996/7/6	50,000,000,000	C 1A10A49B3E 2B82A4404F 9193AD4EB	Bellard	Bellard の公式
1996/10/7	100,000,000,000	9 C381872D27 596F81D0E4 8B95A6C46	Bellard	Bellard の公式
1997/9/22	250,000,000,000	1 0FEE563B92 F0D22962B6 2DFD243160 8547A82	Bellard	Bellard の公式
1998/8/30	1,250,000,000,000	0 7E45733CC7 90B5B5979	Percival	Bellard の公式
1999/2/9	10,000,000,000,000	A 0F9FF371D1 7593E0	Percival	Bellard の公式
2000/9/11	250,000,000,000,000	E 6216B069CB 6C1D3	Percival	Bellard の公式
2010/9	500,000,000,000,000	0 E6C1294AED 40403F56D2 D764026265 BCA98511D0	施子和
2013/3/14	2,000,000,000,000,000	6 53728F1	Ed Karrels	Bellard の公式
2013/5/20	4,000,000,000,000,000	5 CC37DEC	Ed Karrels	Bellard の公式
2016/12/18	100,000,000,000,000,000	A 937EB59439 E485E	高橋大介	Bellard の公式

この公式を用いて $O(n)$ で $n$ ビット目を求める方法を紹介する。 計算方式の方式や計算量の計算では公式を以下の様に分離して考える。

\[ \pi = 4S(8,1) - 2S(8,4) - S(8,5) - S(8,6) \] \[ \text{ただし } S(a,b)=\sum_{n=0}^{\infty}\frac{1}{16^n}\frac{1}{an+b} \]

計算量

とりあえず繰り上がりを考慮せずに考えると、16 進数で $d$ 桁目まで求めるのに必要な各 $S$ の項数 $n$ は

\[ \begin{eqnarray} c\Delta S(a,b) \lt 16^{-d} && \\ n+ \log(an+b) \gt d + \log c && \\ n \gt d + \log c && (\log(an+b)=0 \text{と見なす}) \end{eqnarray} \]

ということになる。概数でいえば 16 進数で数えた桁数と同数程度の項数を計算すれば良い。

計算方式

$S$ の各項 $\dfrac{1}{16^n}\dfrac{1}{an+b}$ を 16 進数で表記する時、 $d$ 桁目以降の数字列は $16^{d-n-1}$ を $an+b$ で割った商の小数点以下の数字列になる。 そしてこの値は「$16^{d-n-1}$ を $an+b$ で割った余り $r$」 を $an+b$ で割った商の小数点以下になる。

\[ 16^{d-1}\Delta S \bmod 1 = \frac{16^{d-n-1}}{an+b} \bmod 1 = \frac{(an+b)q+r}{an+b} \bmod 1 = \frac{r}{an+b} \bmod 1 \quad (0 \leq r \lt an+b) \]

そのため、とりあえず $r$ を求め、$an+b$ で割る。 特に $d-n-1 \ge 0$ では

\[ AB \bmod C = (A \bmod C)(B \bmod C) \bmod C \]

ということを再帰的に利用する。 $a^n$ を $b$ で割った余りを求めるために、まずは $2^m\leq n \lt 2^{m+1}$ を満たす $m$ まで「2 のベキ乗」乗を $b$ で割った余りを求める。

\[ \begin{eqnarray} a^2 \bmod b &=& (a \bmod b)^2 \bmod b\\ a^4 \bmod b &=& (a^2 \bmod b)^2 \bmod b\\ a^8 \bmod b &=& (a^4 \bmod b)^2 \bmod b\\ & \vdots &\\ a^{2^m} \bmod b &=& (a^{2^{m-1}} \bmod b)^2 \bmod b\\ \end{eqnarray} \]

次に指数を足し合わせると $n$ になる $2^n$ の組み合わせのものを掛け合わせて余りを求める。 これは実質的に $n$ の 2 進数表示で 1 になっている箇所の値を掛け合わせることになる。 例えば

\[ \begin{eqnarray} 9 &=& 1001_{\text{[2進数]}} & \rightarrow & a^9 &=& a^8 a\\ 63 &=& 111111_{\text{[2進数]}} & \rightarrow & a^{63} &=& a^{32} a^{16} a^8 a^4 a^2 a\\ 137 &=& 10001001_{\text{[2進数]}} & \rightarrow & a^{137} &=& a^{128} a^8 a\\ 172 &=& 10101100_{\text{[2進数]}} & \rightarrow & a^{172} &=& a^{128} a^{32} a^8 a^4\\ 183 &=& 10110111_{\text{[2進数]}} & \rightarrow & a^{183} &=& a^{128} a^{32} a^{16} a^4 a^2 a \end{eqnarray} \]

というように計算する。

$d-n-1 \lt 0$ では $\Delta S\lt 1$ なので普通に割って結果を求める。 つまり

\[ 16^{d-1}S(a,b) = \sum_{n=0}^{\infty}\frac{16^{d-n-1}}{an+b} = \sum_{n=0}^{N-1} \frac{16^{d-n-1}}{an+b} + \sum_{n=N}^{\infty}\frac{16^{d-n-1}}{an+b} \]

と分け、1 つ目の $\sum$ を上記の方法で、 2 つ目の $\sum$ のうち数項を普通に計算することで $16^dS(a,b)$ の小数点以下の値を求める。 ここで $N$ は $16^{d-N^1} \lt aN+b$ を満たす最小の自然数である。 後半は 1 項進むごとに前の項の 1/16 以下になるので 20 項も計算すれば繰り上げを考慮しても 16 進数で 10 数桁は変わらない形になる。

1996 年、Simon Plouffe が $O(n^3 (\log n)^3)$ の計算量が必要になるものの、任意の基底で $\pi$ の $n$ 桁目だけを求める手法を開発し、 翌 1997 年に Bellard が必要な計算量を $O(n^2)$ まで減少させた。 [FT14] 2003 年、さらに X. Gourdon が別の計算原理を用いて $O(n^2(\log^2 n)^2 \log^3n / \log n)$、もしくは $O(m)$ のメモリを使って $O(n^2\log^3n \log^2n / (m\log^2(n/m)))$ の計算量となる方法を開発した。[FT14]

BBP の公式の一部については $\arctan$ の公式から変換して見つけることができる。 そのために $\log (1+x)$ の Maclaurin 展開を利用する。

\[ \log (1+x) = \sum_{n=1}^{\infty}\frac{-(-x)^n}{n} = x - \frac{x^2}{2} + \frac{x^3}{3} + \cdots \]

一方で arctan の式を以下のように変形する。

\[ \begin{eqnarray} \arctan \dfrac{1}{x} &=& {\rm Im} (\log (x + i) )\\ &=& {\rm Im} \left(\log \left(1 + \dfrac{i}{x}\right) \right)\\ &=& {\rm Im} \left(\sum_{n=1}^{\infty} -\dfrac{(-i)^n}{nx^n} \right)\\ &=& \sum_{n=0}^{\infty} \dfrac{1}{x^{4n}}\left(\dfrac{1}{(4n+1)x} - \dfrac{1}{(4n+3)x^3} \right)\\ &=& \sum_{n=0}^{\infty} \dfrac{(-1)^n}{(2n+1)x^{2n+1}} \end{eqnarray} \] \[ \begin{eqnarray} \arctan \dfrac{1}{x} &=& {\rm Im} (\log (x \pm 1 + \mp 1 + i) )\\ &=& {\rm Im} \left(\log \left(1 + \dfrac{\mp 1 + i}{x\pm 1}\right) \right)\\ &=& {\rm Im} \left(\sum_{n=1}^{\infty} -\dfrac{(\pm 1 - i)^n}{n(x\pm 1)^n} \right)\\ &=& \sum_{n=0}^{\infty} \dfrac{16^n}{(x\pm1)^{8n}}\left( \dfrac{1}{(8n+1)(x\pm 1)}\pm\dfrac{2}{(8n+2)(x\pm 1)^2}+ \dfrac{2}{(8n+3)(x\pm 1)^3}\right.\\ &&\left.-\dfrac{4}{(8n+5)(x\pm 1)^5}\mp \dfrac{8}{(8n+6)(x\pm 1)^6}-\dfrac{8}{(8n+7)(x\pm 1)^7}\right)\\ &=& \sum_{n=0}^{\infty} \dfrac{(-1)^n 4^n}{(x\pm1)^{4n}}\left( \dfrac{1}{(4n+1)(x\pm 1)}\pm\dfrac{2}{(4n+2)(x\pm 1)^2}+ \dfrac{2}{(4n+3)(x\pm 1)^3}\right)\\ \end{eqnarray} \]

最後の変形は $(\pm1-i)^2=\mp2i$ というのを利用している。また、わざと約分していない。

さてこの2つの式変形を利用することで、$x=2^k$ または $x=2^k\pm1$ のパターンのみが含まれる $\arctan$ 公式を BBP 公式に変換することができる。 あらためて $2^k=x, x\pm1$ と置き直すと

\[ \arctan \dfrac{1}{2^k} = \sum_{n=0}^{\infty} \dfrac{(-1)^n}{2^{2kn}}\dfrac{2^{-k}}{2n+1} \] \[ \arctan \dfrac{1}{2^k\pm1} = \sum_{n=0}^{\infty} \dfrac{(-1)^n}{2^{(4k-2)n}}\left( \dfrac{2^{-k}}{4n+1} \mp \dfrac{2^{-2k}}{2n+1} + \dfrac{2^{1-3k}}{4n+3}\right) \]

という形になる。ここで以下のように $x=2^k, 2^k\pm1$ に対して $1+x^2$ の素因数分解 を行い $\arctan$ の素数探索法を適用することで効率的に BBP 系公式を作り出すことができる。

\[ \begin{eqnarray} 1+2^2 &=& 5 &=& 5\\ 1+3^2 &=& 10 &=& 2\cdot5\\ 1+4^2 &=& 17 &=& 17\\ 1+5^2 &=& 26 &=& 2\cdot13\\ 1+7^2 &=& 50 &=& 2\cdot5^2\\ 1+8^2 &=& 65 &=& 5\cdot13\\ 1+9^2 &=& 82 &=& 2\cdot41\\ &\vdots& \end{eqnarray} \]

たとえば

\[ \begin{eqnarray} \pi &=& 4\arctan \dfrac{1}{2^1} + 4 \arctan \dfrac{1}{2^1+1}\\ &=& \sum_{n=0}^{\infty} \dfrac{(-1)^n}{4^n}\dfrac{2}{2n+1} +\sum_{n=0}^{\infty} \dfrac{(-1)^n}{4^n}\left( \dfrac{2}{4n+1}-\dfrac{1}{2n+1}+\dfrac{1}{4n+3}\right)\\ &=& \sum_{n=0}^{\infty} \dfrac{(-1)^n}{4^n}\left( \dfrac{2}{4n+1}+\dfrac{1}{2n+1}+\dfrac{1}{4n+3}\right)\\ &=& \dfrac{1}{4}\sum_{n=0}^{\infty} \dfrac{1}{16^n}\left( \dfrac{8}{8n+1}+\dfrac{4}{4n+1}+\dfrac{4}{8n+3} -\dfrac{2}{8n+5}-\dfrac{1}{4n+3}-\dfrac{1}{8n+7} \right)\\ \end{eqnarray} \] \[ \begin{eqnarray} \pi &=& 8\arctan \dfrac{1}{2^1} - 4 \arctan \dfrac{1}{2^3-1}\\ &=& \sum_{n=0}^{\infty} \dfrac{(-1)^n}{2^{2n}}\dfrac{2^2}{2n+1} - \sum_{n=0}^{\infty} \dfrac{(-1)^n}{2^{10n}}\left( \dfrac{2^{-1}}{4n+1}+\dfrac{2^{-4}}{2n+1}+\dfrac{2^{-6}}{4n+3}\right)\\ &=& \sum_{n=0}^{\infty} \dfrac{(-1)^n}{2^{10n}}\left( \dfrac{2^2}{10n+1} - \dfrac{1}{10n+3} + \dfrac{2^{-2}}{10n+5} -\dfrac{2^{-4}}{10n+7} + \dfrac{2^{-6}}{10n+9} \right)\\ && - \sum_{n=0}^{\infty} \dfrac{(-1)^n}{2^{10n}}\left( \dfrac{2^{-1}}{4n+1}+\dfrac{2^{-4}}{2n+1}+\dfrac{2^{-6}}{4n+3}\right)\\ &=& \dfrac{1}{64}\sum_{n=0}^{\infty} \dfrac{(-1)^n}{1024^n}\left( -\dfrac{32}{4n+1} - \dfrac{1}{4n+3} + \dfrac{256}{10n+1} - \dfrac{64}{10n+3} -\dfrac{4}{10n+5} - \dfrac{4}{10n+7} + \dfrac{1}{10n+9} \right)\\ \end{eqnarray} \]

という形で Bellard の公式も導出できる。 最後の行では $2n+1$ と $10n+5$ の通分をして纏めている。