GHCのIntegerの深淵
前書き
先日の日記のSPOJにおける√2を200万桁求めるという課題で、私はCを用いてそれを達成したのであるが、Haskellで100万桁を求めている人がいて、それがずっと気になっていた。普通に計算すると10万桁で精一杯で、とてもHaskellで100万桁なんて無理のように見える。ところが、ふとしたことから、Haskellで100万桁が可能であるように思えてきて、実際にやってみたら出来てしまった。しかし、そのためにはGHCのIntegerに対する深い理解が必要なのであった。
戦績：
http://www.spoj.pl/ranks/SQRT2/ Haskellで100万桁
http://www.spoj.pl/ranks/EVAL/ ついでにeも100万桁
http://www.spoj.pl/ranks/PIVAL/ さらについでの円周率50万桁
GHCのIntegerにおける乗算
GHCのIntegerはFFTで乗算がなされている。GMPというライブラリがGHCに同梱されているのでおそらくこれが用いられていると思われる。実際のところ乗算はかなり高速に実行される。ところが、GHCに適当に√2を求めるプログラムを書いて入れると余り速くない。また、GHC6.6と6.4ではかなり速度差が見られる。これはどういうことなのだろうか。
show
計算結果を表示させないようにすると、GHC6.6と6.4の速度はほぼ同じになる。つまり、6.4ではshowが遅いということである。なぜshowが遅いのだろう？多倍長整数を表示する処理は、内部表現に依存する。内部表現が10進数なっておれば、その内容を10進数で表示するには、桁ごとにそのまま単に表示していけば良いだけだ。内部表現が16進数になっている場合、10進数で表示するには基数の変換が必要になる。そして、GHCのIntegerにおいては内部表現は16進数になっている。基数の変換は10(あるいはもう少し大きい10の累乗)での単精度除算を繰り返すことによりO(n^2)で行うことが出来る。また、Karatsuba基数変換アルゴリズムなどを用いればO(nlogn)で行うことも出来る。GHC6.4では(おそらく)O(n^2)の実装になっており、桁数が大きくなるとshowは使い物にならなくなる。6.6ではO(nlogn)の実装がHaskellの上でなされており、それなりの速度で変換することが出来る(とはいえそんなに速くない)。
shift
Integerは16進数で保持されている。さらに、IntegerはData.Bitsのインスタンスである。よって、IntegerはBitsのメソッドshiftを用いることによりロジカルシフトが可能なはずなのだが、なぜだか分からないがIntegerに対するshiftはGHC6.4、6.6ともに *(2^x)、`div`(2^x)、によって実装されている。つまり、2の累乗による乗除算はO(n)で出来るはずなのに、GHCではO(nlogn)かかってしまう。
unboxed type
Haskellで多倍長演算を行うには、上記のような問題点を解決しなければならない。そのためには、C言語でIntegerの中身をいじるルーチンを書くか、あるいはHaskellからIntegerの中のデータをいじるかの必要がある。今回の場合、SPOJに通したかったので、Haskell内で完結させる必要がある。HaskellからIntegerのデータをいじる方法をとることにする。
そのためにはunboxed typeというものを使う。GHC.ExtsかGHC.Primあたりのモジュールをimportして-fglasgow-extsオプションを渡すと使うことが出来る。Int#などサフィックスに#が付いた型がそれにあたる。普通のHaskellの値は型情報の付いたオブジェクトになっており、整数などの本来1ワードであらわせるデータも例外でない。unboxed typeが示すのはそのような情報が付随しないデータであり、普通の型とは種も異なっている。プリミティブな型の値を操作するためにこのunboxedな型を用いることが出来る。
data Integer =
S# Int#
| J# Int# ByteArray#Integerのunboxed typeはこのように定義されている。S#は1ワードに収まる範囲の数を表し、J#は多倍長整数を表している。J#の一つ目の引数は桁数と符号をあらわしており、二つ目のByteArray#型の引数が整数のデータにあたる。桁数はByteArray#中のデータのワード数*符号(1 or -1)が格納される。ByteArray#は単なるメモリブロックへのポインタのようなものである。これに対してはHaskellからデータの読み取りをすることが出来る。
効率のよいshiftの実装
intlen = 4 -- 32ビットCPUなら4
extend (J# n# d#) (I# a#) =
case newByteArray# ((n# +# a#) *# intlen#) realWorld# of { (# s#, e# #) ->
case f# e# 0# s# of { s# ->
case unsafeFreezeByteArray# e# s# of { (# s#, ret# #) ->
J# (n# +# a#) ret# }}}
where
f# r# i# s#
| i# ==# (n# +# a#) = s#
| i# >=# a# =
case writeIntArray# r# i# (indexIntArray# d# (i# -# a#)) s# of { s# ->
f# r# (i# +# 1#) s# }
| otherwise =
case writeIntArray# r# i# 0# s# of { s# ->
f# r# (i# +# 1#) s# }
(I# intlen#) = intlenUnboxed typeを用いた左シフトの実装例を示す。簡単のためにワード単位でのシフトしか実装していない。このルーチンは、Integerとシフト数を受け取って、元の配列のサイズ+シフト数のサイズの新しい配列を確保し、元の配列の内容を上位にコピーし、下位は0で埋めるということを行っている。newByteArray#はサイズと状態を受け取り確保した配列と状態を返す。状態というものをここで説明するのは避けるが、ここではCleanにおけるWorldのようなものだと思って差し支えない。操作に順序をつけるのに用いる。上のコードではどこからとも無くrealWorld#というもの(realWorld#が何物なのかはよく知らない)を持ってきてnewByteArray#に渡しているが、これを外部から供給させることも出来て、その場合自分でIOモナドなどを作ることも可能である。caseでバインドしているのは次のStateを同じ変数名でバインドしたいため。letはスコープが再帰的なため全て違う名前にしなければならなくて煩雑になる。newByteArray#はMutableByteArrayというものをを返す。Integerが保持するのはMutableではないByteArray#なので、unsafeFreezeByteArrayを用いてByteArrayに変換する(実際には多分何もやらない)。元のMutableByteArray#への参照が残っていればそれ経由で本来書き換えられないはずのByteArrayが書き換えられてしまうのでunsafeなのであるが、この場合Freeze以降には絶対に書き換えられることはないので、これで問題ない。
基数変換
SPOJに入っているGHCは6.4.1なので、基数変換も自前でやってやらないとSPOJに通すことは出来ない。
(コードは後に添付するプログラムを参照)
ある数を n = a*10^n + b (b<10^n) のように分解すると、10進表記においてはbからaに繰り上がりが発生しない。このような分解は単純にnを10^nで割って商をa、余りをbとすればよい。aとbが同じようなサイズになるように分解すればあとはa,bそれぞれに対してある程度小さくなるまで再帰的にこれを適用すればよい。割り算を高速に行えばこのアルゴリズムはO(nlogn)になる。後に示すコードでは与えられた数をn*(2^-x)の16進小数であるとみなしてそれに対応するために多少のコード変更がなされている。その際にここでもUnboxed typeを操作している。
プログラム
実際に作成したプログラムを示す。
http://fxp.infoseek.ne.jp/tmp/spoj/SQRT2.hs
ニュートン法で√2を求める
http://fxp.infoseek.ne.jp/tmp/spoj/E.hs
テーラー展開によりeを求める
http://fxp.infoseek.ne.jp/tmp/spoj/PI.hs
AGMを用いてPIを計算
http://fxp.infoseek.ne.jp/tmp/spoj/PI2.hs
arctan公式を用いてPIを計算
http://fxp.infoseek.ne.jp/tmp/spoj/PI3.hs
チュドノフスキーの公式を用いてPIを計算
PIの計算速度は3,2,1の順に速い。
