円周率とわし
円周率40億桁求めた。
たった40億桁でこんなこと言うと円周率ガチ勢には笑われるかもしれないが、今生で自作プログラム (Python3で60行ほど、ラップトップでわずか7時間) で円周率をこんなに求められるとは思わなかったなー。
2022年現在のテクノロジでちょちょっと計算して円周率が計算できた、のを、ハマって多桁・高速化した、2022年11月の1週間くらいの記録。
- 基本バージョン・gmpy2バージョン (Github): https://github.com/tarohs/pi (pigl.py, piglgm.py)
円周率とわし
1980年代に高校のとき、『マイコン同好会』の企画でMZ-80 (Z-80の8bitマイコン) で、文化祭で
んで、記事本編に入る前に、円周率.jp とかみながらいろいろ思ったこと3つ。
- その1: 40億桁くらいっていうと、1995年くらいの世界記録と同じ。就職したころだな〜。27年くらいで、世界の最先端くらいの計算機パワーが電気屋で20万円で売られるのかwww
- その2: では高校生だったころはどうかというと。文化祭ではとにかく桁数ほしかったから円周率はやめたんだけど、アルゴリズム的には (計算時間を度外視すると)
\pi - その3: 円周率.jp その他記事でもわかるけど、ここんとこ10年以上『大型計算機』『スーパーコンピュータ』による円周率の記録はなく、『個人用計算機』(つまりPC→クラウド) の記録ばっかり。スーパーな世界は超並列なベクトル演算の世界に移行しているので、並列化が難しい円周率の計算ではなく、ディープラーニングとか流体とかタンパク質折りたたみとか、そういうのでベンチマークして速いようになってる。
ことのはじまり
codewars とかいうプログラミング修練系サイトでいろいろな問題をみていたら、『
Pythonのトレーニングを選んでいたので早速、ガウス--ルジャンドルのアルゴリズムでプログラムを書いてみた。漸化式の項の精度は不変、なんと1回のループごとに桁数が約2倍2倍 (精度2倍ではない) で収束していく (つまりループ8回で250桁とか・16回で6万5千桁とか・24回で1670万桁とか、みたいになる) ありがたい公式である。
Pythonは無限多倍長整数演算があるので、固定小数点で書くのは簡単である。
上記ページに出てくる
んで早速、主要部分が下記のようなプログラムを書いてみた。
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 (数十〜通りのテストデータによるテスト) が通らない。
遅いのだ。
この問題、チート (結果表を持っておいて
結局、上記のアルゴリズムとPythonの無限多倍長演算では、最適化しても400テスト程度で打ち切りになってしまう。掲示板で質問したところでは、1000テスト通らないとパスできないらしい。えぇ? 1万桁の
というわけで、高速化と多桁化でやっきになってしまった(笑)。
高速化への挑戦
Python3の無限多倍長演算は、それなりに速い。調べてみるとどーやら、Karatsuba法を使っているようである。2桁
以前、ESP32だったかSTM32だったかのマイコンで100桁程度のBCD乗算を高速化したくて、C言語でKaratsuba法の乗算ルーチンを書いたことはある。Karatsuba法は、数百〜数千桁くらいならたしかに速いと思える。速いけど、
計算機での試行時間の絶対値ってのは、計算機にとっては何ともないかもしれないが、ひとにとってはバカにならない。いや何を言ってるかというと、例えば1msecかかるものが1秒になっても、この差1000倍だが、ひとは何とも思わない。もちろんGHzクロックで動いている計算機にとって1msecは100万回で1秒は10億回だが、計算機は文句は言わない(笑)。
次の1秒の1000倍は約20分になるが、「ちょっと風呂に入ってくる」くらいの時間ができる。さらに1000倍になると十数日なので、停電対策とか放熱とか止めずに実行したりエラーを起こさない工夫が必要である。そのまた1000倍は1万数千日、もう30年オーバーなので若者は中高年に、中年以降なら人生が終わってしまうかもしれない(笑)。計算機だって、動いてれば文句こそ言わないが、ふつーに言って30年無故障連続稼働してる機械ってのはそうそうないだろう。よく考えたら当たり前のことなのだが、よく考えずにループ回数とか設定すると、人生で終わらないプログラムが簡単に書けてしまう。
だから計算量オーダはもちろん重要で、ベースになる計算機の速さとか掛け算の桁数回数なんかもやっぱり重要なわけで。まあ前述の『MZ-80による……』の時代を考えると、自力 (あまり頭使ってない) で100万桁くらい求められる、ってのは感動なんですけどもね。
なお、どうして1980年代に8bitマイコンで『
閑話休題。最近、競技プログラミングで『〜 (←組み合わせとかの問題) の答を998244353で割った余りの数を答えよ』っていう問題をみたことがある。なんじゃこりゃ？ と思って調べたら、こーいうこと や こーいうこと らしい (驚)。
つまり要点は
- 乗算は、桁をずらして掛けて足す→畳み込み演算なので、FFTを使えば超高速乗算ができる
- 浮動小数点
\sin \cos {\rm mod}\ p
ということなのだが。驚いたのは、わしはFFTもよく使うし有限体上での演算も使う。その世界でのFFT (NTT)もきいたことがあるが、これを超高速乗算に使うとは知らなんだ。不勉強にもほどがある。さらに驚いたのは、競技プログラミングの参加者って素数体上の原始根とかに知見があってて、自前の (
まあそれはともかく。そんなのを実装したのが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
円周率計算はスーパーなコンピューティングが目指す道、並列化・ベクトル化に向いていないかというと、ちょっとそうでもない気もする……が、よくわからない。
漸化式で計算する円周率アルゴリズムは、前の結果が出ないと次の計算ができない。だから順次計算することが必要で、多数の計算機が部分計算を分担する並列計算には向いていない。
用いる演算の中では、掛け算を速くすればよいわけだが(というか加減算を除くと乗算しかない。√はニュートン法で乗除算で求めるし、除算は逆数をニュートン法で求めて掛けるだけだ)、そこをFFT (NTT) 化すれば並列化の余地は残る。
ところで、わしはごく最近話題になっていたのをみて知ったのだが、
要は16が出てくる漸化式にして上の桁を引き算して (違うか？) ※上の桁を整数・下の桁を小数にして整数部分を除去して 、桁をバラしました、ってことみたいなのだが、残念ながら10進数に変換するには
ただ、 円周率.jp によると、現在の
一瞬、それぞれに ※ メモリはたいして必要なく、実際にGPU実装の例もあり。この記事中の『円周率の小数点以下8000兆桁めをGeForceで求める方法』のリンクを参照。この記事でも並列化してみているが、それぞれのプロセスの使用メモリも少なかった。
なお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億桁分の
これに対して、上記実装のChudnovskyのアルゴリズムでは、漸化式の項 (上記ページの
項数が少ないうちは計算桁数が短いのは速い理由でもあるのだろうが、次第にメモリ使用量が増えていって、あるところでどかんと落ちる。
ちなみに、ガウス--ルジャンドル法でも『結果の精度が低いうちはbaseの桁数を伸ばす) を取ってみた。そしたら結果が
多桁化では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する、とかしないと同じく落とされてしまうと思うので、ちゃっちゃっと無改造で実行するのから比べたらかなり敷居が高い。
まあともかく、地球の人口オーダの桁数の
また人生のどこかで、これ以上の
Discussion