ベイリー＝ボールウェイン＝プラウフの公式 （ベイリー＝ボールウェイン＝プラウフのこうしき、英 : 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公式 ）名称はBBP公式に関する論文 (Bailey, Borwein & Plouffe 1997 ) の著者らデイヴィッド・H・ベイリー （英語版 ） ピーター・ボールウェイン （英語版 ） [ 2]
BBP公式は、先行する桁を計算せずに π の十六進法 で n 桁目（二進法 で 4n 桁目）の数を直接求めるスピゴット・アルゴリズム （英語版 ） π の十進法 で n 桁目を計算するものではない。そのような公式はプラウフによって2022年 に発見されている[ 4] アルゴリズム は、分散コンピューティングを使って π の多くの桁を計算するPiHex[ 5] π のような超越数 の n 桁目を直接計算するのは、最初の n 桁を計算するのと同程度に難しいと広く信じられていた。
π 以外の無理数 についても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型公式 ）[ 注 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 = (a 1 , a 2 , ..., am )数列 である。
関数 P はいくつかの解に対して簡潔な表記を与える。例えば、元のBBP公式 は、s = 1b = 16m = 8A = [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型公式の特殊化 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 = 1m > 1b が指数 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]}
この式は、かなり単純な証明によって、π に等しいことが示されている。
以下ではBBP公式 から円周率 π の十六進法 で n 桁目の数を与えるアルゴリズム （BBPアルゴリズム）について説明する。
一般に数 α の b 進法で小数点以下 n 桁目の値を求めるには、α に b n −1n − 1n 桁目だけを必要とするため、小数部の計算は通常の浮動小数点数 の演算で済む。b n −1α α の小数点以下 n − 1n 桁目を計算できる。
b n −1α 剰余 算することでできる。
α がBBP型公式 の一般形で書ける場合、分母 q (k )b (n −1)−k p (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] b n −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 (j 2 )O (n ⋅ (log n )3 ))ショーンハーゲ・ストラッセン法 で行えば乗算の計算量は M (j ) = O (j ⋅ log j ⋅ log(log j ))
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 − 1n + 1
このアルゴリズムは、数千、数百万の桁を持つカスタムデータ型を必要とせずに π を計算する。この方法は、最初の n − 1 桁を計算せずにn桁目を計算し、小さくて効率的なデータ型を使用することができる。
BBPは π の任意の桁の値を直接計算することができ、その間のすべての桁を計算しなければならない数式よりも少ない計算量で計算できるが、BBPは線形的（
O
(
n
(
log
n
)
O
(
1
)
)
{\displaystyle O(n(\log n)^{O(1)})}
n の値が大きくなるほど計算に要する時間が長くなる。つまり、ある桁が「先に」あるほどBBPの計算に要する時間は長くなり、標準の π 計算アルゴリズムと同じようになっていく[ 11]
D. J. Broadhurst はBBPアルゴリズムを一般化し、他の多くの定数をほぼ線形時間および対数空間で計算するために使用できるようにした。カタランの定数 、π3 、π4 、アペリーの定数 ζ(3) 、ζ(5) （ζ(x ) はリーマンゼータ函数 ）、log3 2 、log4 2 、log5 2 および π と log2 の冪乗による様々な積について明示的に結果を与えている。これらの結果は、主に多重対数梯子 を用いて得られるものである。
^ プラウフのサイト .^ ^ “PiHex Credits ”. Centre for Experimental and Constructive Mathematics . Simon Fraser University (1999年3月21日). 2017年6月10日時点のオリジナルよりアーカイブ 。2018年3月30日閲覧。 ^ Weisstein, Eric W. “BBP Formula” . mathworld.wolfram.com (英語). ^ 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 .”
^ Bailey, Borwein & Plouffe 1997 では b の指数は ck となっていて、正整数 c が付け加わっている。^ この mod 演算の対象は整数に限らず、x − [x / y ] y [⋅] は引数の整数部を取り出す関数である。文献によっては引数の小数部を表す関数を {⋅} と定義し、mod の代わりに使用している（x mod y = {x / y }y
^ 浮動小数点数の割り算と足し算も行うが、冪剰余 の計算が支配的となる。
^ 対数の底の異なりは定数としてしか寄与しないため、オーダーの計算では底を書かない。
^ M の引数が log n となっているのは、除数 q に対する冪剰余の計算には q 2 q 2 1 + 2log2 q であり、これは漸近的に log q に等しい。BBPアルゴリズムによる π の計算では除数は高々 n のオーダーになるから、q ≈ n log n 程度と見積もれる。π 以外の例では q の取り得る大きさは異なり、例えば q ≈ n 2
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 . MR 1415794 . 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 . MR 1439159 . Gourdon, Xavier; Sebah, Pascal (12 February 2003). “N-th Digit Computation” . numbers.computation.free.fr . 2023年1月22日閲覧 . {{cite web2 }}: CS1メンテナンス: デフォルトと同じref (カテゴリ )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.![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/af6bc360851499dd2ab2a90bee03fbe2040089d5)
![{\displaystyle \alpha =\sum _{k=0}^{\infty }\left[{\frac {1}{b^{k}}}{\frac {p(k)}{q(k)}}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fa85e25344327bcebd7461413dc9e28dc1a060ce)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/79b71ecf9ae82b1b067d955f13f3c46dc89bb0ed)
![{\displaystyle \pi =P\left(1,16,8,[4,0,0,-2,-1,-1,0,0]\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d5063fb0ed6ff69fe613a4f467d7976346a030c0)
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e42438a4dabc61cbb71c437fff81fc262f3eaac5)
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/61f336c002ef427a75d224981716974747bc62e6)
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ede4c3bf7e9da0c8444c4948c4bd5e8c8f64b815)
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c9b4d1e7c86fa6adcc98b2b80dbbd27e7e43e1b3)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ba24829fa8716fb94e690df022bb7a9e3e48ed23)
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f628c092b1361c97f8f0ad294a1b76441bc6943a)
