ベイリー＝ボールウェイン＝プラウフの公式

ベイリー＝ボールウェイン＝プラウフの公式（ベイリー＝ボールウェイン＝プラウフのこうしき、英: Bailey–Borwein–Plouffe formula）あるいはBBP公式は、1995年にサイモン・プラウフによって発見された円周率 π に関する以下の公式である。

π = ∑ k = 0 ∞ [ 1 16 k ( 4 8 k + 1 − 2 8 k + 4 − 1 8 k + 5 − 1 8 k + 6 ) ] {\displaystyle \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]}
（BBP公式）[1]

名称はBBP公式に関する論文 (Bailey, Borwein & Plouffe 1997) の著者らデイヴィッド・H・ベイリー（英語版）、ピーター・ボールウェイン（英語版）、プラウフに因む。BBP公式は論文公開以前にもプラウフが自身のサイトで紹介していた[2]。

BBP公式は、先行する桁を計算せずに π の十六進法で n 桁目（二進法で 4n 桁目）の数を直接求めるスピゴット・アルゴリズム（英語版）を与える。これは π の十進法で n 桁目を計算するものではない[3]。そのような公式はプラウフによって2022年に発見されている[4]。BBPとBBPに触発されたアルゴリズムは、分散コンピューティングを使って π の多くの桁を計算するPiHex[5]などのプロジェクトで使用されている。BBPアルゴリズムの発見は驚くべきものであった。BBPアルゴリズムの発見以前、π のような超越数の n 桁目を直接計算するのは、最初の n 桁を計算するのと同程度に難しいと広く信じられていた[6]。

π 以外の無理数についてもBBP公式と類似の式が得られている。それらの類似の式はBBP型公式（BBP-typed formulae）[7]またはBBP系公式と呼ばれる。BBP型公式の一般形は以下のように表される：

α = ∑ k = 0 ∞ [ 1 b k p ( k ) q ( k ) ] {\displaystyle \alpha =\sum _{k=0}^{\infty }\left[{\frac {1}{b^{k}}}{\frac {p(k)}{q(k)}}\right]}
（BBP型公式）[8][注 1]

ここで α は無理数、p(k), q(k) は k の整数係数多項式、b は通常、位取り記数法の底であり 2 以上の整数を表すが、より一般には実数でもよい。また q は k の範囲でゼロにならず、p の次数は q より小さい。

ある数 α に対しBBP型公式を（すなわち p(k), q(k), b を）与える系統だったアルゴリズムは知られておらず、ある数に対するBBP型公式は今のところ実験的な方法によってのみ得られている。

特殊化

一般のBBP型公式の特殊な例として、多くの結果を与えたものに以下がある：

P ( s , b , m , A ) = ∑ k = 0 ∞ [ 1 b k ∑ j = 1 m a j ( m k + j ) s ] {\displaystyle 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]}
（BBP型公式の特殊化）

ここで s, b, m は整数、A = (a1, a2, ..., am) は整数の数列である。 関数 P はいくつかの解に対して簡潔な表記を与える。例えば、元のBBP公式は、s = 1, b = 16, m = 8, A = [4, 0, 0, −2, −1, −1, 0, 0] として以下のように書くことができる。

π = P ( 1 , 16 , 8 , [ 4 , 0 , 0 , − 2 , − 1 , − 1 , 0 , 0 ] ) {\displaystyle \pi =P\left(1,16,8,[4,0,0,-2,-1,-1,0,0]\right)}

既知のBBP型公式

BBP型公式のうちBBP公式より以前から知られているもので、BBP型公式の特殊化 P が簡潔な形になるものに、次のものがある。

ln ⁡ 10 9 = 1 10 + 1 200 + 1 3 000 + ⋯ = 1 10 ∑ k = 0 ∞ [ 1 10 k ( 1 k + 1 ) ] = 1 10 P ( 1 , 10 , 1 , [ 1 ] ) {\displaystyle {\begin{aligned}\ln {\frac {10}{9}}&={\frac {1}{10}}+{\frac {1}{200}}+{\frac {1}{3\,000}}+\cdots \\&={\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}}}
ln ⁡ 2 = 1 2 + 1 2 ⋅ 2 2 + 1 3 ⋅ 2 3 + ⋯ = 1 2 ∑ k = 0 ∞ [ 1 2 k ( 1 k + 1 ) ] = 1 2 P ( 1 , 2 , 1 , [ 1 ] ) . {\displaystyle {\begin{aligned}\ln 2&={\frac {1}{2}}+{\frac {1}{2\cdot 2^{2}}}+{\frac {1}{3\cdot 2^{3}}}+\cdots \\&={\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}}}

（実際、恒等式

ln ⁡ a a − 1 = ∑ k = 1 ∞ 1 a k ⋅ k {\displaystyle \ln {\frac {a}{a-1}}=\sum _{k=1}^{\infty }{\frac {1}{a^{k}\cdot k}}}

は a > 1 に対して成り立つ）

プラウフは、逆三角関数の冪級数にも触発された（P 表記は b が整数でない場合にも一般化できる）。

arctan ⁡ 1 b = 1 b − 1 b 3 3 + 1 b 5 5 − ⋯ = 1 b ∑ k = 0 ∞ [ 1 b 4 k ( 1 4 k + 1 + − b − 2 4 k + 3 ) ] = 1 b P ( 1 , b 4 , 4 , [ 1 , 0 , − b − 2 , 0 ] ) {\displaystyle {\begin{aligned}\arctan {\frac {1}{b}}&={\frac {1}{b}}-{\frac {1}{b^{3}3}}+{\frac {1}{b^{5}5}}-\cdots \\&={\frac {1}{b}}\sum _{k=0}^{\infty }\left[{\frac {1}{b^{4k}}}\left({\frac {1}{4k+1}}+{\frac {-b^{-2}}{4k+3}}\right)\right]\\&={\frac {1}{b}}P\left(1,b^{4},4,\left[1,0,-b^{-2},0\right]\right)\end{aligned}}}

新しい等式の探索

BBP型公式の特殊化 P を用いた方法で、π について既知の最も簡潔なものは、s = 1 かつ m > 1 のものである。b が指数 2 または 3、m が指数 2 または他の因数に富む値で、かつ数列 A の項のいくつかが 0 である場合、現在多くの公式が知られている。これらの公式の発見には、個々の和を計算した後、コンピュータでこのような線形結合を探索することが必要である。探索の手順は、s、b、mのパラメータ値の範囲を選び、その和を何桁にもわたって評価し、それらの中間和をよく知られた定数に、あるいはおそらくゼロに足し合わせる数列Aを、整数関係探索アルゴリズム（典型的にはHelaman FergusonのPSLQアルゴリズム）を用いて求めるというものである。

発見と導出

円周率のBBP公式の原型は、1995年にプラウフがPSLQアルゴリズムを用いて見つけたものである。また、BBP型公式の特殊化 P を用いて表現可能である。

π = ∑ k = 0 ∞ [ 1 16 k ( 4 8 k + 1 − 2 8 k + 4 − 1 8 k + 5 − 1 8 k + 6 ) ] = P ( 1 , 16 , 8 , ( 4 , 0 , 0 , − 2 , − 1 , − 1 , 0 , 0 ) ) . {\displaystyle {\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}}}

上式は、以下の等価な2多項式の比に還元される。

π = ∑ k = 0 ∞ [ 1 16 k ( 120 k 2 + 151 k + 47 512 k 4 + 1024 k 3 + 712 k 2 + 194 k + 15 ) ] {\displaystyle \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]}

この式は、かなり単純な証明によって、π に等しいことが示されている[9]。

BBPアルゴリズム

以下ではBBP公式から円周率 π の十六進法で n 桁目の数を与えるアルゴリズム（BBPアルゴリズム）について説明する。

一般に数 α の b 進法で小数点以下 n 桁目の値を求めるには、α に bn−1 を掛けて（n − 1 桁左にシフトして）その小数部分の最上位の桁を計算すればよい。BBPアルゴリズムでは n 桁目だけを必要とするため、小数部の計算は通常の浮動小数点数の演算で済む。bn−1α の整数部分の計算を省略できるなら、α の小数点以下 n − 1 桁までの値を求めずに n 桁目を計算できる。

bn−1α の整数部の計算の省略は剰余算することでできる。 α がBBP型公式の一般形で書ける場合、分母 q(k) で分子 b(n−1)−kp(k) を割った余りだけを計算するよう変形し、最後に和の小数部を取り出す：

{ b n − 1 α } mod 1 = { ∑ k = 0 ∞ b ( n − 1 ) − k p ( k ) mod q ( k ) q ( k ) } mod 1 {\displaystyle \left\{b^{n-1}\alpha \right\}{\bmod {1}}=\left\{\sum _{k=0}^{\infty }{\frac {b^{(n-1)-k}p(k){\bmod {q}}(k)}{q(k)}}\right\}{\bmod {1}}}

特にBBP型公式の特殊化の s = 1 の場合では、以下のように書ける：

{ b n − 1 P ( 1 , m , b , A ) } mod 1 = { ∑ j = 1 m a j [ ∑ k = 0 ∞ b ( n − 1 ) − k mod ( m k + j ) m k + j ] } mod 1 {\displaystyle \left\{b^{n-1}P(1,m,b,A)\right\}{\bmod {1}}=\left\{\sum _{j=1}^{m}a_{j}\left[\sum _{k=0}^{\infty }{\frac {b^{(n-1)-k}{\bmod {(}}mk+j)}{mk+j}}\right]\right\}{\bmod {1}}}

ここで x mod y は x を y で割った余りを表す[注 2]。和の中の剰余算が必要となるのは分子が分母より大きくなる場合のみであり、分子が小さい場合には必要ない。そのため、上述の bn−1P の例では、k ≤ n − 1 となる最初の項だけ剰余算を行う。

上述の手順に従い π の16進法での小数点以下 n 桁目の計算を以下に示す。まずBBP公式を以下のように変形する：

16 n − 1 π = 4 ∑ k = 0 ∞ 16 ( n − 1 ) − k 8 k + 1 − 2 ∑ k = 0 ∞ 16 ( n − 1 ) − k 8 k + 4 − ∑ k = 0 ∞ 16 ( n − 1 ) − k 8 k + 5 − ∑ k = 0 ∞ 16 ( n − 1 ) − k 8 k + 6 {\displaystyle 16^{n-1}\pi =4\sum _{k=0}^{\infty }{\frac {16^{(n-1)-k}}{8k+1}}-2\sum _{k=0}^{\infty }{\frac {16^{(n-1)-k}}{8k+4}}-\sum _{k=0}^{\infty }{\frac {16^{(n-1)-k}}{8k+5}}-\sum _{k=0}^{\infty }{\frac {16^{(n-1)-k}}{8k+6}}}

それぞれの和の小数部の最初の数桁を計算し、次に右辺全体の小数部を計算する。最終的な結果から小数第一位を取り出せば π の小数第 n 位が求まる。

例えば最初の和の小数部は以下のようになる：

{ ∑ k = 0 ∞ 16 ( n − 1 ) − k 8 k + 1 } mod 1 = { ∑ k = 0 n − 1 16 ( n − 1 ) − k mod ( 8 k + 1 ) 8 k + 1 + ∑ k = n ∞ 16 ( n − 1 ) − k 8 k + 1 } mod 1 {\displaystyle \left\{\sum _{k=0}^{\infty }{\frac {16^{(n-1)-k}}{8k+1}}\right\}{\bmod {1}}=\left\{\sum _{k=0}^{n-1}{\frac {16^{(n-1)-k}{\bmod {(}}8k+1)}{8k+1}}+\sum _{k=n}^{\infty }{\frac {16^{(n-1)-k}}{8k+1}}\right\}{\bmod {1}}}

前述の通り、右辺の2つ目の和に関しては適当な精度で値が収束したと見なされるまで計算すればよいが、1つ目の和は n 回の冪剰余[注 3]が必要となる。冪剰余の計算のビット計算量は、M(j) を j ビット整数の乗算に必要な計算量として、O(log n ⋅ M(log n)) で行える[注 4][注 5]。従って全体の計算量は O(n ⋅ log n ⋅ M(log n)) となる。乗算の計算量は乗算アルゴリズムに依存するが、筆算と同じ方法なら M(j) = O(j2) であり、この場合のBBPアルゴリズムのビット計算量は O(n ⋅ (log n)3)) となる。乗算をショーンハーゲ・ストラッセン法で行えばビット計算量は M(j) = O(j ⋅ log j ⋅ log(log j)) となり、この場合のBBPアルゴリズムのビット計算量は以下の通り：

O ( n ⋅ ( log ⁡ n ) 2 ⋅ log ⁡ ( log ⁡ n ) ⋅ log ⁡ ( log ⁡ ( log ⁡ n ) ) ) {\displaystyle O(n\cdot (\log n)^{2}\cdot \log(\log n)\cdot \log(\log(\log n)))}

これは（ショーンハーゲ・ストラッセン法で乗算を実装した場合の）π を n 桁計算するアルゴリズムよりはわずかに遅く、BBPアルゴリズムのビット計算量は log(log(log n)) 倍だけ大きい。

BBPアルゴリズムの結果は n 桁目以降の値も含んでおり、n − 1 桁目ないし n + 1 桁目の計算結果と比較することで計算結果が正しいことを検証できる。

BBPと他の円周率計算方法との比較

このアルゴリズムは、数千、数百万の桁を持つカスタムデータ型を必要とせずに π を計算する。この方法は、最初の n − 1 桁を計算せずにn桁目を計算し、小さくて効率的なデータ型を使用することができる。

BBPは π の任意の桁の値を直接計算することができ、その間のすべての桁を計算しなければならない数式よりも少ない計算量で計算できるが、BBPは線形的（ O ( n ( log ⁡ n ) O ( 1 ) ) {\displaystyle O(n(\log n)^{O(1)})} [10]）であり、n の値が大きくなるほど計算に要する時間が長くなる。つまり、ある桁が「先に」あるほどBBPの計算に要する時間は長くなり、標準の π 計算アルゴリズムと同じようになっていく[11]。

一般化

D. J. Broadhurst はBBPアルゴリズムを一般化し、他の多くの定数をほぼ線形時間および対数空間で計算するために使用できるようにした[12]。カタランの定数、π3、π4、アペリーの定数 ζ(3)、ζ(5)（ζ(x) はリーマンゼータ函数）、log32、log42、log52 および π と log2 の冪乗による様々な積について明示的に結果を与えている。これらの結果は、主に多重対数梯子を用いて得られるものである。

関連項目

出典

[脚注の使い方]
  1. ^ Bailey, Borwein & Plouffe 1997, 1. Introduction, Theorem 1, (1.2).
  2. ^ プラウフのサイト.
  3. ^ Gourdon & Sebah 2003.
  4. ^ Weisstein, Eric W. "Digit-Extraction Algorithm". mathworld.wolfram.com (英語).
  5. ^ “PiHex Credits”. Centre for Experimental and Constructive Mathematics. Simon Fraser University (1999年3月21日). 2017年6月10日時点のオリジナルよりアーカイブ。2018年3月30日閲覧。
  6. ^ Bailey, Borwein & Plouffe 1997, 1. Introduction.
  7. ^ Weisstein, Eric W. "BBP Formula". mathworld.wolfram.com (英語).
  8. ^ Bailey, Borwein & Plouffe 1997, 3. The Algorithm.
  9. ^ Bailey et al. 1997.
  10. ^ Bailey, Borwein & Plouffe 1997.
  11. ^ Bailey, David H. (2006年9月8日). “The BBP Algorithm for Pi”. 2013年1月17日閲覧。 “Run times for the BBP algorithm ... increase roughly linearly with the position d.”
  12. ^ Broadhurst 1998.

注釈

[脚注の使い方]
  1. ^ Bailey, Borwein & Plouffe 1997 では b の指数は ck となっていて、正整数 c が付け加わっている。
  2. ^ この mod 演算の対象は整数に限らず、x − [x / y] y の意味で使う。ここで [⋅] は引数の整数部を取り出す関数である。文献によっては引数の小数部を表す関数を {⋅} と定義し、mod の代わりに使用している（x mod y = {x / y}y）。
  3. ^ 浮動小数点数の割り算と足し算も行うが、冪剰余の計算が支配的となる。
  4. ^ 対数の底の異なりは定数としてしか寄与しないため、オーダーの計算では底を書かない。
  5. ^ M の引数が log n となっているのは、除数 q に対する冪剰余の計算には q2 を表せる長さを持つ整数型が必要となるためである。q2 のビット長は 1 + 2log2 q であり、これは漸近的に log q に等しい。BBPアルゴリズムによる π の計算では除数は高々 n のオーダーになるから、q ≈ n と考えてよく、ビット長は log n 程度と見積もれる。π 以外の例では q の取りうる大きさは異なり、例えば q ≈ n2 などになり得る。ビット計算量の見積もりには顕れないが、ナイーブな方法では浮動小数点数の精度が不足し得る。

参考文献

  • Bailey, David H.; Borwein, Peter B.; Plouffe, Simon (1997). “On the Rapid Computation of Various Polylogarithmic Constants”. Mathematics of Computation 66 (218): 903–913. doi:10.1090/S0025-5718-97-00856-9. MR1415794. 
  • Bailey, David H.; Borwein, Jonathan M.; Borwein, Peter B.; Plouffe, Simon (1997). “The quest for pi”. Mathematical Intelligencer 19 (1): 50–57. doi:10.1007/BF03024340. MR1439159. 
  • Gourdon, Xavier; Sebah, Pascal (12 February 2003). "N-th Digit Computation". numbers.computation.free.fr. 2023年1月22日閲覧。
  • Broadhurst, D. J. (1998). Polylogarithmic ladders, hypergeometric series and the ten millionth digits of ζ(3) and ζ(5). arXiv:math.CA/9803067. 

外部リンク

  • リチャード・リプトン（英語版）, "Making An Algorithm An Algorithm — BBP", weblog post, July 14, 2010.
  • リチャード・リプトン, "Cook’s Class Contains Pi", weblog post, March 15, 2009.
  • “A compendium of BBP-type formulas for mathematical constants, updated 15 Aug 2017”. 2019年3月31日閲覧。
  • デイヴィッド・H・ベイリー（英語版）, "BBP Code Directory", web page with links to Bailey's code implementing the BBP algorithm, September 8, 2006.