🧮

円周率とわし

2022/11/13に公開約11,600字

円周率40億桁求めた。

たった40億桁でこんなこと言うと円周率ガチ勢には笑われるかもしれないが、今生で自作プログラム (Python3で60行ほど、ラップトップでわずか7時間) で円周率をこんなに求められるとは思わなかったなー。

2022年現在のテクノロジでちょちょっと計算して円周率が計算できた、のを、ハマって多桁・高速化した、2022年11月の1週間くらいの記録。

円周率とわし

1980年代に高校のとき、『マイコン同好会』の企画でMZ-80 (Z-80の8bitマイコン) で、文化祭でe = 2.718281\ldots (自然対数の底) を何十時間かかけて (確か) 2万桁求めた。文化祭をみにきてくれた、一般客のおじさんがそのプリントアウトをみて、褒めてくれた。……それも感無量だったのだが、上記の自己レコードも世界記録にははるか及ばないが、そのときの計算時間と桁数を思うとやはり感無量である。

んで、記事本編に入る前に、円周率.jp とかみながらいろいろ思ったこと3つ。

  • その1: 40億桁くらいっていうと、1995年くらいの世界記録と同じ。就職したころだな〜。27年くらいで、世界の最先端くらいの計算機パワーが電気屋で20万円で売られるのかwww
  • その2: では高校生だったころはどうかというと。文化祭ではとにかく桁数ほしかったから円周率はやめたんだけど、アルゴリズム的には (計算時間を度外視すると) \piでも1万桁くらい求められたはず。その当時の世界記録が1千万桁くらい。つまり (最先端ガチ勢:素人わし) の比率が 1000:1 くらい。いまはどうかというと、世界記録が100兆桁、わしの (多分ここに書くことは、高校生程度の知力・財力ならやっぱりなんとかできる) が40億桁で25000:1 くらいに開いてしまっている。
  • その3: 円周率.jp その他記事でもわかるけど、ここんとこ10年以上『大型計算機』『スーパーコンピュータ』による円周率の記録はなく、『個人用計算機』(つまりPC→クラウド) の記録ばっかり。スーパーな世界は超並列なベクトル演算の世界に移行しているので、並列化が難しい円周率の計算ではなく、ディープラーニングとか流体とかタンパク質折りたたみとか、そういうのでベンチマークして速いようになってる。

ことのはじまり

codewars とかいうプログラミング修練系サイトでいろいろな問題をみていたら、『\pik桁目 (k \leq 10000) を求めよ』という問題があった。ちなみにここは、プログラミングを空手かなんかのトレーニングに見立てていて、トレーニング問題には8級〜1級とか級 (kyu) が付けられているのだが、この問題が4級だったか。他の4級問題と比べてみても、かなり無茶振りな4級である (問題の難易度はそこまでではないが、後述・実行制限時間が無茶なので)。

Pythonのトレーニングを選んでいたので早速、ガウス--ルジャンドルのアルゴリズムでプログラムを書いてみた。漸化式の項の精度は不変、なんと1回のループごとに桁数が約2倍2倍 (精度2倍ではない) で収束していく (つまりループ8回で250桁とか・16回で6万5千桁とか・24回で1670万桁とか、みたいになる) ありがたい公式である。

Pythonは無限多倍長整数演算があるので、固定小数点で書くのは簡単である。
上記ページに出てくるa_nとかb_nとかt_nは、0.85とか0.22あたりをうろうろしてるので、それを求めたい精度で整数になる倍を掛けて (例えば5桁なら1万倍して8500と2200とか、10000桁なら10^{10001}倍する) 計算してやればよい。p1, 2, 4, \ldotsのようなのでそのままの整数値でよい。

んで早速、主要部分が下記のようなプログラムを書いてみた。

    base = 10 ** precdig
    gerr = 10 ** errdig
# init: a <- 1, b <- sqrt(2), t <- 1/4, p <- 1 
    a = base
    b = lsqrt(base * base // 2)
    t = base // 4
    p = mpz(1)

    while True:
        an = (a + b) // 2
        b = lsqrt(a * b)
        tt = a - an
        terr = tt * tt // base * p
        t = t - terr
        p = p * 2
        a = an
        if terr < gerr:
            break
    pi = a + b
    pi = pi * pi // t // 4
  • 関数lsqrt()は、ニュートン法で整数の√を求めるルーチンである。ただし引数は必要桁数の2倍の桁数 (base * base のオーダ)
  • precdigは求めたい桁数、baseは整数化の倍率。errdigは末尾の誤差を防ぐため余計に計算する桁数 (8で十分なようだ)

んでもって、codewarsの判定に掛けてみた。もちろんTEST (3通りくらいのテストデータによるテスト) は問題なく通る。

ところがATTEMPT (数十〜通りのテストデータによるテスト) が通らない。

遅いのだ。

この問題、チート (結果表を持っておいてkによって結果を出力するだけ、とか、webから正解を引っ張ってくる、とか) を防止するためのテストも含まれている。もちろんチートはしていないが、ランダムデータによるテストが延々何百回も繰り返されて、それが合計12秒で終わらないとパスしない。

結局、上記のアルゴリズムとPythonの無限多倍長演算では、最適化しても400テスト程度で打ち切りになってしまう。掲示板で質問したところでは、1000テスト通らないとパスできないらしい。えぇ? 1万桁の\piが12msecとかで求まるもんなの?

というわけで、高速化と多桁化でやっきになってしまった(笑)。

高速化への挑戦

Python3の無限多倍長演算は、それなりに速い。調べてみるとどーやら、Karatsuba法を使っているようである。2桁 \times 2桁の筆算の掛け算で4回の乗算が必要なのを3回に落とす (それを再帰的に繰り返す) のがKaratsuba法である。

以前、ESP32だったかSTM32だったかのマイコンで100桁程度のBCD乗算を高速化したくて、C言語でKaratsuba法の乗算ルーチンを書いたことはある。Karatsuba法は、数百〜数千桁くらいならたしかに速いと思える。速いけど、O(n^2)O(n^{1.59})に減る程度である。上記のプログラムで10万桁、100万桁くらいまで増やすと、Python整数演算ではどーにも遅くなってくる。いや遅くなり方は一定なのだが『人間が耐えられない』オーダになってくる。

計算機での試行時間の絶対値ってのは、計算機にとっては何ともないかもしれないが、ひとにとってはバカにならない。いや何を言ってるかというと、例えば1msecかかるものが1秒になっても、この差1000倍だが、ひとは何とも思わない。もちろんGHzクロックで動いている計算機にとって1msecは100万回で1秒は10億回だが、計算機は文句は言わない(笑)。

次の1秒の1000倍は約20分になるが、「ちょっと風呂に入ってくる」くらいの時間ができる。さらに1000倍になると十数日なので、停電対策とか放熱とか止めずに実行したりエラーを起こさない工夫が必要である。そのまた1000倍は1万数千日、もう30年オーバーなので若者は中高年に、中年以降なら人生が終わってしまうかもしれない(笑)。計算機だって、動いてれば文句こそ言わないが、ふつーに言って30年無故障連続稼働してる機械ってのはそうそうないだろう。よく考えたら当たり前のことなのだが、よく考えずにループ回数とか設定すると、人生で終わらないプログラムが簡単に書けてしまう。

だから計算量オーダはもちろん重要で、ベースになる計算機の速さとか掛け算の桁数回数なんかもやっぱり重要なわけで。まあ前述の『MZ-80による……』の時代を考えると、自力 (あまり頭使ってない) で100万桁くらい求められる、ってのは感動なんですけどもね。

なお、どうして1980年代に8bitマイコンで『eを』『2万桁』求めたのかというと、「すげー (※ 当社比・当時比) 桁数の無理数の計算結果を展示したかったから」だ(笑)。ちょっと覚えているのは、eのほうが\piより (いま考えるとかなり収束が悪い) 漸化式に出てくるレジスタ (変数) 数が少ないから、だったと記憶している。たしか32KBメモリ空間で1バイトあたり10進 (BCD) 2桁で演算していたので、3つ (だったか) のレジスタで求められるeだと2万桁いくけど、もっとたくさんのレジスタが必要な\piだと1万桁しか計算できないから、みたいな理由だったと思う。

閑話休題。最近、競技プログラミングで『〜 (←組み合わせとかの問題) の答を998244353で割った余りの数を答えよ』っていう問題をみたことがある。なんじゃこりゃ? と思って調べたら、こーいうことこーいうこと らしい (驚)。

つまり要点は

  • 乗算は、桁をずらして掛けて足す→畳み込み演算なので、FFTを使えば超高速乗算ができる
  • 浮動小数点\sin\cosを使わない、整数 ({\rm mod}\ p) の世界でのFFTはNTT (FMT) という

ということなのだが。驚いたのは、わしはFFTもよく使うし有限体上での演算も使う。その世界でのFFT (NTT)もきいたことがあるが、これを超高速乗算に使うとは知らなんだ。不勉強にもほどがある。さらに驚いたのは、競技プログラミングの参加者って素数体上の原始根とかに知見があってて、自前の ({\rm mod}\ 998244353) 超高速乗算ルーチンとか持ってんのか? ということである。

まあそれはともかく。そんなのを実装したのがGMPらしい。速い・無限多倍長、ってことで、存在は知ってたが、それをPythonで使えるようにするgmpy2って……そもそもPythonの乗算は多倍長だし速いし、誰得? と思っていた。

ところがだ。実際、上のプログラムをgmpy2バージョンに書き換えてみたら、1万桁でPython整数演算の10000倍くらい速いのだ。具体的には、Python3.10の多倍長整数演算を使って数十秒だったのが、0.038秒とかで終わってしまう。ついでに自前の (Pythonループによる) √もクビにしてgmpy2.isqrt()を使った。

import math
import gmpy2
from gmpy2 import mpz

def pi(k):
    global base

    tprecdig = k
    errdig = 8
    precdig = tprecdig + errdig
    base = mpz(10) ** mpz(precdig)
    gerr = mpz(10) ** mpz(errdig)

    a = base
    b = base // 2
    b = gmpy2.isqrt(b * base)
    t = base // 4
    p = mpz(1)
    while True:
        an = (a + b) // 2
        b = gmpy2.isqrt(a * b)
        tt = a - an
        terr = tt * tt // base * p
        t = t - terr
        p = p * 2
        a = an
        print(gmpy2.num_digits(terr), end = ' \r', file = sys.stderr)
        if terr < gerr:
            break

    pi = a + b
    pi = pi * pi // t // 4
    return pi
  • print(gmpy2.num_digits()...)は、収束が目で見てわかるために入れている (イテレーションの前回からのtの差分の桁数)

見ての通り、ほとんどPythonのint型を用いた前出のとコードは同じである。gmpy2.mpz型であると宣言するために、最初にオブジェクトを作るところで余計なコードを書いてやれば、あとはmpz型とmpz型 (int型でもよい) の演算は、勝手にgmpy2に渡されて、超高速演算してくれる (結果オブジェクトもmpz型)。ビバ!! Pythonの型!!

で、肝心のcodewarsを思い出して、上記を書き直してATTEMPTしてみた。

速い。

1万桁の1000回テストも、2秒とかで終わってしまう……のだが、gmpy2をimportするとチート判定になってしまうらしい(笑)。どーせいっちゅーんじゃ。いやまあ、もはやcodewarsはどーでもよい。

多桁化への挑戦

ある程度速くなってくると、今度は多桁を求めてみたくなる。

手元のラップトップ (ARMアーキテクチャのMacbook Pro) では

  • 1千万桁で20秒ちょっと。桁数が10倍になると、およそかかる時間は14倍程度になるようである。

(横軸は求める桁数の桁数w・縦軸秒)

ちなみにこれ、同じARMのMacbook Airでもほとんど遅くならなかった。プロセスモニタでみてても、CPU使用率は10%ほど(笑)。10コアなので、1コアだけしか使われていない、つまりマルチスレッド化されているとは考えにくい。なるほどスーパーコンピュータの世界を、PCも追いかけている。

  • 1億桁で5分ほど。日本人なので『億』って単位になると、ちょっと萌える。日本の人口約1億だし。
  • 10億桁で1時間ちょっと。まあ言うて10億桁ってのは、1ギガ桁である。結果はふつーに1文字1byteのファイルにprint()しているわけだから、1GBの円周率ファイルが出来上がる (あと計算時間と、読みやすくするため小数点以下はpretty printした)。

これなら100億〜1000億桁はちょっとイヤなくらい待たないといけなくても、50億桁くらいは一晩寝ている間に回せば出るに違いない。一晩回してみた。が、なぜか計算は終わっているのに空のファイルが出来上がる (続く)。

\piの10京桁目

円周率計算はスーパーなコンピューティングが目指す道、並列化・ベクトル化に向いていないかというと、ちょっとそうでもない気もする……が、よくわからない。

漸化式で計算する円周率アルゴリズムは、前の結果が出ないと次の計算ができない。だから順次計算することが必要で、多数の計算機が部分計算を分担する並列計算には向いていない。

用いる演算の中では、掛け算を速くすればよいわけだが(というか加減算を除くと乗算しかない。√はニュートン法で乗除算で求めるし、除算は逆数をニュートン法で求めて掛けるだけだ)、そこをFFT (NTT) 化すれば並列化の余地は残る。

ところで、わしはごく最近話題になっていたのをみて知ったのだが、\piの任意n桁目 (ただし16進数) をO(n)で求めるBPPアルゴリズム ってのが最近 (とはいっても30年くらい前……長い\piの歴史の中では最近だ)、発見されたらしい。

要は16が出てくる漸化式にして上の桁を引き算して桁をバラしました、ってことみたいなのだが (違うか?)、残念ながら10進数に変換するにはn桁目より上の桁も必要だ (5で割った余りが伝播する) から「16進数でn桁目が求まっ」ても「十進数でn桁目を求める」のは不可能である。また結局、nのオーダってことは下の桁になるほど時間がかかる、のは避けられない。誰か10が出てくる漸化式作って。いや原理的にできないの? (笑)

ただ、 円周率.jp によると、現在の\pi (全体) の世界記録が100兆桁、それに対してBPPの部分桁記録が1000倍桁も下の10京桁目になっている。なかなか夢がある。

一瞬、n台の計算機にばら撒けばO(n) (っていわねーか、並列化後の実行時間) で\piの全桁が求まる? と思ってしまったが……結局現在のnは百兆とかになってるので、百兆台ものコンピュータを用意するのがまず無理である (1台で複数回実行するにしても、最下桁を割り当てられた計算機とかは掛り切りになっちゃうわけですぜ)。それぞれにO(n)のメモリも必要だろうから、GPUみたいなのにばら撒くわけにもいかないだろう。

なお10進数ではお馴染み『3.141592……』は、16進数では『3.243f6a88……』になるらしい。口に出してみると、なかなかとんがっている(笑)。

Chudnovsky の呪い

円周率.jp をみていると、今年 (2022年3月、岩尾さん) 打ち立てられた100兆桁の記録を始め、最近はChudnovskyアルゴリズム ってのが主役らしい。これはRamanujan系の公式の一種らしい。Ramanujanは、確か円周率の公式だったと思うが、夢に公式をみて起きて確かめてみると正解だった、みたいな伝説があるひとである (これはたぶん都市伝説で、何らかの数学的根拠があったんだと思うけど)。

Rmanujan系、というだけあって、Chudnovskyも謎の13591409とか545140134とか640320とかいう定数が出てくる。素因数分解しても謎な数字だ。

ともかくPython3 + gmpy2で、級数を再帰的に2分割して計算するアルゴリズムで、こんな実装例のページ

なども参考に、とりあえずgmpy2整数演算で書いてみる。

速い。確かに速い。

1億桁で、ガウス--ルジャンドルの5分が78秒ほどになっている。桁数10倍になると時間が20倍に増えるのが、ややガウス--ルジャンドルより悪いオーダみたいなのだが (計算そのものは変わると思えないが、メモリアロケーションとか再帰呼び出しのオーバヘッドだろう)、そもそも4倍くらい速いので1兆桁とかにならない限り負けることはないだろう。

調子に乗って、10億桁に挑戦してみた。

落ちる(笑)。

どうやら、メモリ (VM) をバカ喰いしているので、途中でPythonかOSに落とされているみたいなのだ。

多倍長整数演算だけで作ったガウス--ルジャンドルでは、はじめから終わりまで同じ桁数で計算している。つまり1億桁なら1億桁分のabt分のメモリは最初から終わりまで必要である (オブジェクトとしては生成・消滅を繰り返しているが、必要量のトータルは変わらない)。掛け算をしても2倍、2億桁である。起動してちょっとで落ちなければ、あと最後まで落ちない。

これに対して、上記実装のChudnovskyのアルゴリズムでは、漸化式の項 (上記ページのpqt) の桁数が伸びていく。その上、再帰的に計算しているために、その何倍かの値 (たぶん最大で2倍?) がスタックされていくのだろう。計算順序でなんとかなると思うが、3乗とか出てくると、そのままでは桁数の3倍のメモリが必要である。

項数が少ないうちは計算桁数が短いのは速い理由でもあるのだろうが、次第にメモリ使用量が増えていって、あるところでどかんと落ちる。

ちなみに、ガウス--ルジャンドル法でも『結果の精度が低いうちはabtの桁数を落とす』アプローチ (結果の精度以下の桁数は0にして、どんどんbaseの桁数を伸ばす) を取ってみた。そしたら結果が\piに収束しないのだ(笑)。まあ当たり前か。

多桁化ではChudnovskyは諦めて、ガウス--ルジャンドルに戻ることにした。

Python3.11の呪いとか

そうこうしているうちに、Python3.10→3.11へのバージョンアップが発表された。なんでも実行速度が25%速いとか。嬉しいんだが、まあこのプログラムで時間を喰っているのはgmpy2の (アセンブラで書かれた) 部分で、ループ部分とかほとんど回数がないのだけど (実際、3.11にしても高速化の恩恵は誤差の範囲だった)。

ところが、Python3.11ではint型をstr型に変換する部分で『改悪』がなされた。いやセキュリティのためとかおっしゃるのだが (確かに1億桁のintをstrに変換すると、それだけで何十秒か止まる):

上記ページで紹介がある通り、sys.set_int_max_str_digits(precdig + errdig)とかしておかないと、最後に文字列化 (pretty printのため) する部分で落ちる。いやそれ以前に、print()するだけで落ちる。

さて上記を修正して40億桁に挑戦して (7時間程度で終わる計算になる)、一晩寝てみた。

起きてみたら、やっぱり落ちている。

曰く。sys.set_int_max_str_digits()の引数が不正だとか……? あ、そうか。
この引数はいわゆるふつーの整数型 (符号ありの31ビットなのだろう) なので、10進で20億 (31ビット) までは表わせるが、30億 (32ビット) を超えると、そもそもsys.set_int_max_str_digits()に与えられない……というか (引数がこうなんで内部でもそうだろうから) str型に変換するループのカウンタとかも溢れてしまうもしれない。

この問題は、実はgmpy2.mpz型の整数演算の結果をprint()するときなどはスルーできる。gmpy2は自前で2進→10進文字列変換をもってるのだろう。pretty printしようと思わず、素で (1行につながった) 結果を出すだけなら40億桁まではできたし、なんなら1億桁くらいずつぶった切って (mpz型整数のレベルで商と剰余に分けて) str化すれば、問題なく文字列化できる。面倒なので素で出力したのを、後でCのプログラムで整形したけど(笑)。

今後の展望

さていまのところ自己レコードは40億桁である。なんで40億かというと、結果10進数テキスト4GBでFAT32の1ファイルに収まる(だっけか)……、とかはあまり関係ないけど(笑)。これが労せず (Python3でちゃっちゃっと数十行で書いたレベルで) どのくらい増やせるかは、やはりPCの搭載メモリ量によるようである。

前述、速度はMacbook AirとProで変わらない、と書いたが、実は8GBメモリのAirの方では10億 (1G) 桁の実行に失敗している (ProでChudnovskyしたときと同じく、OSだかPythonだかに落とされている)。gmpy2.mpz型のメモリの使用効率やFFT (NTT)演算のメモリの使用方法はよくわからんのだが、少なくとも16GBメモリのMacbook Proでは、『メモリプレッシャー』 (おそらくmalloc量の微分値) もほとんど赤レベルになることなく40億桁の演算ができているので、物理メモリの量に左右されているのは間違いない。

これ以上はgmpy2を改造して、オブジェクトを分割してswapする、つまりFFT (NTT)乗算の過程で乗数の一部をストレージにswap outする、とかしないと同じく落とされてしまうと思うので、ちゃっちゃっと無改造で実行するのから比べたらかなり敷居が高い。

まあともかく、地球の人口オーダの桁数の\piが今生で自力でちゃっちゃっと求められた、というところで、いちおう満足はしてるのだが。『文化祭をみにきた一般人のおじさん』も驚いてくれるだろう。

また人生のどこかで、これ以上の\piにチャレンジすることがあるだろうか。

Discussion

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