Zenn
🧮

10進100億桁の円周率を2進数にしてみようと思ったはなし

2022/11/29に公開
Python
アルゴリズム
円周率
2進数
10進数
tech

また円周率かよ。

ってかあまり円周率には関係ない (かもしれない) ぷちネタであるが。多倍長精度固定小数点の10進数を2進数に変換しようとしてハマったはなしです。

2進表現の円周率がほしい (べつに要らんけど)

円周率100億桁を自力で求めてみた。まあこれとて、Python + gmpy2の無限多倍長精度の整数演算を使って、固定小数点的に求めてみたわけであるが……最後の結果の確認をするためにも、そして自己満足喜びに浸るためにも、10進数での表現が必須である。すなわち『3.141592...』とあるとひとは誰でも「あ、円周率だ!!」と言ってくれるが、『0x3.243f6a...』を見せられても「なにそれ?」となるだけだ。かっこいいけどね。

んだけど、『円周率100億桁 (10進数)を出力して見てみたい』のネタのときに、その膨大さをひとに話すと、幾人かのひとに「それ、2進数にしたら?」と言われた。

われわれの手が10本指だから10進数が普及したらしいが、10進数が計算機や数千枚のプリンタ用紙の上で冗長なのは、計算機屋なら知っている。10進数1桁の情報量は、2進数にするとおよそ 1 / \log_{10}2 \simeq 3.32ビット程度となる。円周率は乱数と同じなので1ビットのエントロピは1ビットあるとみなせる。(10進) 100億桁のエントロピは332億ビット程度、単純に10進数100億桁をテキストファイル (1桁 == 1文字 == 8ビット) に格納 (10GB) しているとすると、2進数表現なら4.15GB程度になるわけだ。2進数プリントアウトチャレンジのところで書いたように、10進数1桁 (3.3ビット) の1文字が6 × 4ドット (24個) の□■で表現されているのに比べたら、2進数1ビットを□■で表現したら紙の使用量は1/7程度になる(笑)。

面白いことに、この冗長な10GBのテキストファイルをbzip2で圧縮すると4.31GB程度になる。bzip2くんはテキストファイルとか数値表現とか知らないのに、ほぼ情報理論的限界、エントロピまで圧縮してくれている、といえる。

まあ2進数の粒粒で表示した円周率なんて、つまり『\pi = 11.0010010000111111...』なんてのが並んでいても、たとえ0が□で1が■だとか、0が♡で1が👍になってるとかでなかったとしても、ひとはこれを見て円周率だ、と思わないだろう。紙に疑似ランダムな□■をプリントしたとしてもただの砂嵐だし (いや極小文字をプリントアウトしてもほぼ砂嵐だったけど)。まあオブジェとしてはいいのかもしれないけど。

んでもって自分としては、円周率の2進表現なんてとりあえず要らないし、そもそもgmpy2の内部では2進数で計算している (ものを、計算時間の3割くらいかけて最後に10進数変換→プリントアウトしている) んだから、そんなものいつでも取り出せるだろ? とか思っていた。

甘かった。固定小数点なのだった。
※ そもそもmpz型は手続き.to_byte()を持っていないので、『いつでも』バイナリは取り出せない。

つまり「2進表現が欲しければ、\piを求めるプログラムそのもののメモリダンプを追加して計算し直せばいいだろ……」くらいにしか考えていなかったが、そのままメモリダンプしたところでその2進数列は小数部の数列になっておらず、あくまで小数部を10の何乗か倍した整数なのだった (つまり元の\piを求めるプログラムを改造するところで、以下の記事のルーチンを追加しなければいけない。結果のテキストファイルを読み込んで変換するのと、何ら変わらない)。

例えば『3.1415』の小数部の『0.1415』は、精度4桁 (10^4ベース) の固定小数点だと『1415』という整数になる。これは16進数 (2進数)でいうと『0x587 (0b101_1000_0111)』であって、0.1415の2進数表現『0x243(c) (0b0010_0100_0011_11...)』 (※ 10進4桁は2進約13ビットなので、14ビット目は正確ではない) とは似ても似つかない。あれ、ちょっと似てるか(笑)←偶然です。

なぜ似ても似つかないかというと。2進数の小数部はどう表現されるかというのを考えるとわかる。小学校で教わった10進→n進変換の手順を忠実になぞると、上記の例を2進数に変換するには、次のような算法を取る (実はこの算法は、まずい(後述))。

アルゴリズム1.

  1. 0.1415をベース10000でa = 1415とする。最初にc = 5000と取る (これは下1桁目の重み1/2^1をベース倍したものである)
    1. aとcを比較する。aのほうが小さいので1ビット目は0
    2. a \leftarrow a - c \times 0 == 1415のまま
    3. c \leftarrow c / 2でc == 2500になる (次の重み1/2^2=0.25)
    1. aとcを比較する。aが小さいので2ビット目も0
    2. a \leftarrow a - c \times 0 == 1415のまま
    3. c \leftarrow c / 2でc == 1250になる (次の重み1/2^3=0.125)
    1. aとcを比較する。aが大きいので3ビット目は1
    2. a \leftarrow a - c \times 1 == 165
    3. c \leftarrow c / 2でc == 625になる (次の重み1/2^4=0.0625)
  5. ...同様に、4桁の10進数のビット数を超える\lceil 4 / \log_{10}2 \rceil = 14ビット目まで繰り返す (cが0になるので計算停止)

これをコーディングしてしまったばかりに、わしのぷち闘いが始まる。

単純な割り算法

と、その前に、10進数表現された円周率100億桁のテキストファイル (全桁つながっており改行などなし) を変数aに読み込む部分 (以下共通) は、こんなんでできる。

import gmpy2
import sys

with open(INFILE, 'r') as fp:
    s = fp.readline()
ia = int(s[0])
s = s[1:]
a = gmpy2.xmpz(s)
del s

実際のプログラムでは引数でファイル (テキストの入力INFILEとバイナリの出力ファイル) を指定できるようにしていたり、デバッグのため桁数を指定してfp.read()で桁数を切って読めるようにしている。

ファイルの1バイト目は、まあ円周率だから『3』に決まっているが(笑)、真面目にiaに代入して削除している。gmpy2.[x]mpz()は文字列 (10進数) を直接、[x]mpz型に変換できる。すかさずsを削除しているのは、メモリの節約のためだ。

ちなみにこれを書いたときはす〜っかり忘却の彼方だったのだが、100億桁は手許のラップトップ (Macbook Pro) ではメモリが足りず、AWSの計算機を借りて求めていたのだった(笑)。あれ、読み込めるかな? と思っているうちに、案の定落ちた。他のソフト (Chromeとかすべて) 落としてシェル類だけにして動かしたら、なんとか最大VM使用量43GB (前の記事に書いたが、これが50GBを超えると、OSだかPythonだかに叩き殺される) でぎりぎり読み込めるようだった。なお、この読み込み部分の実行だけに25分掛かる(笑)。

さて、前節のアルゴリズム1. の変換部分は以下のようになる。

c = gmpy2.xmpz(10) ** gmpy2.mpz(gmpy2.num_digits(a))
c >>= 1
with open(OUTFILE, 'wb') as fp:
    mask = 0x80
    b = 0
    fp.write(ia.to_bytes(1, 'big'))
    while True:
        if c <= a:
            b = b | mask
            a -= c
        c >>= 1
        mask >>= 1
        if mask == 0:
            fp.write(b.to_bytes(1, 'big'))
            b = 0
            mask = 0x80
            limitbytes = limitbytes - 1
            if c == 0:
                break

※ 本文あちこち、mpz型と書いていますが、演算代入置換したいのでmutableなxmpz型を使っているものもあります。まあ演算過程ではどちらでも同じくらいメモリ喰ってると思うけど (知らんけど)。

で、これ何がまずいかお分かりでしょうか。2つのまずい点があります。

まずひとつは……遅い(笑)。何も考えず上記を動かしたら、上記部分はまあ最大メモリ使用量10GBくらいで余裕で動くのだが、(100億桁 - 100億桁) の引き算 (比較、if c <= a:部分) を1ビットずつ行なうわけですよ。これを332億回 (笑)。翌朝起きて (日曜朝寝したので12時間後くらい)、まだ数KBくらいしかできていない。完成予想サイズ4.15GB分の数KB。つか、ファイルのバッファリングのせいで出力0バイトで、もーいーやと思って止めたから数KB吐き出したんだが、桁数で時間を割ってみたら2.7秒で1回の割でとろとろループを回っている。完成するのに2800年(爆笑)。いくら引き算は乗算・除算に比べ (100億桁だと乗算5分、除算20分のレベル) 速いとはいえ……うかつでした。

なお、上記はcとaの桁数がどんどん減っていくので、どんどん (少ない桁数だとね) 速くなっていく。たぶん1000年くらいで終わるね。それから、ビット演算系はさすがに一瞬で終わる。ただしgmpy2 (Python3) は最適化によるビット演算への読み替えなどはしていないようで、例えばa <<= 3と書くと一瞬だがa *= 8と書くと (乗数が一桁なので、被乗数が100億桁あっても秒で終わるが) オーダとして3桁くらい遅い。

もうひとつまずいのは、そもそも前節のアルゴリズムがまずいのだが。100億桁トライする前に100桁〜10000桁くらいで試してみても、最後 (下) の数ビット〜数バイトが正解とは似ても似つかぬ数字になってしまう。演算誤差がある、っていうことだ。ちなみに一晩で数KBしかできていなかったときは、実はこのエラーは分かっているけど解消の方法がわからなかったので寝たまでである。いずれ止めるの前提で動かしていたので。

どこがまずいかというと。前節アルゴリズムでc (ビット重み) を次々2で割っていくところで、cが整数だと2で割り切れない数になってしまう、ということである。前節の手順の続きを書くと、ステップ7.での除数は312.5 (0.03125)にならなければならないのに、整数なので312になってしまう。

一般化すると、10進n桁精度の固定小数点ではc = 10 ^{n}になるから、n回までは2で割り切れる。だがその後は5のべき乗になってしまうので、絶対に2で割り切れない。除算は (というかループを回る回数は) およそ10進n桁の場合、n / \log_{10}2 \simeq 3.3n回くらい行なうので、後ろの2.3n回くらいのループでは誤差がある数を除数にするので、下の方の桁には結果の誤差が積もり積もってしまう。

まあここが、10進←→2進変換の面倒なツボなのだが、いやはや。人類の手の指が5 (素数) 本だったばかりに……このアプローチではダメである。

加速しないけど、正確

前節の二つの問題点のうち後者、精度に関しては、割と簡単に解決できた。cを1/2に小さくするのでなくaを2倍に大きくすればよい。1が頭にあると桁数が増えそうだがそれは一瞬で、どうせcを引き算したら元の桁数に戻る。つまり

アルゴリズム2.

  1. a = 1415をベース10000で変換する。最初にc = 10000 (1.0)と取る。以後cは不変。
    1. aを2倍する。a \leftarrow 2830
    2. aとcを比較する。aのほうが小さいので1ビット目は0
    3. a \leftarrow a - c \times 0 == 2830のまま
    1. aを2倍する。a \leftarrow 5660
    2. aとcを比較する。aのほうが小さいので2ビット目は0
    3. a \leftarrow a - c \times 0 == 5660のまま
    1. aを2倍する。a \leftarrow 11320
    2. aとcを比較する。aのほうが大きいので3ビット目は1
    3. a \leftarrow a - c \times 1 == 1320
  5. ...同様に、4桁の10進数のビット数を超える\lceil 4 / \log_{10}2 \rceil = 14ビット目まで繰り返す

なんだ、簡単じゃん。プログラムとしては、c >>= 1の行を

        a <<= 1

に替えるだけである。

一瞬、cが10000 (10のべき乗) というキリ番きりの良い数で不変なので、「あれ引き算 (比較) しなくても、ビット調べれば大小関係わかるじゃん!?」と勘違いしそうになったが、10進数できりが良くてもコンピュータ内部は2進数です。たとえば10000は0b10_0111_0001_0000っつー全然きりの良くない数字ですんで、あしからず。

さてgmpy2的 (可変の多倍長整数演算的) には、前述のアルゴリズムでは桁数が小さくなるのでどんどん速くなっていくが、このアルゴリズムではaもcも桁数不変なので一定速度である。1ビットずつの変換をしていたら……前のより時間が掛かってしまう!!

桁数を減らせ

上記のふたつのアルゴリズムは、プログラムにすれば1行しか違わない。だから、簡単にスイッチすることができる。つまり、最初のn回の割り算まではcを2で割っていき、それ以降はaを2倍する。
まあ桁数はスタートの7割くらいまでしか少なくならないが、100億桁の除算・余りを求めるのと70億桁のとでは大違いである。

ループ部分は

    diga = gmpy2.num_digits(a)
    nloop = int(math.ceil(diga / math.log10(2)))
    ...
    keepprec = diga
    for i in range(nloop):
# since (initial) c == 10 ** diga == (2 ** diga) * (5 ** diga),
# shifting c right less than or equals to diga bits
#   (c >> diga) [bits] doesn't lose precision
        if 0 < keepprec:
            c >>= 1
            keepprec -= 1
        else:
            a <<= 1

こんな感じである。

まとめてどん

ビット数を多少減らしたところで計算時間は劇的には変わらない。いや違った、計算時間は2800年が1400年くらいに経ると思うので劇的に変わるが、計算量オーダは変わらない。

んでここで、多桁のまとめ演算を導入する。何ビットか (ここではsとしよう) まとめて2進変換してしまうのだ。これは考え方としては最初の (誤差ダメダメ) アルゴリズムのほうが考えやすい (ただし固定小数点 == 整数であることは、以下ちょっと忘れる)。

1ビットずつ変換する方法では、さいしょに10000 (1.0) の半分の5000 (0.5) がaにいくつ (0または1個) 入るか、を考えていた。ここで例えば4ビットまとめて、625 (1/2^4 = 0.0625) がいくつ (0〜15個) 入るかを考える。最初のaは1415だから2 (0b0010) 個入る。以下 (最初から書くと)

アルゴリズム3.

  1. a = 1415をベース10000で変換する。最初にc = 10000と取る
    1. cを16 (=2^s)で割る。c \leftarrow 625 (0.0625)
    2. (d, r) \leftarrow (\lfloor a / c \rfloor, a \% c)、つまり商と余りを計算する。 d = 2 = 0b0010、 r = 165
    3. a \leftarrow r = 165 (0.0165)
    1. c \leftarrow c / 16 = 39.0625 (0.00390625、ここでは誤差を避けるため実数とする)
      1. (d, r) \leftarrow (\lfloor a / c \rfloor, a \% c) から、d = 4 = 0b0100、 r = 8.75
    3. a \leftarrow r = 8.75 (0.000875)
  4. ...同様に、4桁の10進数のビット数を超える14ビット目を含む4ブロック目まで繰り返す

上記手順3.まででも、2回のループで4 × 2ビット、小数点以下8ビットの0b0010_0100が得られる。このブロックサイズ s を大きくすれば、ループ回数を減らせるはずである。

gmpy2.mpz型を除算して得たs桁の2進ブロックdは、Pythonのint型に変換する (そうしないと.to_byte()できないのでバイナリファイルに書き出せない)。

それから上記のアルゴリズム3. のまま整数演算 (アルゴリズム1. 同様) にするとcを割っていくときの誤差が生じるので、実際はアルゴリズム2.同様にaを大きくするのだが、(ブロックサイズ × ループ回数) が最初のaの桁数 (最初のcに2べきを含む個数) を超えない限り、aを大きくするのではなくcを小さくすることができる (つまりやはり、3割くらい桁数は減らして演算できる)。

本番

実際に本番プログラムを動かしたときは、n = 10進100億桁 (2進約332億ビット) のaに対してブロックサイズはs = 80億ビット == 10億バイト に取ってやり、最初のループで一気に100億桁のa、割ることの76億桁のcまで減らした (cの素因数には100億個の2べきが含まれているので、1ブロック目80億ビット \simeq 10進24億桁ぶんは問題なく減らせて、76億桁となる) 。てきとーな暗算で決め打ちして様子をみたのだが、思ったように桁数が減り、思ったより速く最初の5分 (後述) くらいで落ちずに割り算ができたっぽいので、そのまま続行した。実は桁数を減らすことより、100億桁に対して24億桁 (80億ビット分) ずつ商を求めているのでループ5回目は24億桁のうち4億桁ほどしか有効な数字がなく、もう少し増やせば4回ループで演算が終わっている。

除数cの桁数は徐々に減らすか、きちんと限界を狙って計算 (sがぎりぎり100億弱 == 12.5億バイト弱) すれば、前述の通り 100億 \times (1 - \log_{10}2) \simeq 69.9億桁まで減らせる。これを超えたブロックサイズを設定すると、(1ブロック目でcに含まれる2のべき乗数がnを超えることはできないので) cの右シフトは1回も起こらず桁数が減らなくなってしまう。

これらより、桁数に関わらず s = n / 8 (よりぎりぎり少なく) に取ってやれば、最大の桁数減になり、またループは4回で終わる (\lceil 1 / \log_{10}2 \rceil = 4より)。メモリさえ許せば。

前記事にあるとおり、Macbook Proの環境では120億桁 ÷ 100億桁までの除算は落ちずにできることがわかっているので100億桁 ÷ 76億桁は割と余裕でできる。その後、100億桁 ÷ 76億桁の5回のループで約332億ビット目を含むブロックまで計算した。

計算時間は72分、前にも書いたようにうち25分は文字列の10進100億桁を読んでmpz型に変換する時間、あとのループ5回は1回につき10分程度で回っている。メモリ・計算速度に余裕があるので面倒だから商と余りは別々に書いたが (除算より乗算のほうが軽いので、剰余演算は減算と乗算で書いたほうが速い)、1回10分のループの大半は、100億桁 ÷ 76億桁の除算1回と剰余演算1回 (つまり除算2回) のはずである。100億桁 ÷ 70億桁の4回ループなら、もう少し速くて1時間を切るくらいかもしれなかった。

本番プログラムでは、最後のブロックの無効部分を削除、そして最後のバイトの (精度以下の) 下位ビットの除去 (0マスク) をするところまで書いたが、16進で検算可能ないろいろな桁数 (16進で 〜10000桁) で試してみると、どうもやはり最後の2ビットくらいは演算誤差が出るようである (おそらく5で除算するときの誤差1回分くらい?)。

かくして、小数点以下10進100億桁 (全桁正確) から、2進33219280949ビット程度、バイナリファイルにして4152410118バイトと左の5ビット (うち2ビットは怪しい)の2進数列ファイルが出来上がった。

結論： だからど〜した

かというと、こんな成果に花開いた (※ 開いてない)。

書いてて、その昔『BCDの加算は6をど〜たら』というのがあったのを思い出した。同じようなものだが、プログラミングをたしなむひとのうち、ほんの一部のひとしか出会わない、ムダニッチな技術かもしれない。

Discussion