🧮

e (自然対数の底) を100億桁求める〜ひとはなぜ円周率に惹かれるのか?

2022/12/08に公開

Python

またまた円周率絡みのぷちネタ。

円周率を自作プログラムで100億桁求めた。アルゴリズムも3つくらい試したし、バイナリダンプもした。全部は無理だがいちおうかなりの無理筋桁を一度にプリントアウトもしてみた。の続き。

ほかにもなんか、計算して面白い数字はないものか? ということで、40年ぶりくらいに、自然対数の底 (その他) を計算してみた。

まあネタとしては、binary splitting法を理解すること (※ ちゃんとしてないが) と、AWSを借りることについての考察の続きも。

ルート2に人権はないのか?

わしの計算機には、円周率 $\pi$ が100億桁転がっている。

Twitterで流行ったコピペ的に言うと

円周率100億桁が家にあると、ちょっと嫌なことがあっても「まあ家に帰れば円周率100億桁あるしな」ってなるし仕事でむかつく人に会っても「そんな口きいていいのか？私は自宅で円周率100億桁とよろしくやってる身だぞ」ってなれる。戦闘力を求められる現代社会において円周率100億桁と同棲することは有効

ということになるが、真面目な話『唯一無二で揺るぎがない数列』を100億桁ストレージに飼ってると、わしの計算機はメモリは貧弱でもグラボはなくても、世界的に有名なアレが大量にある、という勝った感がある。いやそれは言いすぎだが、一見乱数にみえるこの数列も世界で唯一無二の絶対的真実、求め方を変えても同じもんが出てくる、と思うと、なんか愛着が湧いてきたんだよね最近。さらにこのページで言いたいのは、計算して面白い数という意味でも、やっぱり $\pi$ は断然魅力的だ。

ランキング的にいうと、つまらないほうのトップがまず整数。「整数1を小数点以下100億桁計算しました」とか言ったら、頭湧いてんじゃないかと疑われる。

次に、有理数。まあこれも例えば「 $1/7$ を100億桁計算しました」というひとがいないのは、『0.142857』に続いて、ず〜っと『142857 142857』と繰り返すからだ。ということは、小学生でも気が効いたやつなら知っている。Pythonではこうすれば固定小数点で100億桁くらい計算できる (※ 計算してないが)。

import gmpy2

s = gmpy2.mpz('142857' * (100000000002 // 6))

まあでもこんなことしてると、固定小数点演算のデバッグのためにまじめに $1/7$ を100億桁計算してみたりするんですけどね(笑)。

じゃあ無理数は面白いかというと。あまり「 $\sqrt{2}$ を100億桁計算しましたのでギネス記録」とかいう新記録はきいたことがないし、ぐぐっても『 $\sqrt{2}$ を1万桁掲載!!』とかいうページはみつからない。ちょっと数学の心得があるひとなら、「 $\sqrt{2}$ は $x^2 - 2 = 0$ の解だから超越数じゃない」とかおっしゃるだろうが、まて。超越してると面白いのか?

たしかに、自然対数の底 $e$ (ネイピア数) ってのは超越しているので面白いのか『自然対数の対数の底1万桁』なるページも、ちゃんと見つかるし、100万桁を書籍にしたものも存在する。つかわし自身、ここで書いたように1980年代の8ビットマイコンで2万桁計算しプリントアウトしたものを高校の文化祭に展示し、ひとびとの瞠目を一身に集めた (※ 集めてない。おじさんがひとり、感心してくれただけだ)。

だけど同じ超越数でも、円周率には一歩、引けを取る。まあ $e$ は漸化式が簡単で、多桁求めようとすると円周率ほどレジスタ (メモリ) が必要ないから……というのは大人の後知恵である。思うに『引数があるから』良くないんじゃないだろうか? $\sqrt{2}$ や、 $e$ だって $\exp(1)$ という関数の値 (まあ $e^1$ ってことだけど) である。

つまり、「 $\sqrt{2}$ を100億桁計算(略)」とか「 $\exp(1)$ を100億(略)」とか威張っても、「じゃあ俺は $\sqrt{3}$ でギネス記録www」「 $\exp(2)$ はまだ〜?」と突っ込まれるのがオチである (※ いや突っ込みませんて)。 $\sin(x)$ とか $\cos(x)$ もたいていの (※ 当社比) $x$ では超越数になるけど、「 $\sin(1)$ を100億桁(略)」ってひとがなかなかいないのは、「 $\sin(2)$ はまだ〜?」と言われるからだろう。そこんとこ、 $\pi$ には引数がない(笑)。 $\pi(2)$ とかいうのが考えられない唯一無二性があるからである。ほんとかよ?

あと、 $\pi$ はけっこーいろいろな ~~変態な~~ 数の根幹である。

このページのあとのほうの話題で $e$ を求め終わったあと、「もっとbinary splitting methodが使えそうなおもろい関数は居らぬのか!? 誰か!? 誰かおらぬか!?」「殿、落ち着きなされ」「リーマンゼータ関数とつ $\zeta(いくつか) = \displaystyle \sum _{n = 1} ^{\infty} \frac{1}{n ^{いくつか}}$ はどうじゃ?」なんて叫びながらぐぐったら、3秒で『 $\zeta(2) = \pi^2/6$ 』とかいう式がみつかっちまった。

つまり、100億桁の $\pi$ が家にあると、 (いままでのはなしの流れで、固定小数点表示で $10^{100億}$ 倍されて整数になっているとすると)

import gmpy2

pi = gmpy2.mpz(314159267...(100億1桁))
base = gmpy2.mpz(10) ** gmpy2.mpz(10000000000)
zeta2 = pi * pi // base // 6

gmpy2でこんなコードを書けば、一瞬で (※ うそ。メモリが十分にあったら10分くらいかな?) $\zeta(2)$ は固定小数点で100億桁求まってしまうのだ。死にそうになりながらbinary splittingする漸化式とか考えるより、速い。

あちなみに。 $\sqrt{2}$ を100億桁ですが、これも固定小数点でよければ

import gmpy2

sqrt2 = gmpy2.isqrt(2 * gmpy2.mpz(10) ** gmpy2.mpz(20000000000))

とかで、一瞬で求まってしまうので念の為(笑)。

√2 100億桁 (頭と尻の100桁)

Macbook Pro (ARM)で計算15分、2進→10進変換 & print 55分くらい(笑)でした。

1.
4142135623 7309504880 1688724209 6980785696 7187537694
8073176679 7379907324 7846210703 8850387534 3276415727
... 
1633199632 9587931035 4586453811 9433834694 3776113126
5250836258 6378469245 0310196801 1178369494 4997699614

このVM量では (200億桁 ÷ 100億桁) の除算はできなくても、200億桁のisqrt()はできるらしい。

あとわしも本文ではいろいろ偉そうに書いているが、頭と終わりの1000桁くらい head & tail したらやっぱり、あと消してしまった。ストレージがもったいないんだもん(笑)。

おまけ: 黄金数100億桁 (頭と尻の100桁) とプログラム

黄金数 $(\sqrt 5 + 1) / 2$ 。

1.
6180339887 4989484820 4586834365 6381177203 0917980576
2862135448 6227052604 6281890244 9707207204 1893911374
... 
5818634698 9002785068 1347162340 2609198093 2411241189
4643236373 7662543421 0515442930 4312640839 3761355984

まあこれもPython電卓みたいなもんですけどね (上は100億2桁計算して2桁切ってます)。計算は19分弱。

import gmpy2
digs = 10000000000
golden = gmpy2.isqrt(5 * gmpy2.mpz(10) ** gmpy2.mpz(digs * 2))
golden = golden + gmpy2.mpz(10) ** gmpy2.mpz(digs)
golden >>= 1
print(golden)

まあそんな小難しい理由ではなく、 $\pi$ は全世界のひとが知ってるくらい (そして無限に続くことも知ってるくらい) 有名だ、ってことかもしれませんけどね。嫁 (※ 哲学科卒) に「きょうさ〜自然対数の底を100億桁求め……」言うたら「なにそれ?」と一蹴された(泣)。

binary splitting法を理解してないまま $e$ が求まった

Chudnovskyの公式を使って $\pi$ を求めた際、binary splitting法というのが出てきた。これは円周率.jpのBinary Splitting法のページによると、いろんな級数の和 ( $e$ とかChudnovskyの $\pi$ の公式とか、 $\sin$ や $\cos$ など広く) を

\displaystyle \frac{P(0, N)}{Q(0, N)} = \sum_{n = 0} ^{N - 1}\frac{a(n)}{b(n)} \frac{p(0)\cdots p(n)}{q(0)\cdots q(n)}

のかたちにしてやる。 $P(a, b)$ とか $Q(a, b)$ を求めたければ、 $a, b$ の真ん中へんの数 $m$ を決めてやって

\begin{array}{rcl} P(a, b) &=& P(a, m) Q(m, b) + T(a, m) P(m, b)\\ Q(a, b) &=& Q(a, m) Q(m, b)\\ T(a, b) &=& T(a, m) T(m, b) \end{array}

と2分割→それを再帰的に計算できる、ってものだ。 $T()$ ってなに!? 級数を『のかたちにしてやる』ってどうやって!?

まあ、ChudnovskyのアルゴリズムについてはC++による実装やBASICによる実装も参考にしたので、な〜んとなく上記が理解……はできてないのにプログラムは書けた。アルゴリズムって凄い。

再帰はちょっとだけ知っている()。前半はまた前半の前半、前半の後半……と分けることにより、最後 $a$ と $b$ の差が1になるまで再帰的に計算できる (論文では5くらいがよいとか)。また前の演算結果に無関係な複数の演算に分割できるので、並列処理向きだ。しかも統合の最後では2つの整数の分数になる。整数になるってことは、固定小数点ならそのまま (10の何乗か倍してからだが) 割り算すれば、固定小数点表現の答が得られる。

そして円周率.jpの一般化された解説ページによれば、いろいろな級数の和に応用できるという。

んじゃ $e$ を求めてみよう。たしか

e = 1 + \frac{1}{1!} + \frac{1}{2!} + \cdots

だから、そんなに難しくないだろう。

と、上記からリンクされている論文 ``Fast multiprecision evaluation of series of rational numbershttps'' をみてたら

$e$ はこの方法が最もかんた〜んに適用できる例であり (※ 超訳)

とかある。が〜ん。

で $e$ の場合は、 $a(n) = 1, b(n) = 1, p(n) = u, q(n) = nv$ ……って $u$ とか $v$ ってなに？と思ったら、これは $\exp(x)$ を求めるときに $x = u / v$ という整数比で表わされるとする、ということらしい。いま求めたいのは $\exp(1)$ だから $u$ も $v$ も1でおっけーで、つまり $q(n) = n$ しか出てこない。たしかにかんた〜ん。

んでもって、binary splittingのプログラム自体は、すでにChudnovskyで円周率で (もっと難しいのを) 作っているので、こんなんでいいの? とか疑いつつ、寝る前に歯を磨きながら (まあPythonのintで固定小数点で) 10分くらいで書きました (ほんとです)。結果……

171828182845...

( $e \simeq 2.71828182845\cdots$ )

っておい、最後に1足すの忘れてる!! ( $N = 0$ から計算すれば必要ないような……まああまり深く考えない)。けど、アルゴリズムって凄い!!

さらに翌日、Chudnovskyで行なったのとまっっったく同じ流れで、固定小数点 (Python int型) 化→gmpy2で高速化→concurrentで並列プロセス化、を30分で書くことができた。

binary splitting法を理解したい

なんでこんな方法で計算できるのか、翌日追いかけてみた。

4項だとわかりやすいだろう。 $(0, 4)$ → $(0, 2)$ と $(2, 4)$ → $(0, 1)$ と $(1, 2)$ 、 $(2, 3)$ と $(3, 4)$ と分割される。ちょうど、上記の $e$ の公式で精度1桁に取ってやると、4項になるので、途中の式にprint()を入れてみた結果

> 0 2 4
> 0 1 2
t = 2 * 1 + 1 = 3
q = 1 * 2 = 2
> 2 3 4
t = 4 * 1 + 1 = 5
q = 3 * 4 = 12
t = 12 * 3 + 5 = 41
q = 2 * 12 = 24
27

つまり $e$ (精度1桁) は $41/24 (+1)$ (3組あるt =、q =の最後) ですよってことなのだが……? 分母の24はわかる。( $1\times 2\times 3\times 4)$ だろ。分子は……? つまり、 $\displaystyle \frac{1}{1} + \frac{1}{(1 \times 2)} + \frac{1}{(1 \times 2 \times 3)} + \frac{1}{(1 \times 2 \times 3 \times 4)}$ のすべての項を $(1\times 2\times 3\times 4)$ で通分すると…… $(2\times 3\times 4) + (3\times 4) + 4 + 1$ だから、『 $(2 \times 1 + 1) \times (3\times 4)$ 』 (さいしょのt掛ける2番目のq) 足すことの『 $(4 \times 1) + 1$ 』 (2番めのt) になる!!

binary splitting法をちょっとだけ理解した

その後もよくわからなかったが、英語版Wikipediaのbinary splittingの説明から参照している解説のページのおかげで、うすぼんやりとわかってきた。ただし、このページ自体は文字が化けている (IBM PCキャラクタグラフィックス???) ので、ページトップにあるPostscirptファイルをダウンロードするとまだわかりやすい。

さてそのPSファイルも、円周率.jpの解説ページ (および元論文) ほど一般化されていないが、もろ $e$ のケースがある。こういうことらしい (項がやや整理されていないので書き直す):

\begin{array}{rcl} Q(a, b) &=& (a + 1)(a + 2)\cdots b,\\ P(a, b) &=& \left\{ \begin{array}{llll} \{(a + 2)\cdot&(a + 3) \cdot \cdots \cdot&(b - 1)\cdot&b\} +\\ \{ &(a + 3) \cdot \cdots \cdot&(b - 1)\cdot&b\} +\\ & \vdots\\ \{ & &(b - 1)\cdot&b\} +\\ \{ & & &b\} + 1 \end{array} \right. \end{array}

上記の計算例の $(a, m, b) = (0, 2, 4)$ のときは

P(a, b) = P(0, 4) = 2\cdot 3\cdot 4 + 3\cdot 4 + 4 + 1

に対して、2分割したものは

\begin{array}{rcl} P(a, m) &=& P(0, 2) &=& 2 + 1\\ P(m, b) &=& P(2, 4) &=& 4 + 1\\ Q(a, m) &=& Q(0, 2) &=& 1\cdot 2\\ Q(m, b) &=& Q(2, 4) &=& 3\cdot 4 \end{array}

より

\begin{array}{rcl} P(a, b) &=& P(a, m) Q(m, b) + P(m, b)\\ &=& (2 + 1) (3\cdot4) + (4 + 1)\\ &=& 2\cdot 3\cdot 4 + 3 \cdot 4 + 4 + 1 \end{array}

おお!!

一般化して、 $P(a, m) Q(m, b) + P(m, b)$ は

\left\{ \begin{array}{llll} \{(a + 2)\cdot&(a + 3) \cdot \cdots \cdot&(m - 1)\cdot&m\} +\\ \{ &(a + 3) \cdot \cdots \cdot&(m - 1)\cdot&m\} +\\ & \vdots\\ \{ & &(m - 1)\cdot&m\} +\\ \{ & & &m\} + 1 \end{array} \right\}\\ \times\\ (m + 1)(m + 2)\cdots b\\ + P(m, b)\\ =\\ \left\{ \begin{array}{llll} \{(a + 2)\cdot&(a + 3) \cdot \cdots \cdot&(b - 1)\cdot&b\} +\\ \{ &(a + 3) \cdot \cdots \cdot&(b - 1)\cdot&b\} +\\ & \vdots\\ \{ & &(b - 1)\cdot&b\} +\\ \{ & & &b\} + (m + 1)\cdots b \end{array} \right\}\\ +\\ \left\{ \begin{array}{llll} \{(m + 2)\cdot&(m + 3) \cdot \cdots \cdot&(b - 1)\cdot&b\} +\\ \{ &(m + 3) \cdot \cdots \cdot&(b - 1)\cdot&b\} +\\ & \vdots\\ \{ & &(b - 1)\cdot&b\} +\\ \{ & & &b\} + (b + 1)\cdots b\\ \end{array} \right\}\\ =\\ \left\{ \begin{array}{llll} \{(a + 2)\cdot&(a + 3) \cdot \cdots \cdot&(b - 1)\cdot&b\} +\\ \{ &(a + 3) \cdot \cdots \cdot&(b - 1)\cdot&b\} +\\ & \vdots\\ \{ & &(b - 1)\cdot&b\} +\\ \{ & & &b\} + 1 \end{array} \right\}\\ = P(a, b) (ふぅ)

となるわけだ!!

誰が考えたんでしょうね〜頭良い!!

ちなみにまだわたしが理解できていないのは

より一般化されたケース ( $T$ が入ったり、さらに原論文で $B$ とか $C, D$ が入っている場合)
$e$ みたいな簡単なケースでなくて、どーやって元の級数を $P$ とか $Q$ とか ( $A$ とか $B$ とか) にするのか

あたりである。

何項まで計算するか

さて、 $e$ は超越しているので(笑)、10000桁くらいまでの答え合わせは簡単にできる。gmpy2のisqrt()は精度が良いので、余分なエラー桁をみなくてもよさそうだったが、まあ2桁くらい取っている。

さて100桁くらいまで正しいことを確認して、あれ? 100項まで計算してるけど、こんなに必要? とか気づく。

試しに80項くらいに落としても……100桁目まで合っている。50項だとダメ。試すと70項までは落とせるようだ。

まあこれは $1/x!$ が有効桁数内で0、つまり $x!$ の桁数が求めたい桁数を超えたときに計算を打ち切ればよいのだが、binary splittingをするにあたっては初めから項数を決める必要がある。 $x!$ の桁数って? ぐぐったら簡単にみつかった。インターネット便利(笑)。

だがしかし、dig $=\displaystyle x\log_{10}\frac{x}{e}$ を $x$ について解く、なんていう数学的素養は持ち合わせていない。おそるおそる、ニュートン法で解くプログラムを作ってみると……ちゃんと数回の繰り返しで求まる!! アルゴリズムって、凄い!!

x!の桁数 (プログラムリスト)

# digit of factorial:
#   dig == floor(x log10(x/e))
# get x for given dig by newton method
#   f(x) = x log10(x/e) - dig
#   => x_{n+1} = x_n - (x_n log10(x_n / e) - dig) / (log10(x_n / e) + 1)
#      (x_1 = dig)

def digfact(dig):
    xn = dig
    epsilon = 10 ** (-math.log10(dig) - 2)
    print(epsilon, file = sys.stderr)
    while True:
        xerr = (xn * math.log10(xn / math.e) - dig) / \
            (math.log10(xn / math.e) + 1)
#        print(xn, xerr, xn * math.log10(xn / math.e))
        if abs(xerr / xn) < epsilon:
            break
        xn = xn - xerr
    return int(math.floor(xn))

epsilonは求めたい $e$ の桁数に応じて設定しているが、まあこれも誤差を見込んで2桁くらい余計に求めることにしている。おいおい、 $e$ の値を求めるために、 $e$ の値が必要だよ!!(笑)

これで計算時間の最小化ができ、余計な桁数の計算によるメモリ不足は抑えられる、たぶん。知らんけど。

$e$ の100億桁を求める

時間は前後する。 $\pi$ を計算したときハマったのは

計算のかなりの部分は再帰でマルチプロセスに分けられるから速くなる……が、最後に整数2個に統合されたところでは親プロセスが延々、超多桁の割り算をするので、結局時間がかかる
最後の整数の割り算は、たとえば (200億桁 ÷ 100億桁) を求めようとしたら、gmpy2では50GBオーバーのメモリ (VM) が必要→AWSを借りれば金で解決
そのあと (これもまたシングルプロセスで) 延々1時間掛けて、100億桁を10進数に変換

といったあたりで、今回の自然対数の底でも同じことでハマる (というか、すでに回避することは諦めて、同じ罠にハマることにする)。

eをbinary splitting法で求める (プログラムリスト)

import math
import gmpy2
from gmpy2 import mpz
from concurrent.futures import ProcessPoolExecutor as Executor
import sys


errdig = 2
#errdig = 0
maxproclev = 3
#maxproclev = -1

# (x!の桁数を求める digfact() は前出)

#
# get e by binary splitting method
#
# at exp(1) = \sum_x (1 / x!),
# u = v = 1, P(n) = 1, Q(n) = n, B(n) = 1
#   => Q = Ql Qr = Q(n1, m) Q(m, n2),
#      T = Br Qr Tl + Bl Pl Tr = Q(m, n2) * T(n1, m) + T(m, n2)
#
# シングルプロセス版
def ebsiter(n1, n2):
    if n2 - n1 == 1:
        t = 1
        q = n2
    else:
        m = (n1 + n2) >> 1
#        print('>', n1, m, n2)
        t1, q1 = ebsiter(n1, m)
        t2, q2 = ebsiter(m, n2)
        t = q2 * t1 + t2
        q = q1 * q2
#        print('t =', q2, '*', t1, '+', t2, '=', t)
#        print('q =', q1, '*', q2, '=', q)
    return t, q

# マルチプロセス版
def ebsitermp(n1, n2, lev):
#    print('>', n1, m, n2)
    if n2 - n1 == 1:
        t = 1
        q = n2
    else:
        m = (n1 + n2) >> 1
        if (lev <= maxproclev):
            with Executor() as exec:
#                future1 = exec.submit(ebsitermp, n1, m, lev + 1)
                future2 = exec.submit(ebsitermp, m, n2, lev + 1)
                print('exec future', lev + 1, '(2)', file = sys.stderr)
#                t1, q1 = future1.result()
                print('exec future', lev + 1, '(1)', file = sys.stderr)
                t1, q1 = ebsitermp(n1, m, lev + 1)
                t2, q2 = future2.result()
                print('done future', lev + 1, file = sys.stderr)
        else:
            t1, q1 = ebsiter(n1, m)
            t2, q2 = ebsiter(m, n2)
        t = q2 * t1 + t2
        del t1, t2
        q = q1 * q2
        del q1, q2
    return t, q


def ebs(dig):
    prec = dig + errdig
    nterm = digfact(prec)
    print('prec', prec, 'digs, nterm', nterm, file = sys.stderr)
## Python int version, single process
#    t, q = ebsiter(0, nterm)
#    base = 10 ** prec
## gmpy2 version, single process
#    t, q = ebsiter(mpz(0), mpz(nterm))
#    base = mpz(10) ** mpz(prec)
## gmpy2 version, multiprocess
    t, q = ebsitermp(mpz(0), mpz(nterm), 0)
    base = mpz(10) ** mpz(prec)
#    print(t, q)
    return base + base * t // q


if __name__ == '__main__':
    if len(sys.argv) != 2:
        print('usage: ebs DIGITS', file = sys.stderr)
        sys.exit(1)
    dig = int(sys.argv[1])
    e = ebs(dig)
    print(e)

ebs()のコメントアウトで、Python intのシングルプロセス版・gmpy2のシングルプロセスおよびマルチプロセス版を切り替えられる。

ebsiter()、ebsitermp()に与える初期値をmpz()で囲ってやるかどうかだけで、Python int型かgmpy2.mpz型かが決定される ←サブルーチン内の加算乗算で、一方がmpz型だと結果もmpz型なので、演算全体にmpz型が波及する。便利。

桁数と計算時間は、こんな感じになる。

紫: Macbook Pro 16プロセス並列、緑: 同8プロセス並列、青: Macbook Air並列なし、橙: AWS g3.4xlarge 16プロセス並列 (以上gmpy2)、黄: MBA Python int

基本的なPython int (Karatsuba乗算) に比べてgmpy2 (FMT/NTT乗算) が速いのはいまだに驚くが、今回の $e$ のbinary splitting法では、並列数を増やしたことによる時間の減少が少ない……というか、Macbook Proは10コアなので多めの16並列にするとムダがなかろ、と思っているのだが、どういうわけだかかなりの桁数までは同じくらいである。例によって、半額以下のMacbook Airとの計算そのものの速度差はあまりないが、したがってコア数が多少多いのも活かせない。メモリが8GBと16GBの差は、桁数の限界に影響するだろうが……それとてProでも足りなくてAWSを借りるんだから(次節)、あまり関係ない(笑)。

またまたAWSを借りて計算

Chudnovskyで $\pi$ を求めたときにすでに、(200億桁 ÷ 100億桁) の割り算は、手許のMacbook ProではVM不足によりできない、ということは分かっていた。binary splitting法ではどうあがいても、最後は100億桁の整数どうしの割り算になるので、同じところで落ちるはずだ。試したら、やはり25分 (10コアで16プロセス並列) で落ちた。

今回はAWSでtopを観ていたので、除算部分の最高は62.5GB程度のメモリを喰っているようだった (アーキテクチャは異なるがGMPの内部表現なのだからほぼ同じくらいのメモリを喰っている)。

さて前回AWSを借りたときは頭に来ててろくに選びもせず、とあるプランを借りたが、今回はいろいろ吟味して

オハイオでなく東京ノードで借りる(笑)
メモリは例によって128GBあれば足りるだろう
コア数は増えてもあまり関係ない (最後のシングルコア部分がボトルネック)

ということで、r6g.4xlarge (16vCPU、128GiB) というのを借りてみた。USD0.9728/時間、激安!!

このプランではアーキテクチャがx86ではなくARMなので、gmpy2が動くまでに余計な手順が必要だった。

AWS (x86、ARM) でgmpy2を動かすまで

前記事にも書いたが、まずはpipがないので

wget https://bootstrap.pypa.io/get-pip.py
python3 get-pip.py
export PATH=$PATH:~/.local/bin

が必要である (https://docs.aws.amazon.com/elasticbeanstalk/latest/dg/eb-cli3-install-linux.html)。さらにARMでは

error: command 'x86_64-linux-gnu-gcc' failed with exit status 1

とか出るので

sudo apt-get install build-essential libssl-dev libffi-dev python3-dev
sudo apt install libmpc-dev

が必要である。

その後

pip3 install gmpy2

でgmpy2がめでたく入った。

ところが。少ない桁数でテストしてみたら、これが遅いのなんの(笑)。1億桁で試したら、Macbook Air/Proの5倍くらい遅かった。さらに桁数が増えたときの時間の増加も、なんか知らんが凄い。

こちらはさっそく停止して、『高速コンピューティング』プラン g3.4xlargeを借りてみることにした。上記とほぼ同じスペックの16vCPU・122GiBだが、USD1.58/時間である。さっそくインスタンスを作成しようとしたら……『vCPU数が制限を超えている』旨のエラーが。こないだChudnovskyで円周率計算をしたときは64vCPU使えたのに?

ダッシュボードの『制限』でvCPUで検索してみると、『standard (A, C, D, H, I, M, R, T, Z)』では64vCPUとなっているが、それ以外は0vCPUである。gで始まっている『高速コンピューティング』には制限があるようだ。知らんけど。

手順に従ってGタイプのvCPU制限を16にしてくれ、とリクエストを出してみる。人間が処理しているようで、サポートにメッセージが送られるようだ。知らんけど。

使用目的なんぞ書かないといけないので、てけとーに『to execute more mathematical simulation』とか書いてみる。なんだろ、原爆とか人類に反乱するAIでも設計するとか思われてるのかな? (笑)

こりゃ実行は明日かな〜、と思っていたが (※ 実は時間順序が逆だが、Gタイプが借りられないので申請してから上記の超遅い激安Rタイプを借りてみている。Rが遅すぎて諦めて、ふとダッシュボードを見ると) 30分ちょっと? でリクエストが通っている。深夜に働いているAWS使用目的審査員とかいるんかな? 乙です。知らんけど。

ということでさくっとマシンを作ってプログラムを流し、テストランしてみた。まあ1億桁で (例によって) Macbook Proの2倍くらい遅いけど、まあ許容範囲だ。体感では、Chudnovskyで5時間弱で $\pi$ を求めたとき (オハイオのc6a.16xlarge・USD2.448/時間) と同じくらいな気がする (コア数が違うけど、1コアのボトルネック部分の計算速度とか)。

計算時間は3時間19分弱。さすがに並列部分は50分ほど? で終わり、 $\pi$ より ~~超越度が小さい~~ 計算量が少ないんだな、とか感じる。例によって除算→2進10進変換・表示部分は同じだけ時間がかかって、うっかり寝る前に実行したのだが、終わらないので寝る。たまたまトイレに起きたら実行が終わっていたので、即止める(笑)。

ダウンロードは翌日やったのだが (実は失敗して2度ダウンロードしている)、オハイオの1時間半よりやはり東京は断然速く、30分弱で10GBほどのファイルは落とせた。しめて4時間弱 (寝たり失敗したりは除く) で、ちょっとだけ円高になってきたし800円弱!! やっぱり $\pi$ (1700円)より $e$ の方が安い!! (※ 慣れてきたからだろ)

AWSのプランは、もう少し使用目的←→性能←→料金についてのおべんきょが必要そうである。このあたりをきちんと知っているのが『AWSの資格』ってことか?

まあ世の中、AWSはレンタルサーバとして、低性能でいいから24時間稼働する、みたいな使われ方が多く、わしの今回一連のような、自分でマシン組むのが面倒 & わずかのピーク性能だけ欲しいので、『昔の計算機センタ』のような使い方 (秒いくらでぶん回す) というのはマイナーなのかな、とも思うが (しかしdeep learning流行りなのでそういう使い方も多いような気もするが)。

AWSのプランはよく出来ている。金を積むと計算が速く終る。

ま〜結局、世の中金ですわ。

結論: eの100億桁目は

6です。

自然対数の底 100億桁 (頭と尻の100桁)

e=2.
7182818284 5904523536 0287471352 6624977572 4709369995
9574966967 6277240766 3035354759 4571382178 5251664274
... 
3067538682 6276859582 4437326336 8509746850 9182383460
0420394247 5819248992 9739871342 5699481693 1502679526