🧮

円周率 (16進数) の8304820223桁目から16桁を求めたはなし

2022/12/03に公開約19,700字

これも割とぷちネタ。

タイトル通り、円周率の途中だけ求めるBBPの公式ってのを実装して実行してみた。BBPの公式を利用したアルゴリズムは凄い。少ない一定量のメモリで (時間は桁数に比例するが) 円周率を途中から求められる。結構面白くてハマったので自分記録。2022年12月頭の (足掛け) 3日間くらいの記録。

途中から失礼します〜

円周率を10進で100億桁、自分で書いたプログラムで求めてみたが、2進数表現が欲しくなり100億桁の10進数→332億桁の2進数変換をした。だがその末尾数ビットか数バイトは正しい自信が持てなかった。

いろいろやっているときに、円周率を途中だけ求める『BBPの公式』、というのがあるのに気がついた。知らなかったのもそのはず、わたしが学生時代には誰も知らなかった公式である。円周率の長い歴史の中ではごく最近の1995年に、BさんとBさんとPさんが気づいた。

ただしこれは、16進数 (2進数) で任意の桁を求める、という方法である。わしが (10進数) 100億桁ぽっち求めて喜んでいるのは、あわよくば紙かなんかに100億の数字を出力してやれ、とのアピールを目論んでいるからであって、2進数の01や16進数で『0x3.243f6a88...』とか出ても嬉しくもなんともない。なかった。

ちなみに、10進数の\piを求めたくてメモリが足りなくて悩んでいるときにここで「誰か (円周率の途中求めるのを) 16進数バージョンじゃなくて10進数で作って〜!!」と書いたが、いろいろ漁ってみてたらBPPさんたちの発見はすぐに、10進数どころかn進数で求められるように発展 (1997) していたのだった。

まあそれはともかく、10進 (100億桁) →2進 (33219280957ビット) 変換した円周率を検算したくなったので、BBPアルゴリズムを実装してみることにした。

*BBP系公式 (円周率.jp)
*BBPアルゴリズム
*円周率の小数点以下8000兆桁めをGeForceで求める方法~GTC 2013の「π」セッションレポート (4gamers.net)

世の中にはこんなに親切に解説してくださっているページもあるのだが、どーも頭に入って来ない。つか、総じて他のアルゴリズム (arctan公式やモンテカルロ法からガウス--ルジャンドル、Chudnovskyまで) に比べて、BBPには「やってみました」「こうすればいいですよ」系のページが少ない。

世の、円周率何十兆桁とか求めているいわゆる『円周率ガチ勢』におかれてましては、この公式で「ささっと検算して合ってました」的な使い方をされているので (もちろん途中の何京桁目だけ求めるガチ勢もいたり)、そ、そんなにみんな簡単にできるの!? と引け目を感じた。

ところが上の3番めのEd KarrelsさんのGPUを用いた計算方法のプレゼンテーションのスライドを何枚かみていたら、「あ!! そういうことか!!」とやっとわかった (ちなみにきちんと説明してくれている上の1番目・2番めのページでわしが混乱していたのは、小数点以下の桁数の取り方とか、ノーテーションとか高速指数計算 (後述) の引数に何を突っ込むかとか、であった)。

しくみ

ひとの褌で偉そうに解説しようっと(笑)。

BBP公式は、ふつーに\piを求める漸化式である。漸化式は、ふつーは途中がすべて必要である。途中ってのは、項数 (繰り返し回数) だけではなく計算結果を保持する桁数についてもである。出力だけに関して言えば『もうこれ以上値が変わらない桁はさっさと出力しちゃう』アルゴリズムもあるが、保持はしておかないとその後の計算に影響する (知らんけど)。

ところが、BBP公式をあれこれすると、途中の計算結果は劇的に少ない桁数を保持すれば十分、になる。VMが たった50Gしかない (実メモリ16GB) のMacbook Proでも安心である。漸化式の中に1/16^nが出てくるのがミソである。

どーいうことか。上述のEd Karrelsさんの言葉を借りると「小数点以下(N + 1)桁目 (16進数) を求めたければ16^Nを掛けてやればよい ((N + 1)桁目が1桁目にやってくる)」「我々は整数部分には興味がない」ということである。つまり\pi16^Nを掛けても16^N \piが求まるだけであるが、整数部分を無視できるとすればそれは(N + 1)桁目 (16進数) になっている。たとえば10進数で考えると、\pi = 3.141592...の3桁目 (から) が欲しければ、10 ^2を掛けてやると314.1592...となるので、整数部を無視できるなら見事、1(592...)が取り出せる。

この式の右辺は\displaystyle \frac{なんとか}{8n + かんとか}みたいな形をしているので、16^N倍すると結局、\displaystyle \frac{16^{N-k}}{8(N-k) + かんとか}の割り算の小数部分だけをお手軽に計算できればよい (なんとか倍したりとかそれぞれを足したりするのは、小数部分が求まったらそれだけで行なう)。

ここでk0〜N〜\inftyだけど、とりあえず0〜Nくらいの範囲が難しいのでそれを考える。k = N〜\inftyつまり(N - k) = 0〜-\inftyの範囲 (実際-\inftyはない。-20くらいで十分である) は、右辺は1未満の小さな数になるので、そんなに苦労しないで計算できる (このことは次節『基本的な実装』で蒸し返す)。

そうすると、(N-k) = N〜0くらいのすべての数について、指数計算をして、割り算して、整数部分を無視して、足し合わせて、を繰り返さないといけないわけだが。

ここに高速指数計算法 (とわたしは覚えているが、いろいろな名前があるようだ) という、『指数計算をしつつ割り算をしつつ整数部分を無視する』という都合の良いアルゴリズムがある。いや正確に言うと『指数計算を({\rm mod}\ n)の剰余体の上で行なう』もっというと『a ^b {\rm \%}\ c』 (ab乗をcで割った余り) を求めるものである (※ 以下modulo演算 (余り) は、Pythonにならって『%』と記す)。

上の16^{N-k}というのは、でかい (※ 小学生かよ)。どのくらいでかいかというと、最大で16^N。まあ16進83億桁の円周率を求めるんだから、83億桁、とかでしょうねー (※ 投げやり)。でもって(8k + かんとか)はどのくらいかというと、 最大でN = 求めたい16進数の桁数 = 83億は16進数でも10桁くらいだから、10桁くらいである。だから、16^{N-k}(8k + かんとか)で割った余りは10桁以内で (割る数より小さい。だって余りなんだもの)、そして余り10桁を割る数(8k + かんとか) 10桁で割った商は1より小さいと来ている (余りは割る数より小さい。だって余りなんだもの)。んでそれは、元の数 (83億桁) を割る数(8k + かんとか)で割った小数点部分と同じだ。

これが納得できない場合は、10進数ですげー大きな数、たとえば100億 (※ 小学生かよ) を3で割った場合を考えてみるとよい。(100億 ÷ 3) の小数部分、つまり『3333333333.33333...』の小数部分は、 (1 ÷ 3) の小数部分と同じで『.33333...』だ。なぜなら、100億と1は、3を法として合同だから……簡単にいうと、1は100億を3で割った余りだから、である。

問題は『a ^b\ {\rm \%}\ c』 (16^k(8k + かんとか)で割った余り) をどうやって高速で計算するか、であるが、これはあちこちに書いてあるので見てね。簡単にいうと、(a\ {\rm \%}\ c)\times(b\ {\rm \%}\ c) \equiv (a \times b)\ {\rm \%}\ c、すなわち『acで割った余りとbcで割った余りを掛けたものは、abを掛けてからcで割った余りに等しい』ということを利用すれば、aの1乗、2乗、4乗、8乗、……を作っていって(2乗かける2乗は4乗で、4乗かける4乗は8乗だ)、bの2進数表現のビットが1のところだけ掛けて答とすればよい (たとえばb = 5 =0b101だったらa^4 \times a^1)。その過程でやっぱり数が大きくなりそうだったら、こまめに『{\rm \%}\ c』してやる。


BBPアルゴリズムをざっくり言うと

  • \pi16^{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進数の\piの正解とかはあるので、とりあえず出来てきた小数部分を16^{てきとーな桁数}倍して整数化してhex()print()してみればよい。

そのへんで気づいたのだが、今回gmpy2の出番はない (※ 実は後で蘇る)。前にgmpy2を使ったのは整数の加減乗除が速いからであって、まあその部分を固定小数点で書いても、今回は大した桁数でない。高速指数計算法が桁数をびしばし削ってくれるので、ほんとーに出てくる整数といったらデバッグ段階では10進数桁 (10000桁程度なら4桁とか8桁とか) くらいの数しかない。実数部分はさすがに10000回とか積算するのである程度の精度は重要だが、Pythonのfloatで十分である。固定小数点は (頭が混乱するひとつの原因なので) すぐにやめた。

試みに数回のループ計算の生の値をみてみる。小数点シフト量0にして、1回だけでループ計算してみる。いきなり『3.133333』が得られて (0.13は2/16の2、\pi=0x3.243f...の2と近い) 「おぉ」とか呟いてしまう。

\pi =0x3.243f... (16進数) に対して

  • 0回シフト: 小数部 = .1333... = 2.133.../16
  • 1回シフト: 小数部 = .2664... = 4.262.../16
  • 2回シフト: 小数部 = .2463... = 3.941.../16
  • 3回シフト: 小数部 = .9627... = 15.40.../16

っていうか、\pi=3.1415...と3.133333...の差ってなに? と、このあたりで気づいたのだが、どの求め方でも『N〜\infty項はふつーに計算する』とある。ここをまだプログラムに書いていなかった。『ふつー』って何よ? 床屋さんが一番困るのは「ふつーにしてください」だよ? と思ったが、まあふつーにfloatで書け、ってことですよね。N項くらいまでを『ふつーでない』、すなわち16^N倍した小数部分を積算しておいてやったから、あと、16^N倍されてもくそでかくならない (1とかそれ未満) 項は、ふつーに計算して足せ、と。

そして、『ふつーでない』積算が終わった段階では、小数点以下N桁目『だけ』を求めるならその値 (いまは小数点以下1桁目) は正しい (実は『まあまあ正しい』) けど、その次のN+1, N+2, \ldots桁も欲しければ、ここからはふつーにfloatでN+1, N+2, \ldots項目を足し込んでやってくれ、っていうことだ。実は『まあまあ』、と書いたのは、小数点以下1桁目 (N桁目) だけにしても、下の桁 (N + 1桁目とか) からの繰り上がりもあるから、『ふつーでない』計算結果だけじゃ正しいとは限らないよ? 下の桁も計算しないと。ということだ。

欲しいN桁目からどのくらい連続で計算して精度を保証するか、ちょっと考えた。はじめの目的としても、検算したいのは8304820223桁目 (16進) あたりなのだが、8304820222桁目だって正しいかどうかわからない。できれば(略)4820214桁目あたりから連続で10桁くらい欲しい……だがまてよ、(略)20223桁目を正しくする保証をするには(略)0224桁目、……以下からの繰り上がりを考えて、(略)233桁目くらいまでは計算しないといけないんじゃ?

ちなみにある桁 (16進) が繰り上がるかどうかは、単純に考えて確率1/16くらいである。余計にその下のx桁を計算しておけば、N桁目が誤っている確率は1/16^xくらいなので、まあ10桁くらい余分に計算しておけば限りなく誤り確率は減らせるかな、と思った。

ここで『任意の』余計な桁数を『ふつーに』(浮動小数点で) 計算する部分を書こうと思ったので、なんだかんだ丸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))、a ^b \ {\rm \%}\ cを計算してくれることが分かった。おぉ……C言語とかで「可能な限り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回の項が進むときの係数のオーダとして1024^Nが出てきている。つまり1024進数っていうか、どちらにしても2進数とはソリが良いわけではあるが、収束は1項あたり10ビット、16進の4ビットの2.5倍速いわけなのだろう。

すでにBBPは正しく動いているので、コピーして書き換えてみる。頭がこんがらがってきたのは、BBPの公式の4項に対して7項もある……だけではなく、ループ1回毎に \pm が反転している (つまり振動しながら収束するようなイメージなのだろう)。

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を足せばよいだけである。整数部分があまりに大きく (またはマイナスに大きく) なると、限られた浮動小数点の精度の中では下の桁に悪影響があるが、まあまあ(-1. \ldots 0 \ldots +1.)の範囲に収まるのであれば、途中でなにか整数を足してからmodf()するとかの配慮は要らないだろう。

さて話が逸れたが。Bellardのアルゴリズムはそれ以外にも胃に悪い(笑)。1024 = 2^{10}、10ビット単位で求まる、ということは、1024進数的N (ここから見たい) に奇数を指定すると (表示は16進なので) 16進数的には全然似ていない値になる。ここは最終的には『欲しい16進の桁数を指定して→16進と一致する (5の倍数 (の剰余)) 桁数に丸め (プログラムの都合で、具体的には16進3, 8, 13, 18, ...桁を指定すれば1024進数の切れ目と一致する)→ループ回数を決定する』という安易なアプローチに落ち着けた。

その上、N = 0のときにはうまくいかなかったり、64で割るところで何かがまずいらしく下の方のビット (計算誤差で食い違う) だけではなく上の方のビットがおかしかったり (整数化して切り捨てる点がおかしくて切り捨て残ったゴミが混ざる) した。

それでも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) 以上に分けられない。もうひとつは、ループ回数 (k) を複数に分割することである。ループ1回1回についてプロセスを起動→結果集めをしていたら、オーバヘッドのほうが大きくなってしまうだろうから、ある程度ループする領域を大きめに分割してマルチコアに分けたほうがよい。

なんにせよ、『並列化が簡単』すなわち『前の結果に依存しない』漸化式になっている、というのはありがたい。

もうgmpy2 (mpfr精度を指定できる) バージョンまで出来ていたので、それをまずは簡単な上記前者 (公式の項) ごとに分けてみた。方法は、以前使ったのと同じconcurrent.futuresProcessPoolExecutorを使う方法である。

  • リストは本番プログラム (後出) 参照

あれ? なんか結果がおかしい。小さいNでやっているが、下の何ビットかが狂っている。これは明らかに計算誤差である。

といっても、なぜ計算を分割したときだけ計算誤差が発生しているのかわからない。分割点がおかしい (ループ項数に重複とか抜けがある) とか、分割バージョンでmpfr()を宣言し忘れているとかPythonのfloatで計算してからmpfr型に代入しちゃって精度を落としているとか、親子プロセス間のデータの受け渡し (ピクルス化するらしい) とか、いろいろな疑いでプログラムを見直した。

せめて、分割されたそれぞれの項の和の正解がわかれば……と思ったので、まずはループの項数kで8プロセスに分割するプログラムを書いてみた。やはり最終結果に誤差がある。次に、マルチプロセスに分けないプログラムで同じ点で8分割しそれぞれでの積算を表示するプログラムを書いてみた。こちらは最終結果が正しいので『正解』といえるだろう。

比較してみると、ある桁数以下に誤差がある……やっぱりPythonのfloatで計算している?

おやっと思って、マルチプロセスで呼び出されたサブルーチン内でのmpfrの精度を表示させてみると……52ビット (初期値) に戻っている。つまり、mainで計算桁数に応じたビット数 (90ビットくらい) まで精度を上げているのに、子供が起動されるとその中では52ビット計算になっているわけだ。

なるほど……concurrent.futuresの子プロセスの起動の方法を思い出したら、forkした子供がスクリプトを読んで関数を実行する (っぽい) のだった。mainは読み飛ばす。いくらグローバルな属性であっても、子供は精度の設定をしていない。そうか、では子プロセスに入ったところでもう一度精度を設定してやらなければならない。

引数にそれを追加して設定……こんどは、プロセス分割してもエラーがない。公式の4項で分割しても、ループ回数で分割してもエラーがない。そこで、両方で分割して8コア (決め打ち) にすることにした。どーせ10コアMacbook Proで実行しても推定1時間ちょっとである、64コアとかのAWSを借りるまでもない(笑)。

公式の4項それぞれで前半・後半に分割で8つ……前後半の最適な分割点はどこ (欲しいNの何倍) だろう? さいしょは、桁数が多くなると計算に時間がかかるだろう、と思ったので、0.7倍 (N^2オーダの計算量のときここが半分になるはず) にしてみた。と、子供のほうが終わるのに時間が掛かっている。

何通りか試してみた。プロセスは#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コアの計算が終わったようだ。親はまだ終わっていないので、おそらく子供のループ前半部分が先に終わったのだろう。やはり桁数が多いと、kが大きいところで多少時間が掛かるようだ。

結局親 (公式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

ログインするとコメントできます