2007-05-06GHCのIntegerの深淵Haskell前書き先日の日記のSPOJにおける√2を200万桁求めるという課題で、私はCを用いてそれを達成したのであるが、Haskellで100万桁を求めている人がいて、それがずっと気になっていた。普通に計算すると10万桁で精一杯で、とてもHaskellで100万桁なんて無理のように見える。ところが、ふとしたことから、Haskellで100万桁が可能であるように思えてきて、実際にやってみたら出来てしまった。そのためにはGHCのIntegerに対する深い理解が必要なのであった。ついでの円周率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の上でなされており、それなりの速度で変換することが出来る(とはいえそんなに速くない)。shiftIntegerは16進数で保持されている。IntegerはData.Bitsのインスタンスである。*(2^x)、`div`(2^x)、によって実装されている。2の累乗による乗除算はO(n)で出来るはずなのに、GHCではO(nlogn)かかってしまう。typeHaskellで多倍長演算を行うには、上記のような問題点を解決しなければならない。そのためには、C言語でIntegerの中身をいじるルーチンを書くか、あるいはHaskellからIntegerの中のデータをいじるかの必要がある。今回の場合、SPOJに通したかったので、Haskell内で完結させる必要がある。HaskellからIntegerのデータをいじる方法をとることにする。typeというものを使う。GHC.ExtsかGHC.Primあたりのモジュールをimportして-fglasgow-extsオプションを渡すと使うことが出来る。Int#などサフィックスに#が付いた型がそれにあたる。普通のHaskellの値は型情報の付いたオブジェクトになっており、整数などの本来1ワードであらわせるデータも例外でない。typeが示すのはそのような情報が付随しないデータであり、普通の型とは種も異なっている。プリミティブな型の値を操作するためにこのunboxedな型を用いることが出来る。typeはこのように定義されている。S#は1ワードに収まる範囲の数を表し、J#は多倍長整数を表している。J#の一つ目の引数は桁数と符号をあらわしており、二つ目のByteArray#型の引数が整数のデータにあたる。ByteArray#は単なるメモリブロックへのポインタのようなものである。これに対してはHaskellからデータの読み取りをすることが出来る。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に通すことは出来ない。ように分解すると、10進表記においてはbからaに繰り上がりが発生しない。このような分解は単純にnを10^nで割って商をa、余りをbとすればよい。aとbが同じようなサイズになるように分解すればあとはa,bそれぞれに対してある程度小さくなるまで再帰的にこれを適用すればよい。割り算を高速に行えばこのアルゴリズムはO(nlogn)になる。後に示すコードでは与えられた数をn*(2^-x)の16進小数であるとみなしてそれに対応するために多少のコード変更がなされている。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.hsAGMを用いてPIを計算http://fxp.infoseek.ne.jp/tmp/spoj/PI2.hsarctan公式を用いてPIを計算http://fxp.infoseek.ne.jp/tmp/spoj/PI3.hsチュドノフスキーの公式を用いてPIを計算PIの計算速度は3,2,1の順に速い。備考Windows版GHCで大きな桁数のPIを計算するとなぜか落ちる。(GHCのバグだろうか)64bitCPU上のLinuxにおいて円周率1億桁を15分ほどで計算できた。HaskellのIntegerの性能は(うまく使えば)かなり良い。HaskellHaskellを現場言語にするために、こんなものを作ってみました。writeup月別アーカイブはてなブログをはじめよう!tanakhさんは、はてなブログを使っています。ブログを報告するスターをつけました引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。