円周率を巡る冒険の続き
円周率を100億桁求めたい。
という『円周率とわし』の続きであるが、手許のラップトップ (Macbook Pro) だけで求めたいと半ば意地になった、ある意味バカバカしい記録である。
自分自身ではいろいろためになることもあったが、公開する記事としてはどうか? (笑)
- AWSで実行したバージョン: (Github): https://github.com/tarohs/pi (pichgmmt.py)
- Chudnovskyの公式 (円周率.jp): Ramanujan系公式
- 実装について: Chudnovsky の公式を用いた円周率の計算用メモ
- 実装例: C++ - 円周率計算(Chudnovsky の公式使用)!
Python3で並列化
前記事の『Chudnovskyの呪い』の続き。
こちらの『C++ - 円周率計算(Chudnovsky の公式使用)!』の記事を参考に、ChudnovskyのアルゴリズムをPython3で実装してみた (というかちょちょっと試した時点では、Python3 + gmpy2へのほとんど単純書き換え)のだが、その後まずは並列化してみた。
前記事にもちょっと書いたが、円周率計算は並列化には向いていない。乗算をFFT (NTT)化すれば並列処理はできるだろうが、そこは良くできた子のgmpy2に任せちゃっているので、そもそもgmpy2がマルチプロセス化されていないので、CPUのマルチコアを活用できていない。いちおうMacbook Proだって10コア持っているのだからぶん回したいのだが、計算が宴もたけなわになったときでさえload averageは (1コアの) 60%ほど (だからアクティビティモニタ『CPU使用率』にすれば6%とか) である。涼しい顔でとろとろ計算している。
マルチプロセス化しようと思えば、Python3にはconcurrentというパッケージがあるらしい。ただしプロセス間でのデータの受け渡しにpickleを使っているので、『オブジェクト (変数) がpickle化可能でないとダメ』らしい。なんのこっちゃ。
調べてみると幸い、gmpy2のmpz型その他はpickle化に対応しているようである。さっそく、Binary Splittingのページ中t
、q
、p
を求める部分):
def comppqt(n1, n2, lev):
if n1 + 1 == n2:
p = (2 * n2 - 1) * (6 * n2 - 1) * (6 * n2 - 5)
q = c3_24 * n2 * n2 * n2
t = (a + b * n2) * p
if n2 % 2 == 1:
t = -t
else:
m = (n1 + n2) // 2
p1, q1, t1 = comppqt(n1, m, lev + 1)
p2, q2, t2 = comppqt(m, n2, lev + 1)
q = q1 * q2
p = p1 * p2
t = t1 * q2 + p1 * t2
return p, q, t
を、異なるコアに振ってみた。
from concurrent.futures import ProcessPoolExecutor as Executor
def comppqt(n1, n2, lev):
if n1 + 1 == n2:
p = (2 * n2 - 1) * (6 * n2 - 1) * (6 * n2 - 5)
q = c3_24 * n2 * n2 * n2
t = (a + b * n2) * p
if n2 % 2 == 1:
t = -t
else:
m = (n1 + n2) // 2
if lev <= maxproclev:
with Executor() as exec:
# future1 = exec.submit(comppqt, n1, m, lev + 1)
future2 = exec.submit(comppqt, m, n2, lev + 1)
# p1, q1, t1 = future1.result()
p1, q1, t1 = comppqt(n1, m, lev + 1)
p2, q2, t2 = future2.result()
q = q1 * q2
del q1
p = p1 * p2
del p2
t = t1 * q2 + p1 * t2
else:
p1, q1, t1 = comppqt(n1, m, lev + 1)
p2, q2, t2 = comppqt(m, n2, lev + 1)
q = q1 * q2
p = p1 * p2
t = t1 * q2 + p1 * t2
return p, q, t
lev
は再帰の段数である。はじめ、うっかりこれを考えず走らせたら……あっという間にプロセスを、(2の再帰の段数)乗だけ作り始めて (ちなみに再帰の段数は、桁数maxproclevel == 3
)。
コメントアウトしているところは、親スレッドでは子供二つにforkしなくとも、子供ひとつをforkして自分も計算をはじめればよい、ということである。future2 = exec.submit()
したところでforkして子供が計算をはじめるので、あと自分が次の段数を呼び出し、future2.result()
で合流するようである (順番を間違えると子供が帰ってきてから親が計算を始めてしまう(笑))。
それから当初、mainは地べたに書いていたのであるが、エラーが起きる。エラーメッセージをみると、『きちんとif __name__ == '__main__':
の中に書け』とある。は? と思ったが、そうすると確かにエラーは解消する。どうやら仕組みとして、fork()した子供は、同じPythonスクリプトを読んで目的の関数を実行するみたい。行儀悪くmainを地べたに置いていると、初期化その他が子供にも読まれておかしくなる。随分疎な、コストがかかるforkであるが、頭良い。
実行してみると、ぶん回る。ファンがぶんぶん回るのだ。しかも、1億桁1分!! 最後の
実はこのグラフは、次のピクルス化した後のデータである (本質的にはディスクアクセスの時間が加わっているだけで、オーダは変わっていない)。グラフの紫・緑・青はそれぞれ、再帰による多項式の完了実時間・最後の
ピクルスを作る
いちおーこれで速くなったのであるが、『桁数を増やすとメモリ不足でどかんと落ちる』現象はいかんともしがたい。
macOS12.6でメモリ16GB搭載のMacbook Proであるが、『アクティビティモニタ』でみていると、Python3のメモリ使用量が50GBを超えると、OSだかPythonだかに落とされる (Chromeとか他の文房具ソフトはすべて落とし、数枚のTerminalとFinderだけ開いた状態)。これを超えずに計算できるのは20億桁くらいまで。前記事にあるように、ガウス--ルジャンドルのアルゴリズムでひと晩かけて40億桁の自己レコードを叩き出しているので、速いだけではあまり意味がない。
ここで1泊2日の入院が入る。大腸カメラによるポリープ手術で、まあ日帰り手術に毛が生えたようなものだが、一昨年も同じ手術で入院して病院は暇なことがわかっているので、ラップトップPCなどを持ち込んでみた。ぼ〜っとしていられる時間なのでプログラムをいじったりはしなかったが……ふと考えたのは、変数を酢漬けにして『ピクルス化』できるのであれば、ピクルスにしてディスクに書き込んだり、読んで生野菜に戻したりしながら桁数更新できないか?」ということである。Pythonのpickleというのは、変数 (オブジェクト) にpickle化する (または戻す) 手続き、『直列化』してファイルに読み書きできるバイト列に変換する手続きを定義しておけば、ディスクに読み書きしたり他のプロセスとやり取りしたり (前節のはこれ) できる、ということらしい。
さっそく退院後にやってみた。
汎用化のためにこんなのを作る。
def pkldump(var, fname):
if isinstance(var, str):
return var
fn = pklprefix + fname + '.pkl'
print('dump to {}'.format(fn), file = sys.stderr)
with open(fn, 'wb') as fp:
pickle.dump(var, fp)
return fname
def pklload(var):
if isinstance(var, str):
fn = pklprefix + var + '.pkl'
print('load from {}'.format(fn), file = sys.stderr)
with open(fn, 'rb') as fp:
var = xmpz(pickle.load(fp))
return var
pkldump()
は、mpz変数とラベル (ファイルのパス名の一部) を渡すとpickle化して書き出す。pklload()
はラベルを渡すとピクルスを読み込んで生野菜に戻して返してくれる。だから
a = pkldump(a, 'a')
...
a = pklload(a)
のようにすれば、巨大なmpzオブジェクトa
が不用不急、すなわち今 (...の部分) 使わないとき、短い文字列であるラベルにしといたり、mpzオブジェクトに戻したりできる、という仕組みである。あと、使わなくなった巨大mpzオブジェクトは、まめにdel
する。
具体的なアプローチとしては、再帰的に項を計算する部分で、return
するp
かq
かt
が巨大であればピクルス化して、文字型のラベル (ファイルパスの一部) を返してやる。巨大かどうかの判定は、gmpy2.num_digits()
でmpz型の桁数を取れば一瞬でわかる。実際q
が一番大きくてt
が同じくらいの桁数、p
は小さい。呼んだ側ではpklload()
で元に戻し (文字型かどうか判定しているのはそのためである。もしそうでなければpklload()
はそのままで返す) 計算する。
元に戻すのもいっぺんに行なうのではなく、q
(p
) の計算では上の段のq1
とq2
(p1
とp2
) が必要なのでそこだけ戻し、t
の計算ではt1
とq2
、p1
とt2
だけ必要なのでそこだけ戻す (乗算で必要な2項だけを戻して結果をまたピクルス化、加算のときはまた戻して使う)。
計算の主要部分は下記のようなものを追加する。comppqt()
の第4引数になっているub
という文字列は、項の前半を計算しているのか後半を計算しているのかに応じてu
・b
を渡している。親の段からもらったラベルに付け加えるので、例えばp
の再帰4段目・後半の前半の後半の後半をピクルス化したファイル名は『/PATH/TO/BASENAME-pbubb.pkl
』 のようになる。ピクルスを格納するテンポラリファイル名にこれがないと、並列化したときに同じ用途 (ラベル) のテンポラリファイルにアクセス競合が起きてしまう。
def comppqt(n1, n2, lev, ub):
...
p1, q1, t1 = comppqt(n1, m, lev + 1, ub + 'b')
p2, q2, t2 = comppqt(m, n2, lev + 1, ub + 'u')
...
if not (isinstance(q1, str) or isinstance(q2, str)):
q = q1 * q2
del q1
p = p1 * p2
del p2
t = t1 * q2 + p1 * t2
else:
# either of (p1, q1, t1) or (p2, q2, t2) is pickled
q = pklload(q1)
qq = pklload(q2)
q *= qq
del qq
q = pkldump(q, ub + 'q')
p = pklload(p2)
pp = pklload(p1)
p *= pp
del pp
p = pkldump(p, ub + 'p')
tt = pklload(t1)
q2 = pklload(q2)
tt *= q2
del q2
tt = pkldump(tt, ub + 'tt')
p1 = pklload(p1)
t = pklload(t2)
t *= p1
del p1
tt = pklload(tt)
t += tt
del tt
t = pkldump(t, ub + 't')
if not isinstance(q, str) and (\
# digits_pickle < gmpy2.num_digits(p) or \
digits_pickle < gmpy2.num_digits(q) or \
digits_pickle < gmpy2.num_digits(t)
):
# p, q, r is still mpz type, but too long
p = pkldump(p, ub + 'p')
q = pkldump(q, ub + 'q')
t = pkldump(t, ub + 't')
return p, q, t
このようにして、50億桁まで達成した (総計算時間は1時間47分)。
元のChudnovskyの10億桁に比べたら大きな進歩だが、ガウス--ルジャンドルの40億桁に比べてわずか10億桁のレコード更新だ。調子に乗って100億桁に行こうとしたが……落ちる。60億桁でも落とされる。
アセンブラっていうか、電卓
なんとか桁数が、地球の総人口を超えられないものか(特に意味はなし)。
上記のピクルス化したデータは、使用済みになってもディスクに残している。といっても桁数が莫大になった最後の再帰の数段分だけであるが。上記digits_pickle
は21億桁にしている。10億桁まではストレートに求まったのだから、10億桁どうしの乗算・除算はできるだろう、あまり桁数が少ないとピクルス化・戻す時間がもったいない、ということからである。
100億桁の計算をしてメモリ不足で落ちた後のピクルスをみると、最後の4段分のファイルが残っている。落ちる直前のファイルをみると、再帰の大元の親の2段目 (子) から4段目 (やしゃご) までは、結果を残せているようだ。2段目の
つ・ま・り。2段目のファイルをロードして、上の40億桁を削って精度を保持するために必要な100億桁に戻してやれば、なんとか最後まで計算がいくはずである。
上記のプログラムをコピって改造して、2段目の
そして
前節のプログラム、「ピクルスを戻して、計算して、またピクルス化して……」のコードをみても、まるでレジスタが2つしかないCPUのアセンブラである(爆)。ロードして、計算して、またメモリに戻して。しかもこのアセンブラ、1ステップの実行に1分とか5分とか15分とか掛かる。レジスタは10進100億桁の長大なレジスタだ。
ということで『最後の部分の計算プログラム』は汚いし長いしつまらないので載せない。
このステップはさくっとできたわけではなく、実行しているとエラーで止まる (mpz型と、ピクルス == 文字列型を間違って演算してみたり) ので、エラーが起きたところを修正して、そこまでをコメントアウトして、最後に残っているピクルスをロードする部分を追加して再実行、なんてことを仕事の合間に内職でやっていて1日掛かってしまった……。
だが。
最後の最後の1ステップの、分子を分母で割れば (固定小数点100億桁の)
万事休す?
gmpy2のメモリ効率
後知恵になるが、VMが50GBを超えない範囲でできることはどうやら
- 130億桁
130億桁 → 260億桁\times - 120億桁 ÷ 100億桁 → 20億桁
- 200億桁 ÷ 40億桁 → 160億桁
が限界であった。ストレートな実装 (漸化式に出てくる変数を全部保持した状態) ではメモリを喰うし、必要以外をピクルス化したとしても上記演算といえども、左辺の2項と右辺の1項の合計3つの変数は同時に保持できなければならない。100億桁だって、言うて10進10G桁 == およそ2進3GBくらいのはずであるが、さらに乗算 (FFT) は係数を開いた状態にして log2(n)倍位 2〜数倍くらいのメモリを喰うのだろう (※ 並列処理しない限り、半分にしたのをさらに再帰的に半分にして……の和は2倍程度ですね)。除算に至っては、おそらくニュートン法で逆数を求めているだろうから、200億桁が被除数にくると200億桁どうしの乗算が必要なくらいのメモリが要るのだろうと推測する。
除算は乗算より条件が厳しくなる。(
長い除算
これも入院中の産物(笑)。
いちばん最後のフェーズも、もちろん (一番短い変数) 100億桁に切り詰めました。精度としては、分母が100億桁(ちょっと)、分子が100億桁(ちょっと)あれば、100億桁の固定小数点の
ところが、上のコラムに書いたように、(200億桁 ÷ 100億桁)ができない (固定小数点でストレートに(100億桁 ÷ 100億桁)を計算すると、
ここまでいろいろやってきたので、(上のコラムに書いたように) 被除数120億桁・除数100億桁の除算はできることは知っている。調べたらKnuth先生の頭文字D、いやアルゴリズムDなんてのが出てきたが、引き戻し法 (2進の除算) を効率よくするようなものらしく、そもそも何桁かの除算が使える今回の用途には向いてない。
う〜ん? と考えているうちに、『筆算のアルゴリズム』でいいじゃん、ということに気づいた。といっても、1桁ずつ、除算して、おろして、除算して、おろして、…… (10進ベース) のではなく、20億桁ずつおろしながら (
被除数: 12345678: 8桁 (200億桁)
除数: 789: 3桁 (100億桁)
前提: 5桁÷3桁はできる
(120億桁÷100億桁はできる)
→ 5 - 3 = 2桁ずつおろす
(120億 - 100億 = 20億桁ずつおろす)
15|64|72
------------
789)12345|67|8
11835 (== 789 x 15)
--------- <- 12345 // 789 => 64 ...510
510|67
504 96 (== 789 x 64)
--------- <- 57180 // 789 => 72 ...372
5 71|80 <- 0をpadding
5 68 08 (== 789 x 72)
_______
3 72
1.2345678 / 7.89 = 0.156472471482889733840...
これまた汚な目のプログラムなので省略するが
- まず被除数200億桁を作る
- それを上位120億桁と、下位80億桁に分ける
- これも200億桁を20億桁ずつ割っていき (かろうじて (200億桁 ÷ 20億桁) は落ちずに実行できる) 上位120億桁を作る。それに
を掛けて元の数から引いて、余りを『下位』にする10^{80億}
- これも200億桁を20億桁ずつ割っていき (かろうじて (200億桁 ÷ 20億桁) は落ちずに実行できる) 上位120億桁を作る。それに
- 以下を、『下位』が0桁になるまで繰り返す
- 120億桁を100億桁で割るのは落ちずに実行できるので、商の20億桁と、余り100億桁を得る
- 『下位』をさらに、上位20億桁と (新しい)『下位』に分ける
- 余りに
を掛けて上位20億桁を足す (120億桁)10^{20億}
つまり、10進数なら『下位』の末尾1桁ずつ『おろして』付け加えて次の桁の除算を行なうのを、
これまた、アルゴリズムがこんがらがって発狂しそうになりながら、土曜1日掛かってデバッグしましたね。
そもそも「〜億桁」と呟いていると気が狂いそうになるので、しまいには「200桁を120桁と80桁に分けて除算、余りに20桁付け加えて……」とか呟くという、逆『大阪のおばちゃん』状態。
で、プログラムの実行は、寝る前にセットして3時間ほど掛かりましたけど。
寝る前に頭の20億桁が出たときに、ちょっと覗いてみると『31415...』とおなじみの数列が出ているので、「やった!!」と喜び勇んで寝て、翌朝みてみると……暗記している50桁までが、すでに違う。これ、
まあ40億桁までは (そこそこ) 信用できる値が出ているので途中までの検算はできたわけだけど、そもそも50桁くらいでおかしくなってるって……早めに気づいてよかったね。ってか?
おそらく原因は……
- 140億桁を100億桁に丸めたりとか、途中で精度を落としながら実行したため? →これ、以前にガウス--ルジャンドルアルゴリズムの途中だったかで試してみて、やっちゃダメじゃん? と思ったことがあったのだが。単純な除算・乗算なら、有効数字100桁 × 100桁は200桁精度、有効数字200桁 ÷ 100桁は100桁精度なのだが、漸化式の途中でそれをやると精度はがたがたに落ちる
- 途中の『億桁アセンブラ電卓』での計算式を誤った?→人為的ミスですな
- gmpy2の
num_digits()
が怪しい →これ、試してみるとプラマイ1くらいの誤差がある。例えば、gmpy2.num_digits(gmpy2.mpz(6))
は正しく1
になるが、gmpy2.num_digits(gmpy2.mpz(8))
は2になる。原因は 取ってるからなんだろうけど、その精度? ともかく、『億桁アセンブラ電卓』では桁の分割とかにこれを使っているので、境界で誤差が出たら影響している可能性がある\log
原因は不明だけど、頭にきたので、ここで低スペックなラップトップでの実現はやめて、金に物言わせる。
AWSに逃げる・ごめんよMacbook
そもそもこれ、『Chudnovskyの公式・マルチプロセス版』くらいができた時点で、メモリやCPUコアがふんだんにあるマシンに移せばよかったのです。
いや……移しかけてはみました。自分ではここんとこ、ハイスペックなぶん回しマシンは組んだりしていなかったので、研究室の学生が使っているLinuxデスクトップを「ちょっと貸して」って借りてはみたものの……。
CPUコアが8!! メモリが16GB!! それだけならSSDにスワップ取ればいいが、freeしたらswapが2G!! (意味あんのか?(笑)) swapを128GにしようとしたらSSDの空きが70G!! (512GBあるはずなのに……プレインストールのWindowsをケチ臭く消してない) あとまあsshは使えないしftpも使えないし……。
計算機の能力増加と根性の減少
いちお、その学生ってdeep learningとかやってるんですが……(GPUはそこそこ)。なんか、いまの学生の研究って、「2週間ぶん回す」とか「メモリ (VM) 100GB取る」とかいうのがなくて、手許の『ぱーそなるなコンピュータ』でちょちょっと組んで長くても3分くらいで実行して、という規模に陥りがちなのかな、と思いましたね。まあ学生研究なんてリミットが短いわけですけど、そんな中でもワークステーションが遅かった90年代とか、4桁下のパラメータで5分くらいで念には念を入れて検証して、あと2週間ぶん回す (その間に停電が起きたら泣く)、なんてことがままあったわけです。手許の『ぱーそなるなコンピュータ』の性能がその当時のワークステーションの4桁倍くらいになってるけど、計算の複雑度とか要求される結果とか (まさにdeep learningのSoTAな成果とか) も4桁倍くらい増えてる気がするので、そういう技術力および根性でいいのかな、と。
さてそうなると、AWSに逃げます当然(笑)。逃げる、っつーか、当初の『非力なラップトップPCでどこまで行くか』はどこへやら。もはや、100億の目を見るまで、撤退できませぬ。
まずはささっと無料の1コア16GBメモリとかのインスタンスで、マルチプロセス版のChodnovskyが動くことを確認。
さ〜て、どんだけコア数・メモリがある爆速マシンを借りてやろうか? (笑) 10分とかで計算が終わっちゃったら、USD20/時間とかの高級マシンでも、元が取れちゃう♡ (※ 意味分からん)
という幻想ははかなく。
借りたのは、c6a.16xlargeってインスタンスです。オハイオのマシンでUSD2.448/時間。64coreでメモリが128GB、Python3.10が入っていました。なんかpipがなくてgmpy2のインストールができなかったので、こことか参考にしました。
再帰は6段・64プロセスまで分岐したんだけど……まず1億桁とかで実行してみても……Macbook Proとそれほど変わりない!! つか、遅いくらい!! (泣)
コンテナがどうとかバーチャルマシンがどうとかのオーバヘッドがあるのかもしれないけど、う〜んちょっとこれはがっかりだな……というか、Macbook Proを『非力』と思っていたけど、結構なスペックだったんだね。ごめんよ。
まあ桁数-時間の増加予測だと、4時間ほどで実行が終わる。オハイオのマシンというのはデフォルトの無料のがそれだったので、「いけねいけね」と東京のにしようとしたんだけど、東京のほうが時間あたり料金は高いんだな (50セントくらい)。というか、世界中で北米が最安。まあ電気代とか土地代なんで、当たり前か。
とはいえ、メモリが128GBあるので、なんとか100億桁は実行できました。総実行時間が4時間37分(だいたい予想通り)。ちなみに再帰・64プロセス並列の漸化式計算の部分はさすがに速くて5501秒 (1時間31分)ほど、その後の
後2者はシングルプロセスシングルコアで実行されていたわけで、実際AWSダッシュボードでCPU使用率が1.57% (笑) から上がらないのを悶々とみていたわけですが。ならもっとしょぼいプランを借りれば……とお思いでしょうが、メモリ128GB使えるのは、64コアのプランから、ということになっているみたいなので、しょうがない。あとswapを増やしたら……とも思ったけど、できないかもしれないし、できたとしても (自分のマシンじゃないから) そうする手間も惜しいし、ディスク余計に借りてお金取られたら元も子もないし。
で、頭の100桁くらい(笑)をちらと確認して、手許のマシンにダウンロード〜 (ファイルも10GBも置いておいたらお金取られそう)。……と、遅い(笑)。
しまった。オハイオからだとファイル転送が遅い。結局、時間あたり1ドルとかケチってしまったばっかりに、scpする1時間半で余計にお金が掛かってしまいました (東京からだとどのくらいの速度だったかしらんけど)。
しめて7時間ほど、17ドルちょっと。円安のきょうび2400円。
円周率100億桁の値段としてこれが、高いのか安いのか、よくわからん(笑)。
結論: 円周率100億桁目は
0です (たぶん)。
円周率100億桁・最初の1000桁と最後の1000桁 (たぶん)
elapsed time 1: 5501.65093588829 [sec]
elapsed time 2: 12606.938476324081 [sec]
pi(10000000000):
3.
1415926535 8979323846 2643383279 5028841971 6939937510
5820974944 5923078164 0628620899 8628034825 3421170679
8214808651 3282306647 0938446095 5058223172 5359408128
4811174502 8410270193 8521105559 6446229489 5493038196
4428810975 6659334461 2847564823 3786783165 2712019091
4564856692 3460348610 4543266482 1339360726 0249141273
7245870066 0631558817 4881520920 9628292540 9171536436
7892590360 0113305305 4882046652 1384146951 9415116094
3305727036 5759591953 0921861173 8193261179 3105118548
0744623799 6274956735 1885752724 8912279381 8301194912
9833673362 4406566430 8602139494 6395224737 1907021798
6094370277 0539217176 2931767523 8467481846 7669405132
0005681271 4526356082 7785771342 7577896091 7363717872
1468440901 2249534301 4654958537 1050792279 6892589235
4201995611 2129021960 8640344181 5981362977 4771309960
5187072113 4999999837 2978049951 0597317328 1609631859
5024459455 3469083026 4252230825 3344685035 2619311881
7101000313 7838752886 5875332083 8142061717 7669147303
5982534904 2875546873 1159562863 8823537875 9375195778
1857780532 1712268066 1300192787 6611195909 2164201989: 1000
...
5587569231 1916741447 7702337702 0571791321 1247929546
7195629658 6043092157 1132233849 8808957834 8826655316
5152051372 2330549534 2497091035 5641285479 9020178849: 9999999000
8736018427 8022939085 3632263602 9164960724 9253951680
8239044524 1046598281 9714069883 1849164728 6813215342
8422627907 8451199996 8436917914 5365789613 4202662817
8056580523 9636009121 2256412876 3747009292 0821011851
8780197373 9779761749 1898503651 0263510453 7058348249
9629926346 4929940352 3407787818 4555300243 2561520362
0695010065 2100590569 0146270860 9337228712 5496764316
8804232392 2474677334 6899098720 9641263945 0718828928
9195140928 9334808059 2171247681 8782703796 3876325002
4173880437 4375094501 8520796690 8513958393 8586030364
8332574998 0260967993 9347499516 2191340520 5406737144
6647748614 2058925629 6645192075 8448056490 0980574725
9617692648 9854554620 4750476046 2415063924 3013401091
3189682535 5636077322 8765858346 0774423594 7826580142
4261033692 7895337802 8468073307 2238712577 3387670997
0958650098 8285428578 1732594351 2311990632 4044487901
7204315997 9122170646 1686387699 1811484107 5539249505
7462161487 1660343893 3047927634 9327722831 0815563680
9763261541 1423749758 2083180752 2573977719 9605119144
9403994581 8580686529 2375008092 3106244131 4758821220: 10000000000
Discussion