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

BBP系公式

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

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

公式

BBPの公式

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

この流れの公式が数々発見される元になった公式. 発見者3名(David Bailey, Peter Borwein, Simon Plouffe) の頭文字を取って BBP の公式と呼ばれるようになった.[FB02]

この公式は発見されたのはごく最近 (1995年9月19日[JB02])ではあるが, 式が正しい事の証明は高校数学程度の知識で十分可能である. 別ページにその流れを記した.

その他の公式

\[ \pi = \sum_{k=0}^{\infty} \frac{(-1)^k}{4^k} \left( \frac{2}{4k+1} + \frac{2}{4k+2} + \frac{1}{4k+3} \right) \quad{\rm \ (by\ Adamchik\ and\ Wagon)} \] [FB02] \[ \pi = \frac{1}{64} \sum_{k=0}^{\infty} \frac{(-1)^k}{1024^k} \left( -\frac{32}{4k+1} - \frac{1}{4k+3} + \frac{256}{10k+1} -\frac{64}{10k+3} - \frac{4}{10k+5} - \frac{4}{10k+7} + \frac{1}{10k+9} \right) \quad{\rm \ (by\ F. Bellard)} \] [FT10]

部分計算の記録

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

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

桁位置 結果 達成者 公式 備考
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 の公式

計算方式

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

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

計算量

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

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

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

計算方式

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

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

そのため,とりあえず $r$ を求め,$ak+b$ で割る. 特に $d-k-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^k$ の組み合わせのものを掛け合わせて余りを求める. これは実質的に $n$ の 2 進数表示で 1 になっている箇所の値を掛け合わせることになる. 例えば

9 = 1001[2進数] $a^9=a^8 a$
63 =111111[2進数] $a^{63} = a^{32} a^{16} a^8 a^4 a^2 a$
137 = 10001001[2進数] $a^{137} = a^{128} a^8 a$
172 = 10101100[2進数] $a^{172} = a^{128} a^{32} a^8 a^4$
183 = 10110111[2進数] $a^{183} = a^{128} a^{32} a^{16} a^4 a^2 a$

というように計算する.

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

\[ 16^{d-1}S(a,b) = \sum_{k=0}^{\infty}\frac{16^{d-k-1}}{ak+b} = \sum_{k=0}^{K-1} \frac{16^{d-k-1}}{ak+b} + \sum_{k=K}^{\infty}\frac{16^{d-k-1}}{ak+b} \]

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

$\pi$ 以外の公式

これらの応用で,$\pi$ 以外の値を求める公式も求められた.

\[ \ln \left(\frac{9}{10}\right) = - \sum_{k=0}^{\infty}\frac{1}{10^k}\frac{1}{k} \]

任意基底

1996 年,Simon Pluffe が $O(n^3 (\log n)^3)$ の計算量が必要になるものの,任意の基底で $\pi$ の $n$ 桁目だけを求める手法を開発し, 翌 1997 年に Bellard が必要な計算量を $O(n^2)$ まで減少させた. [FT14]

arctan 公式から BBP 公式への変換

$\log$ の Taylor 展開

\[ \log(1+x) = \sum_{k=1}^{\infty} \dfrac{-(-x)^k}{k} \]

を利用することで、容易に arctan 公式を BBP 公式に変換することができる。 変換式としては

\[ \arctan\left(\dfrac{1}{a}\right) = \Im \log\left(1 + \dfrac{i}{a}\right) = \sum_{k=0}^{\infty} \frac{(-1)^k}{a^{2k+1}}\frac{1}{2k+1} \] \[ \arctan\left(\dfrac{1}{a\pm1}\right) = \Im \log\left(a\pm1 + i\right) = \Im \log\left(1 + \dfrac{\pm1+i}{a}\right) = \sum_{k=0}^{\infty} \frac{(-4)^k}{a^{4k+3}} \left(\frac{a^2}{4k+1}\mp\frac{a}{2k+1}+\frac{2}{4k+3}\right) \]

という形になっており、特に $a$ が 2 の冪乗になる式が利用しやすい。 具体的な数を挙げると

\[ \begin{eqnarray} \arctan\left(\dfrac{1}{2}\right) &=& \sum_{k=0}^{\infty} \dfrac{(-1)^k}{2^{2k+1}}\dfrac{1}{2k+1} \\ \arctan\left(\dfrac{1}{3}\right) &=& \sum_{k=0}^{\infty} \dfrac{(-1)^k}{2^{2k+2}}\left(\dfrac{1}{4k+1}-\dfrac{1}{2k+1}+\dfrac{1}{4k+3}\right)\\ &=& \sum_{k=0}^{\infty} \dfrac{(-1)^k}{2^{6k+5}}\left(\dfrac{8}{4k+1}+\dfrac{2}{2k+1}+\dfrac{1}{4k+3}\right)\\ \arctan\left(\dfrac{1}{7}\right) &=& \sum_{k=0}^{\infty} \dfrac{(-1)^k}{2^{10k+8}}\left(\dfrac{2^5}{4k+1}+\dfrac{2^2}{2k+1}+\dfrac{1}{4k+3}\right)\\ \end{eqnarray} \]

といった式になる。ここから例えば

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

という形で Bellard の公式が導かれる。

ちなみにこの arctan からの変換系では 2 の冪乗を基底に取る BBP 公式しか導けないことが証明されていて、 arctan に限らないとしても恐らく他の基底では不可能だろうと予想されている。 [FT14]