円周率 (16進数) の8304820223桁目から16桁を求めたはなし
これも割とぷちネタ。
タイトル通り、円周率の途中だけ求めるBBPの公式ってのを実装して実行してみた。BBPの公式を利用したアルゴリズムは凄い。少ない一定量のメモリで (時間は桁数に比例するが) 円周率を途中から求められる。結構面白くてハマったので自分記録。2022年12月頭の (足掛け) 3日間くらいの記録。
- 本番ソース (Github): https://github.com/tarohs/pi (bbp.py)
途中から失礼します〜
円周率を10進で100億桁、自分で書いたプログラムで求めてみたが、2進数表現が欲しくなり100億桁の10進数→332億桁の2進数変換をした。だがその末尾数ビットか数バイトは正しい自信が持てなかった。
いろいろやっているときに、円周率を途中だけ求める『BBPの公式』、というのがあるのに気がついた。知らなかったのもそのはず、わたしが学生時代には誰も知らなかった公式である。円周率の長い歴史の中ではごく最近の1995年に、BさんとBさんとPさんが気づいた。
ただしこれは、16進数 (2進数) で任意の桁を求める、という方法である。わしが (10進数) 100億桁ぽっち求めて喜んでいるのは、あわよくば紙かなんかに100億の数字を出力してやれ、とのアピールを目論んでいるからであって、2進数の01や16進数で『0x3.243f6a88...』とか出ても嬉しくもなんともない。なかった。
ちなみに、10進数の
まあそれはともかく、10進 (100億桁) →2進 (33219280957ビット) 変換した円周率を検算したくなったので、BBPアルゴリズムを実装してみることにした。
*BBP系公式 (円周率.jp)
*BBPアルゴリズム
*円周率の小数点以下8000兆桁めをGeForceで求める方法～GTC 2013の「π」セッションレポート (4gamers.net)
世の中にはこんなに親切に解説してくださっているページもあるのだが、どーも頭に入って来ない。つか、総じて他のアルゴリズム (arctan公式やモンテカルロ法からガウス--ルジャンドル、Chudnovskyまで) に比べて、BBPには「やってみました」「こうすればいいですよ」系のページが少ない。
世の、円周率何十兆桁とか求めているいわゆる『円周率ガチ勢』におかれてましては、この公式で「ささっと検算して合ってました」的な使い方をされているので (もちろん途中の何京桁目だけ求めるガチ勢もいたり)、そ、そんなにみんな簡単にできるの!? と引け目を感じた。
ところが上の3番めのEd KarrelsさんのGPUを用いた計算方法のプレゼンテーションのスライドを何枚かみていたら、「あ!! そういうことか!!」とやっとわかった (ちなみにきちんと説明してくれている上の1番目・2番めのページでわしが混乱していたのは、小数点以下の桁数の取り方とか、ノーテーションとか高速指数計算 (後述) の引数に何を突っ込むかとか、であった)。
しくみ
ひとの褌で偉そうに解説しようっと(笑)。
BBP公式は、ふつーに
ところが、BBP公式をあれこれすると、途中の計算結果は劇的に少ない桁数を保持すれば十分、になる。VMが たった50Gしかない (実メモリ16GB) のMacbook Proでも安心である。漸化式の中に
どーいうことか。上述のEd Karrelsさんの言葉を借りると「小数点以下
この式の右辺は
ここで
そうすると、
ここに高速指数計算法 (とわたしは覚えているが、いろいろな名前があるようだ) という、『指数計算をしつつ割り算をしつつ整数部分を無視する』という都合の良いアルゴリズムがある。いや正確に言うと『指数計算を
上の
これが納得できない場合は、10進数ですげー大きな数、たとえば100億 （※ 小学生かよ） を3で割った場合を考えてみるとよい。(100億 ÷ 3) の小数部分、つまり『3333333333.33333...』の小数部分は、 (1 ÷ 3) の小数部分と同じで『.33333...』だ。なぜなら、100億と1は、3を法として合同だから……簡単にいうと、1は100億を3で割った余りだから、である。
問題は『0b101だったら
BBPアルゴリズムをざっくり言うと
-
\pi 16^{83億桁} - だがその計算では、整数部分は無視できる
- だから83億回くらいの繰り返しで、小数部分だけを積算してやる
- でもやっぱ、83億回くらいの繰り返しが必要だ
ということである。
基本的な実装
さて、ここまでわかったら、あとはPython3で書くだけである。高速化するなら後でgmpy2を使えばいいので、とりあえずPythonのintで書いて、小さな数で試してみよう。以前と同じで、割り算の部分は固定小数点 (整数演算) で誤差を避けて……と。
基本的な実装 （プログラムリスト）
def bbp(digstart):
digstart = digstart - 1
overdig = 13
rss = 0.
for k in range(0, digstart):
rs = 4. + exp16mod(digstart - k, 8 * k + 1) * 4 / (8 * k + 1) -\
exp16mod(digstart - k, 8 * k + 4) * 2 / (8 * k + 4) -\
exp16mod(digstart - k, 8 * k + 5) * 1 / (8 * k + 5) -\
exp16mod(digstart - k, 8 * k + 6) * 1 / (8 * k + 6)
k8 = k << 3
rs = 4. + (pow(16, digstart - k, k8 + 1) << 2) / (k8 + 1) -\
(pow(16, digstart - k, k8 + 4) << 1) / (k8 + 4) -\
(pow(16, digstart - k, k8 + 5) ) / (k8 + 5) -\
(pow(16, digstart - k, k8 + 6) ) / (k8 + 6)
rss += rs
rss = rss - int(rss)
rss, _ = math.modf(rss)
# print(k, ':', rss)
#print('---')
for k in range(digstart, digstart + overdig):
rs = 16 ** (digstart - k) * (4 / (8 * k + 1) -\
2 / (8 * k + 4) -\
1 / (8 * k + 5) -\
1 / (8 * k + 6))
rss = rss - int(rss)
rss, _ = math.modf(rss)
# print(k, ':', rss)
return rss
出来てはみたものの、デバッグどうしよう? 状態である。まあ16進数のhex()でprint()してみればよい。
そのへんで気づいたのだが、今回gmpy2の出番はない (※ 実は後で蘇る)。前にgmpy2を使ったのは整数の加減乗除が速いからであって、まあその部分を固定小数点で書いても、今回は大した桁数でない。高速指数計算法が桁数をびしばし削ってくれるので、ほんとーに出てくる整数といったらデバッグ段階では10進数桁 (10000桁程度なら4桁とか8桁とか) くらいの数しかない。実数部分はさすがに10000回とか積算するのである程度の精度は重要だが、Pythonのfloatで十分である。固定小数点は (頭が混乱するひとつの原因なので) すぐにやめた。
試みに数回のループ計算の生の値をみてみる。小数点シフト量0にして、1回だけでループ計算してみる。いきなり『3.133333』が得られて (0.13は0x3.243f...の2と近い) 「おぉ」とか呟いてしまう。
0x3.243f... (16進数) に対して
- 0回シフト: 小数部 = .1333... = 2.133.../16
- 1回シフト: 小数部 = .2664... = 4.262.../16
- 2回シフト: 小数部 = .2463... = 3.941.../16
- 3回シフト: 小数部 = .9627... = 15.40.../16
っていうか、
そして、『ふつーでない』積算が終わった段階では、小数点以下
欲しい
ちなみにある桁 (16進) が繰り上がるかどうかは、単純に考えて確率
ここで『任意の』余計な桁数を『ふつーに』（浮動小数点で） 計算する部分を書こうと思ったので、なんだかんだ丸1日ハマってしまった。面倒になったので『任意の』はやめて、Pythonのfloatの桁数の有効精度だけ (正確にはときどき整数部が大きくなるので、それマイナス10進数桁) 余計に計算して、使えるところだけ表示しよう、というアプローチに出た。
実際には (後でわかったが) Pythonのfloat型は52ビットの精度がある。つまり16進なら13桁。16進10000桁の表で、100桁、1000桁、10000桁などで答え合わせをしてみると、9〜10桁くらいは信用できる値である。それに合わせてとりあえず、表示を8桁にした。どこで4桁くらい失っているかというと、おそらく10000回 (16進で4桁回) の累積誤差が4桁未満、あと整数部分に1桁程度はみ出してしまう (すぐに切り捨てて大きな数にならないようにするが)、それらのせいなのだろう。
基本の実装では、桁数を増やすとほぼリニアに時間が増える (桁数10倍で時間11倍程度)。1000万桁 (16進数) で1分程度なので、計算してみると目的の83億桁くらいではほぼ丸1日弱、ということになる。もう少し速くしてみてから本番に入ろう （ぜったいそれまで丸1日以上掛かるけどねっ）。
- 計算時間については後述『本番』の節参照。
自前の関数と組み込み関数と
前記の高速指数計算法の関数は、まあわしは生涯で10回くらいは書いたことがあるので11回目もさくっと書いた。
高速指数計算法 (16のa乗 mod b) (プログラムリスト）
def exp16mod(a, b):
p = 16
r = 1
while 0 < a:
if a & 1 == 1:
r = (r * p) % b
p = (p * p) % b
a >>= 1
return r
だがいろいろ調べているうちにPythonの基本の組み込みpow()関数に3引数を与えると (pow(a:int, b:int, c:int))、pow()関数は避けましょう」な発想になっていたのだが、実は高速指数計算が組み込みで実装されていたとは。さっそく、自前の関数をクビにしてこちらを使ってみると、かなり速かった (当たり前である。ループをPythonで書いているのがアセンブラになるんだから)。
あと小数部分を取り出すのにa - int(a)みたいにしていたが、これもmathの中のfrac, _ = math.modf() (タプルで頼んでない整数部も返るので、捨てる) を使ってみたら、ほんのり速くなったので、そちらを使ってみた。
gmpy2の浮動小数点は、速いの?
ここまできたら、gmpy2のmpfr型を使ったら、速くなるの? と思った。ちなみに以前、円周率を求めたときにmpz型がPythonのint型より速かったのは、計算が数億桁とかの整数 (固定小数点実数) だったからで、今回は前述のように整数部分の桁数はほぼない。ふつーのひとが日常で行なう数桁〜 数100桁 (ぉぃ) 未満の計算なら、GMPの乗算のFFT (NTT)法どころか、Karatsuba法 (まあPythonのint乗算はこれらしいが) でさえ通常の乗算よりオーバーヘッドで遅くなってしまう。
まずはmpfr型の実力拝見。あちこちの初期値をmpfr()で囲んでやるだけでgmpy2化できるので、やってみた……ら、遅くなった(泣)。大差ないのだが1〜2割くらい、Pythonのint型より遅い。
ちなみに、gmpy2では高速指数計算はpowmod()、小数部分のみを取り出すのはfrac()という関数がある。
後でmpfr型は復活するので、『1〜2割でよかったね♡』というところである。
Bellardのアルゴリズム
なにやら、ガチ組の途中桁数計算では『Bellardの公式』を使っているらしい。オリジナルのBBPが16進数向けであるのに対し、この公式ではループ1回の項が進むときの係数のオーダとして
すでにBBPは正しく動いているので、コピーして書き換えてみる。頭がこんがらがってきたのは、BBPの公式の4項に対して7項もある……だけではなく、ループ1回毎に
Bellardの公式 （プログラムリスト）
def bellard(digstart):
digstart = int((digstart - 3) / 5) * 2
print('dig', digstart * 2.5 + 3, '... (digstart', digstart, ')')
overdig = 6
rss = 0.
for k in range(0, digstart):
rs = \
- pow(1024, digstart - k, 4 * k + 1) * 32 / (4 * k + 1) \
- pow(1024, digstart - k, 4 * k + 3) / (4 * k + 3) \
+ pow(1024, digstart - k, 10 * k + 1) * 256 / (10 * k + 1) \
- pow(1024, digstart - k, 10 * k + 3) * 64 / (10 * k + 3) \
- pow(1024, digstart - k, 10 * k + 5) * 4 / (10 * k + 5) \
- pow(1024, digstart - k, 10 * k + 7) * 4 / (10 * k + 7) \
+ pow(1024, digstart - k, 10 * k + 9) / (10 * k + 9)
if k % 2 == 0:
rss += rs
else:
rss -= rs
rss, _ = math.modf(rss)
for k in range(digstart, digstart + overdig):
rs = (1024 ** (digstart - k)) * (
-32 / (4 * k + 1) \
- 1 / (4 * k + 3) \
+ 256 / (10 * k + 1) \
- 64 / (10 * k + 3) \
- 4 / (10 * k + 5) \
- 4 / (10 * k + 7) \
+ 1 / (10 * k + 9))
if k % 2 == 0:
rss += rs
else:
rss -= rs
rss, _ = math.modf(rss)
rss, _ = math.modf(1. + rss)
return rss
Bellardの公式では明示的にマイナスが出てくる。modf()では (gmpy2.frac()も) 結果が負数になったら小数部分としてマイナスの値を返す。そのまま計算していると、わけがわからなくなってしまう。BBPのときも公式の項から項を引くとマイナスになったりして (いやそのままではならんが、小数部だけになると大小関係がどうなるかわからないので)、負数に整数部があまりに大きくならないような適当な数を足してからmodf()したりしていた。『各 (公式の) 項の最大値』『ループ1回に積算される値の最大値』とかいろいろ悩みながらやっていた。
Bellardの公式でいろいろやってみた結果、負数 (の小数部) については、それほど神経質に考えなくてよいようだ。そのまま最後まで計算して、最終の積算結果がマイナスだったら1を足せばよいだけである。整数部分があまりに大きく (またはマイナスに大きく) なると、限られた浮動小数点の精度の中では下の桁に悪影響があるが、まあまあmodf()するとかの配慮は要らないだろう。
さて話が逸れたが。Bellardのアルゴリズムはそれ以外にも胃に悪い(笑)。
その上、
それでも10倍速いとかなら我慢して赤く塗るが、BBPの2〜3割程度しか速くならない。収束は速いが、それ以上に各項の計算が複雑化しているのが災いしているようだ。
- 計算速度のグラフは『本番』の節参照。
欲しい桁数を自由に指定できないのと、そうすると同じ精度では欲しい範囲の正しい桁数が減っちゃうことから、やっぱりBBPに戻ることにした。
N
ここまではPytho3のint型の演算精度が52ビットなので、16進数8桁程度は『信頼できる』として値を表示していた。さて本番実行が近くなると、やっぱり連続した16進16桁とかくらいは欲しくなる。
つまり、たかだか10000桁 （16進数） くらいで答え合わせをしている段階では、10000回くらいしかループは回っていないので、除算のエラーが積算されても最悪14ビット未満くらいだろう。積算が整数部にはみ出すことがあるので4〜8ビットの損失があるとすると、52ビット精度から引き算すると、8桁の16進数 （4ビット） 32ビットがやっとである。本番では83億回くらい積算するので、積算誤差は32ビットくらいか。もっと精度が欲しくなる。
gmpy2では、浮動小数点の精度が指定できる。ただし、プログラム全体として （変数ごと、ではない） 精度をビット単位で指定できる、ということである。後にここでまたハマることになる。2度ボツにしたgmpy2を、またまた復活してきた。
結局、表示して欲しいビット数 + 求めたい桁を表わすビット数 + 20ビットくらいの余裕をみた精度に設定することにした。
- リストは本番プログラム （後出） 参照
さて並列化
もちろん以前の経験から、さいしょからこいつはマルチプロセス・マルチコアで並列化で速くしてやろう、と目論んでいた。mainもはじめから行儀よく書いてある（笑）。プログラムを作って動かしてみているMacbook Proは文房具ではあるが10コアある。足りなければ、AWSでまた64コアとか借りてやってもよい(笑)。
BBPアルゴリズムで並列化のポイントはふたつある。ひとつは、公式の中の項 （4つある） ごとの計算である……これは4 （Bellardでは7） 以上に分けられない。もうひとつは、ループ回数 (
なんにせよ、『並列化が簡単』すなわち『前の結果に依存しない』漸化式になっている、というのはありがたい。
もうgmpy2 (mpfr精度を指定できる） バージョンまで出来ていたので、それをまずは簡単な上記前者 （公式の項） ごとに分けてみた。方法は、以前使ったのと同じconcurrent.futuresのProcessPoolExecutorを使う方法である。
- リストは本番プログラム （後出） 参照
あれ? なんか結果がおかしい。小さい
といっても、なぜ計算を分割したときだけ計算誤差が発生しているのかわからない。分割点がおかしい (ループ項数に重複とか抜けがある） とか、分割バージョンでmpfr()を宣言し忘れているとかPythonのfloatで計算してからmpfr型に代入しちゃって精度を落としているとか、親子プロセス間のデータの受け渡し (ピクルス化するらしい) とか、いろいろな疑いでプログラムを見直した。
せめて、分割されたそれぞれの項の和の正解がわかれば……と思ったので、まずはループの項数
比較してみると、ある桁数以下に誤差がある……やっぱりPythonのfloatで計算している?
おやっと思って、マルチプロセスで呼び出されたサブルーチン内でのmpfrの精度を表示させてみると……52ビット （初期値） に戻っている。つまり、mainで計算桁数に応じたビット数 （90ビットくらい） まで精度を上げているのに、子供が起動されるとその中では52ビット計算になっているわけだ。
なるほど……concurrent.futuresの子プロセスの起動の方法を思い出したら、forkした子供がスクリプトを読んで関数を実行する （っぽい） のだった。mainは読み飛ばす。いくらグローバルな属性であっても、子供は精度の設定をしていない。そうか、では子プロセスに入ったところでもう一度精度を設定してやらなければならない。
引数にそれを追加して設定……こんどは、プロセス分割してもエラーがない。公式の4項で分割しても、ループ回数で分割してもエラーがない。そこで、両方で分割して8コア （決め打ち） にすることにした。どーせ10コアMacbook Proで実行しても推定1時間ちょっとである、64コアとかのAWSを借りるまでもない（笑）。
公式の4項それぞれで前半・後半に分割で8つ……前後半の最適な分割点はどこ （欲しい
何通りか試してみた。プロセスは#0〜#6 (公式の項1〜4のループ前半、1〜3の後半） を起動してから、親は#7 （公式の項4のループ後半） を実行する。終わった後、#0〜6の順に子供を迎えに行くので
- 迎えに行ってもしばらく子供が帰ってこなければ、前半が重い
- 親が迎えに行ったら全員揃っているようであれば、後半が重い
ことになる。
試してみると、0.5くらいのときにほぼ親子が同時に終わることがわかった。つまり……桁数が増えても、あまり計算時間は変わっていない。
本番プログラム （プログラムリスト） （8コア決め打ち並列化・gmpy2使用・mainも = 名前が『bbpmpgmp()』ってw）
def bbpmpth(startk, endk, digstart, digdisp, shift, mod8):
fpprec = len(str(digstart)) + 4 * digdisp + 20
gmpy2.get_context().precision = fpprec
rs = mpfr(0.)
for k in range(startk, endk):
k8 = k << 3
dk = mpz(digstart - k)
rs += mpfr(gmpy2.powmod(mpz(16), dk, k8 + mod8) << shift) / (k8 + mod8)
rs = gmpy2.frac(rs)
return rs
def bbpmpgmp(digstart, digdisp):
digstart = digstart - 1
overdig = digdisp + 10
rss = mpfr(0.)
digdiv = int(digstart * .5)
futures = []
pm = [1, -1, -1, -1, 1, -1, -1, -1]
with Executor() as exec:
futures.append(exec.submit(bbpmpth, 0, digdiv, digstart, digdisp, 2, 1))
futures.append(exec.submit(bbpmpth, 0, digdiv, digstart, digdisp, 1, 4))
futures.append(exec.submit(bbpmpth, 0, digdiv, digstart, digdisp, 0, 5))
futures.append(exec.submit(bbpmpth, 0, digdiv, digstart, digdisp, 0, 6))
futures.append(exec.submit(bbpmpth, digdiv, digstart, digstart,
digdisp, 2, 1))
futures.append(exec.submit(bbpmpth, digdiv, digstart, digstart,
digdisp, 1, 4))
futures.append(exec.submit(bbpmpth, digdiv, digstart, digstart,
digdisp, 0, 5))
# futures.append(exec.submit(bbpmpth, digdiv, digstart, digstart,
# digdisp, 0, 6))
print('start future #0-', len(futures))
rs = bbpmpth(digdiv, digstart, digstart, digdisp, 0, 6)
print('done future #', len(futures))
rss = mpfr(8.) - rs
for i in range(len(futures)):
rss += pm[i] * futures[i].result()
print('done future #', i)
rss = gmpy2.frac(rss)
for k in range(digstart, digstart + overdig):
rs = 16 ** (digstart - k) * (mpfr(4.) / (8 * k + 1) -\
mpfr(2.) / (8 * k + 4) -\
mpfr(1.) / (8 * k + 5) -\
mpfr(1.) / (8 * k + 6))
rss += rs
rss = gmpy2.frac(rss)
return rss
if __name__ == '__main__':
global fpprec
if len(sys.argv) != 2:
print('usage: bbp DIGSTART', file = sys.stderr)
sys.exit(1)
digstart = int(sys.argv[1])
digdisp = 16
# precision requires at least
# (displayed_precisions + loop_error + redundant)
fpprec = len(str(digstart)) + 4 * digdisp + 20
gmpy2.get_context().precision = fpprec
print('prec:', gmpy2.get_context().precision)
rss = bbpmpgmp(digstart, digdisp)
irss = int(gmpy2.floor(rss * (16 ** (digdisp))))
print(hex(irss)[2:])
本番
さて表示する16桁が、バイナリファイルで求めている （小数点以下） 8304820223〜8304820238桁 （ちなみに最終桁の最後3ビットはあてにならないので、0でマスクしている） に一致しているかどうか。
これには予想で1時間15分ほどかかるはずなので、台湾ドラマ『翻糖花園』 （日本放映は2012年） の最終回を観ながら実行した。まだ25歳くらい? の簡嫚書が超可愛い♡ 『翻糖花園』は、森七菜主演の日本のTBSドラマ『この恋あたためますか』 （2020年） とどことなくシチュエーション （ケーキ作り美少女がイケメン社長に恋を……とか） やエピソード似ているのだが、台湾・韓国のイケメン男優がどろどろの性格に豹変したり、その結果主役の美少女がとんでもなく悲惨な目に遭う （かわいそうな表情の簡嫚書が、また良い（笑）） のが、低予算のおバカドラマにしてはめちゃくちゃ重くて、良い。いずれにしても美少女がケーキを作るドラマに、はずれはないな……。
緑色: 基本のBBP、青: Bellardの公式、赤・黄: mpfr52ビットと90ビットくらい、紫: 本番 （mpfr約90ビット・8コア並列化）
と、円周率と関係ないようだが （やっぱりない）、こないだ数年ぶりにこのDVDを引っ張り出して以来数週間、24話 （なにしろ1話正味45分あるのが24話分である） 観続けている間と時期を同じくして、わしの円周率は最初のガウス--ルジャンドル→Chudnovsky並列化→プリントアウト・2進化→BBPと進化してきた。面倒くさいアルゴリズムを考えているときにこういうおバカっぽいドラマやアニメは、良い（笑）。
さて予定通りコア8つで800%で、ファンもぶんぶん回っている。予想近くの1時間10分すぎで、プロセスモニタによると4コアの計算が終わったようだ。親はまだ終わっていないので、おそらく子供のループ前半部分が先に終わったのだろう。やはり桁数が多いと、
結局親 （公式4項目のループ後半……つまり一番計算量が多そう） が16分遅れで、80分ほどで全計算が終わった。ドラマの番薯も記憶を取り戻し米恩と仲良くなった （意味不明） ところで、答があっていればわしも幸せである……。
10進100億桁→2進化 の最後の部分のダンプ
15 f8 f4 35 f0 f0 bf 8f d6 48 60 3e 3f a4 43 cc
d5 96 37 ae bc 82 f7 cd c2 67 dd a1 d0 ba 57 a5
c2 cc 5a 85 3d 28 44 5f 63 89 03 b9 8d 5e 02 d3
1d 12 06 40
BBPでの8304820223〜（略)0238桁目
time python3 bbp.py 8304820223
8d5e02d31d120646
python3 bbp.py 8304820223 34287.70s user 42.03s system 718% cpu 1:19:37.43 total
米恩と記憶を取り戻した番薯
ということで、最後のバイトの下の (0クリアした） ビットを除けば、一致していた!! 幸せ!!
なんちて、両方同じ狂い方してたらわからんけどな〜（笑）。
おまけ
まだAIには勝てるな……。
Discussion