π = k = 0 \sum \infty \frac{1}{1 6 ^{k}} (\frac{4}{8 k + 1} - \frac{2}{8 k + 4} - \frac{1}{8 k + 5} - \frac{1}{8 k + 6})

BBP の公式

Bailey-Borwein-Plouffe formula

望遠鏡のように、遠くの桁だけをのぞけないか。

16 進法なら、円周率の n 桁目を、前の桁を 1 つも求めずに計算できる級数。1995 年に計算機の探索で見つかり、桁の記録の検算に使われている。

年代: 1995年(発見)・1997年(論文)
人物: サイモン・プルーフ、デヴィッド・ベイリー、ピーター・ボールウェイン、ヘラマン・ファーガソン、ファブリス・ベラール、コリン・パーシヴァル、高橋大介

数学情報科学工学

逸話

遠くの桁をのぞく望遠鏡

1995 年、カナダのサイモン・フレーザー大学で、1 つの問いが持ち上がった。ある定数の小数で、決めた位置から先の数字を、そこまでの数字をすべて計算するより速く求められないか。共同研究者のデヴィッド・ベイリーは 2021 年の追悼文で、問いを立てたのは構成・実験数学センターの所長だったピーター・ボールウェインだと書いている。ベイリーはこの問いを「望遠鏡のように、数字の列の遠くをのぞけるか」と言い換えた。

当時は、π の $d$ 桁目だけを求めるのも、前の $d$ 桁をすべて求めるのと手間は変わらない、と広く思われていた。1997 年の論文の冒頭にもそうある。それまでの円周率の計算法は、どれも桁数とともに増える大きなメモリを要した。

手がかりは、昔から知られた $lo g 2 = \sum_{k = 1}^{\infty} 1/ (k 2^{k})$ だった。 $- lo g (1 - x) = x + x^{2} /2 + x^{3} /3 + \dots$ というテイラー展開に $x = 1/2$ を入れた式である。この級数を使うと、2 進法で書いた $lo g 2$ の決めた位置からの数字が、多くの桁を扱う特別なプログラムなしに、ふつうの計算機の数だけで求まる。だれが気づいたかは記録によって違う。ベイリーの 2006 年の解説はボールウェインとサイモン・プルーフの 2 人とし、2021 年の追悼文はボールウェインがすぐに見つけたと書く。プルーフ本人は 2003 年の投稿で、π の $n$ 桁目を求める式を 1974 年から探していたとし、1995 年 8 月 7 日ごろに自分が気づいたと書いている。

では、円の面積と円周率の π ではどうか。同じ形の式は、π にはまだ知られていなかった。追悼文によれば、何か月かの探索ののち、プルーフのプログラムがある関係を見つけ、式を整えると冒頭の式になった。

1997 年、ベイリーとボールウェイン兄弟とプルーフの 4 人は、解説「π の探求」に、この式にもその証明にも難しいところは何もないと書いた。こうした方法がこれほど長く見つからずにいたのは率直に言って驚きで、たとえばオイラーが見つけていてもおかしくなかった。今の研究者に有利な点は、進んだ計算機の技術だけだ、と続けている。ベイリーは 2006 年に、この式は計算機のプログラムが見つけた π の重要な新しい公式の、ほぼ確実に最初の例だと書いた。

起源

整数の関係を探す

BBP の公式は、証明から見つかったのではない。1997 年の論文は、この式を「ひらめきの混じった当て推量と、整数関係アルゴリズム PSLQ による大がかりな探索」で見つけたと書き、自分たちの証明について「発見の道筋をすっかり隠している」と添えた。

整数関係アルゴリズムは、実数 $x_{1}, \dots, x_{n}$ を高い精度で与えると、 $a_{1} x_{1} + \dots + a_{n} x_{n} = 0$ となる整数 $a_{1}, \dots, a_{n}$ (すべてが 0 ではない)を探す手順である。PSLQ はヘラマン・ファーガソンとベイリーがまとめた方法で、2 人の報告書(RNR-91-032)は 1992 年 7 月 14 日の日付を持つ。スティーブ・アーノを加えた 3 人の解析は、1999 年に『Mathematics of Computation』に載った。

$lo g 2$ と同じように桁を取り出せる形の定数を並べ、π との間に整数の関係があるかを調べる。関係が見つかれば、それを π について解くと、桁を取り出せる π の式になる(ベイリーの追悼文)。たとえば $S_{j} = \sum_{k = 0}^{\infty} 1/ {1 6^{k} (8 k + j)}$ を高い精度で求めると、π との間に次の関係が成り立つ。

π = 4 S_{1} - 2 S_{4} - S_{5} - S_{6}

だれがどの道具で見つけたかも、記録によって違う。ベイリーの 2006 年の解説と追悼文では、プルーフが、ベイリーの書いた PSLQ のプログラムで探索して見つけた。プルーフは 2003 年 6 月のニュースグループへの投稿で、PSLQ のプログラムを書いたのはベイリーとファーガソンだとしたうえで、発見をもたらしたのは手元にあった、数式処理ソフトの Pari-GP と Maple をつなぐプログラムだと書いた。同じ投稿で、ベイリーが加わったのは発見の 2 か月後で、ボールウェインとベイリーをこの式とアルゴリズムの共同の創始者として受け入れたことが人生で最大の誤りだったと述べている。3 人の名前で出た論文は、PSLQ による発見としている。どちらの記述が細部まで正しいかを決める記録は、確かめられなかった。

発見は 1995 年のうちに世に出た。プルーフの回想によれば、1995 年 10 月にはカナダの新聞『グローブ・アンド・メール』に記事が出た。論文は、1996 年 4 月の NASA エイムズ研究センターの技術報告(NAS-96-016)を経て、『Mathematics of Computation』の 1997 年 4 月号(66 巻 218 号)に載った。計算はエイムズ研究センターのワークステーションの空き時間で行い、16 進法の 100 億桁目からの 14 桁を求めた。

論文によれば、発表のすぐ後に、ファーガソンやドナルド・クヌースら 6 人の同僚が、同じ型の別の π の式を寄せた。それらはどれも、BBP の公式に、0 に等しい次の級数の有理数倍を足したものだった。

0 = i = 0 \sum \infty \frac{1}{1 6 ^{i}} (- \frac{8}{8 i + 1} + \frac{8}{8 i + 2} + \frac{4}{8 i + 3} + \frac{8}{8 i + 4} + \frac{2}{8 i + 5} + \frac{2}{8 i + 6} - \frac{1}{8 i + 7})

証明

積分で確かめ、桁を取り出す

$S_{j} = \sum_{k = 0}^{\infty} \frac{1}{1 6 ^{k} ( 8 k + j )}$ とおくと、示す式は $π = 4 S_{1} - 2 S_{4} - S_{5} - S_{6}$ である。筋は 1997 年の論文と解説「π の探求」に沿う。

$0 \leq x \leq 1/ 2$ では $x^{8} \leq 1/16$ なので、等比級数の和の式が使え、この範囲で一様に収束する。
$\frac{1}{1 - x ^{8}} = k = 0 \sum \infty x^{8 k}$
両辺に $x^{j - 1}$ を掛け、項ごとに $0$ から $1/ 2$ まで積分する。 $(1/ 2)^{8 k} = 1/1 6^{k}$ である。
$\int_{0}^{1/ 2} \frac{x ^{j - 1}}{1 - x ^{8}} d x = k = 0 \sum \infty \frac{( 1/ 2 ) ^{8 k + j}}{8 k + j} = \frac{S _{j}}{2 ^{j /2}}$
$j = 1, 4, 5, 6$ の 4 つに $2^{j /2}$ と係数を掛けて足す。 $4 \cdot 2^{1/2} = 42$ 、 $2 \cdot 2^{2} = 8$ 、 $2^{5/2} = 42$ 、 $2^{3} = 8$ である。
$4 S_{1} - 2 S_{4} - S_{5} - S_{6} = \int_{0}^{1/ 2} \frac{4 2 - 8 x ^{3} - 4 2 x ^{4} - 8 x ^{5}}{1 - x ^{8}} d x$
$y = 2 x$ と置き換える。 $d x = d y / 2$ 、 $1 - x^{8} = (16 - y^{8}) /16$ で、分子は $2 (4 - 2 y^{3} - y^{4} - y^{5})$ になる。
$= \int_{0}^{1} \frac{16 ( 4 - 2 y ^{3} - y ^{4} - y ^{5} )}{16 - y ^{8}} d y$
分母と分子を因数分解する。 $(y^{2} + 2) (y^{2} + 2 y + 2)$ が共通である。
$16 - y^{8} = - (y^{2} - 2) (y^{2} + 2) \cdot (y^{2} + 2 y + 2) (y^{2} - 2 y + 2) 4 - 2 y^{3} - y^{4} - y^{5} = - (y - 1) (y^{2} + 2) (y^{2} + 2 y + 2)$
約分して、部分分数に分ける。右辺を通分すると、分子は $4 y (y^{2} - 2 y + 2) - (4 y - 8) (y^{2} - 2) = 16 y - 16$ になる。
$\frac{16 ( y - 1 )}{( y ^{2} - 2 ) ( y ^{2} - 2 y + 2 )} = \frac{4 y}{y ^{2} - 2} - \frac{4 y - 8}{y ^{2} - 2 y + 2}$
1 つ目は対数になる。 $0 \leq y \leq 1$ で $y^{2} - 2 < 0$ である。
$\int_{0}^{1} \frac{4 y}{y ^{2} - 2} d y = [2 lo g (2 - y^{2})]_{0}^{1} = - 2 lo g 2$
2 つ目は、分子を $2 (2 y - 2) - 4$ と分け、 $y^{2} - 2 y + 2 = (y - 1)^{2} + 1$ と $\frac{d}{d y} arctan (y - 1) = \frac{1}{( y - 1 ) ^{2} + 1}$ を使う。
$\int_{0}^{1} \frac{4 y - 8}{y ^{2} - 2 y + 2} d y = [2 lo g (y^{2} - 2 y + 2)]_{0}^{1} - [4 arctan (y - 1)]_{0}^{1} = - 2 lo g 2 - π$
7 から 8 を引くと、 $lo g 2$ が消えて π だけが残る。これで $π = 4 S_{1} - 2 S_{4} - S_{5} - S_{6}$ が示された。
$(- 2 lo g 2) - (- 2 lo g 2 - π) = π$
ここから桁の取り出しである。16 進法で小数点から $d + 1$ 番目から先の数字は、 $1 6^{d} π$ の小数部分 ${1 6^{d} π}$ の数字である(図示の 1 枚目)。小数部分は、整数を足し引きしても変わらない。
${1 6^{d} π} = {4 {1 6^{d} S_{1}} - 2 {1 6^{d} S_{4}} - {1 6^{d} S_{5}} - {1 6^{d} S_{6}}}$
$1 6^{d} S_{j}$ の項を $k \leq d$ と $k > d$ に分ける。 $k \leq d$ の項 $1 6^{d - k} / (8 k + j)$ は整数の部分を捨ててよいので、分子を $8 k + j$ で割った余りに替えられる。
${1 6^{d} S_{j}} = {k = 0 \sum d \frac{1 6 ^{d - k} mod ( 8 k + j )}{8 k + j} + k = d + 1 \sum \infty \frac{1 6 ^{d - k}}{8 k + j}}$
余り $1 6^{e} mod m$ は繰り返し二乗法で求める。2 乗して $m$ で割ることを、 $e$ を 2 進法で書いた桁に従って繰り返す(図示の 2 枚目)。掛け算は $2 lo g_{2} e$ 回より少なく、途中の数は $m^{2}$ を超えない。合同式で書くと次のとおりである。
$1 6^{29} = 1 6^{16} \cdot 1 6^{8} \cdot 1 6^{4} \cdot 16 \equiv 16 \cdot 37 \cdot 18 \cdot 16 \equiv 18 (mod 41)$
$k > d$ の残りは、1 項ごとに 16 分の 1 より小さくなる。倍精度の数の刻みより小さくするには、14 項ほどで足りる。
$k = d + 1 \sum \infty \frac{1 6 ^{d - k}}{8 k + j} < \frac{1}{8 d + 8 + j} i = 1 \sum \infty \frac{1}{1 6 ^{i}} = \frac{1}{15 ( 8 d + 8 + j )}$
割り算と和は浮動小数点の数で行うので、丸めの誤差がたまる。論文は、1 回の演算の誤差が最後の 1 ビットまでなら、壊れるのは多くても $lo g_{2} (2 n)$ ビットと見積もった。倍精度の数の 53 ビットで $n = 1 0^{6}$ なら、控えめに見ても 32 ビット、16 進法で 8 桁が残る。
$53 - lo g_{2} (2 \times 1 0^{6}) = 53 - 20.93 \dots \approx 32$
余りの計算は正確でなければならない。途中の数が倍精度の数で正確に表せる $2^{53}$ 未満に収まる位置までしか、ふつうに書いたプログラムでは扱えない。ベイリーの 2006 年の解説は、これをおよそ $1.18 \times 1 0^{7}$ としている。
$(8 d + 6)^{2} < 2^{53} ⟺ d \leq 11863282$
結果は、位置を 1 つずらしてもう一度計算して確かめる。すべての演算が変わるので、数字が 1 つずれて重なれば、正しいと判断できる。論文の表は、こうして 2 回ずつ計算された。
${1 6^{d} π} = {16 {1 6^{d - 1} π}}$

応用

和としては遅い。 そのまま足すと、BBP の公式は 1 項ごとに約 1.2 桁ずつ正しくなる( $lo g_{10} 16 = 1.204 \dots$ )。 $k = 0$ の項だけでは $4 - 1/2 - 1/5 - 1/6 = 3.1333 \dots$ で、8 項では $3.14159265358897 \dots$ となり、小数第 11 位まで合う。1 項で 14 桁ずつ進むチュドノフスキーの公式には遠く及ばないので、何兆桁の記録を頭から計算するのには使われない。

100 万桁目の数字。 16 進法で 100 万桁目からの 14 桁は 26C65E52CB4593 である。BBP の論文の表にある数字で、取り出しで求めても、頭から 100 万桁あまりを計算しても、同じ 14 桁になる。取り出しに要るのは 400 万回の累乗の計算である。ベイリーは 2006 年に、100 万桁目からの 10 桁は当時のパソコンで 5 秒ほどで求まると書いた。1000 万桁目からの 17AF5863EFED8D も、論文の表にある。倍精度の数だけで取り出してみると、100 万桁目では 10 桁、1000 万桁目では 8 桁までが正しく、その先は丸めの誤差で崩れる。

記録の検算。 何兆桁の記録では、最後の桁が正しいかを BBP 型の式で別に確かめる。2009 年 12 月 31 日に約 2 兆 7000 億桁を発表したファブリス・ベラールは、別の式で全部を計算し直す代わりに、結果の最後の 2 進の桁を、BBP の方法と自分の式で求めて照らし合わせた。2010 年の近藤茂とアレクサンダー・イーの 5 兆桁では、BBP の公式とベラールの式を 2 台の計算機で同時に走らせ、16 進法で 4 兆 1524 億 1011 万 8610 桁目で終わる 32 桁を求めた。66 時間と 64 時間かかった。最後の桁が合っても、途中の桁が機械の誤りで変わっている心配は残る。そこでどちらの記録も、最後の掛け算を、ある素数で割った余りで照らす検査を組み合わせている。

遠くの桁の記録。 1997 年 1 月、ベラールは、1 項ごとに 2 進法で 10 桁ずつ進む同じ型の式を見つけた。BBP の公式より約 43 % 速い。ベラールはこの式で、2 進法の 1 兆桁目からの 152 ビットを、20 台あまりのワークステーションを使って 12 日で求めた。1998 年には、サイモン・フレーザー大学の 17 歳の学生コリン・パーシヴァルが、計算を世界中のパソコンに分ける PiHex を始めた。2000 年 9 月には、56 か国の 1734 台の空き時間を使い、2 進法の 1000 兆桁目が 0 であることを突き止めた。2010 年には、ヤフーのクラウド計算のチームが、1000 台の計算機の空き時間を 23 日使い、2000 兆桁目も 0 だと発表した。2018 年には、筑波大学の高橋大介が、ベラールの式で 16 進法の 10 京桁目からの数字を求めた論文を出している。

10 進法では。 10 進法で同じことができる式は見つかっていない。1997 年の論文は、16 を 2 から 128 までの数に替えた同じ形の式を PSLQ で探し、係数の大きさが $1 0^{12}$ 以下の範囲には、16 のほかに無いことを確かめた。前の桁なしで 10 進法の $n$ 桁目を求める方法は、プルーフが 1996 年 11 月に書いた論文で示され、手間は $n^{3} (lo g n)^{3}$ に比例した。ベラールは 1997 年 1 月に、これを $n^{2}$ に比例するところまで改良した。どちらも 16 進法の取り出しより重い。

図示

π の 16 進法の数字に 16⁴ を掛けると、数字はそのままで小数点が 4 桁右へ動く。整数の部分を捨てた小数部分の先頭が、5 桁目の数字になる。

16²⁹ を 41 で割った余りを、2 乗して割ることの繰り返しで求める。29 = 16 + 8 + 4 + 1 に当たる 4 つを掛け合わせ、7 回の掛け算で 18 を得る。

BBP 型の式で取り出した、π の数字の位置の記録(縦は 2 進法の桁で、対数の目盛り)。破線は、頭から全部を計算した最高の桁数を 2 進法に直したもの。

1 枚目は証明の手順 10 にあたる。2 枚目の 41 は、 $d = 34$ の $k = 5$ 、 $j = 1$ の項の $8 k + j$ である。3 枚目の値は、BBP の 1996 年の報告(16 進法の位置を 2 進法に直した)、ベラール・PiHex・ヤフーの発表、高橋の論文の題から取った。取り出した位置は、頭から全部の記録の数百倍先まで届いている。

シミュレーション

つまみで位置 $n$ を選ぶと、BBP の公式でその位置から 8 桁だけを計算する。4 つの小さな文字盤は和 ${1 6^{d} S_{j}}$ の小数部分を指し、それらを組み合わせた大きな文字盤の針が指す目盛りが $n$ 桁目の数字になる。上の列はマチンの公式で頭から 10050 桁を計算したもので、一致した桁を縦の線で結ぶ。それより遠くは「1 つ前からも計算」で、位置 $n - 1$ から取り出した列が 1 つずれて重なるかを確かめる。

広がる学問

数学 BBP 型の式は $π^{2}$ や $(lo g 2)^{2}$ などにもあり、1997 年の論文に並んでいる。式を計算機の探索で見つけ、あとから証明する「実験数学」の代表例でもある。ベイリーによれば、この式があることは、π が正規数か(同じ長さのどの数字の並びも同じ割合で現れるか)という未解決の問いにつながり、その筋の考えから、π ではないが具体的な数の族が正規数であることが証明された。
情報科学 位置 $n$ が大きくなっても、要るメモリはほとんど増えず、時間はほぼ $n$ に比例して増える。論文は計算量の言葉で、2 進法の π が、メモリは $lo g n$ の多項式、時間はほぼ $n$ に比例する計算の組(SC*)に入ると述べた。16 進法の数字を 10 進法に直すことが同じ条件でできるかは、論文の時点で分かっていなかった。
工学何か月も続く計算は、機械の誤りが起きるものとして組み立て、BBP 型の式で照らして確かめる。5 兆桁の記録の報告は、これほど長く多くの部品を使う計算では、故障は起こりうるのではなく必ず起こると書く。ヤフーの計算は、位置を 4 ビットずらした 2 回の計算が重なる部分だけを結果とし、計算機や通信の誤りを見逃さないようにした。

参考文献

D. H. Bailey, P. B. Borwein, S. Plouffe, "On the rapid computation of various polylogarithmic constants", Mathematics of Computation 66 (1997), 903-913(NASA Ames Research Center, NAS Technical Report NAS-96-016, 1996 年 4 月)
D. H. Bailey, J. M. Borwein, P. B. Borwein, S. Plouffe, "The Quest for Pi", Math. Intelligencer 19(1) (1997), 50-57
D. H. Bailey, "The BBP Algorithm for Pi"(2006 年 9 月 17 日); D. H. Bailey, "Peter Borwein: A philosopher's mathematician"(2021 年 9 月 9 日)
S. Plouffe, "The story behind a formula for Pi", sci.math.symbolic(2003 年 6 月 22 日付); S. Plouffe, "On the computation of the n^th decimal digit of various transcendental numbers", arXiv:0912.0303(1996 年 11 月 30 日の未発表の論文の再録)
H. R. P. Ferguson, D. H. Bailey, "A Polynomial Time, Numerically Stable Integer Relation Algorithm", RNR Technical Report RNR-91-032(1992 年 7 月 14 日); H. R. P. Ferguson, D. H. Bailey, S. Arno, "Analysis of PSLQ, an integer relation finding algorithm", Math. Comp. 68 (1999), 351-369
F. Bellard, "A new formula to compute the n'th binary digit of pi"(1997 年 1 月 20 日)、"Computation of the n'th digit of pi in any base in O(n^2)"(1997 年 1 月 10 日)、"Computation of 2700 billion decimal digits of Pi using a Desktop Computer"(2010 年 2 月 11 日第 4 版); T.-W. Sze, "The Two Quadrillionth Bit of Pi is 0! Distributed Computation of Pi with Apache Hadoop", arXiv:1008.3171(2010)
A. J. Yee, S. Kondo, 5 兆桁の記録の詳細(https://www.numberworld.org/misc_runs/pi-5t/details.html 、2026 年 10 月 24 日参照); D. Takahashi, "Computation of the 100 quadrillionth hexadecimal digit of π on a cluster of Intel Xeon Phi processors", Parallel Computing 75 (2018), 1-10