きっかけ

http://d.hatena.ne.jp/itchyny/20120304 この記事を見て、私も昔円周率を計算したことがあるのを思い出しました。 http://d.hatena.ne.jp/tanakh/20070506#p1 公式は同じくチュドノフスキーを使用、GHCのIntegerをそのまま用いて円周率を計算、 当時（5年前）のマシン（おそらくCore2あたり）で1億桁が15分ほどだったようです。

SPOJというサイトのコンテスト向けに、 $\sqrt{2}$ 200万桁計算の高速化 をCで書くとこから始めて、 ナイーブ乗算、カラツバ法、そしてFFT乗算の実装にいたり、円周率へ。 最初AGMを実装し、 古くて新しい方法、Arctan系の公式の実装を経て Binary Splitting Methodを知り、 チュドノフスキーの公式へ至りました。 さらにHaskellでの実装、当時起こった色々な問題への対処、 いろいろ思い出されます。

ところで、これを書いてる時に見つけた http://円周率.jp/ が内容充実しててすごいです。 円周率を計算したい方にはとても参考になると思います。

そんなわけで

円周率熱が再燃してきたので、再びちょこっと書いてみました。

https://github.com/tanakh/pi

Haskellで書いてあります。Haskellには円周率のような計算を行うにおいて、優れた点がいくつかあります。

  • 多倍長整数の組み込みサポート
  • 多倍長整数の演算がgmpによって実装されており、高速
  • 効率のよいコンパイラの存在
  • 並列化が容易

https://github.com/tanakh/pi/blob/master/PI.hs が実際に円周率を計算しているコードです。 とてもシンプルです。シンプルに書いて、それなりに速いです。 Corei5-2400のマシンで100万桁2.2秒、1億桁で629秒程度でした。 並列化が容易と書きましたが、並列化はしていません。 というのも、並列化しても（おそらく）メモリバンド幅がボトルネックで あまり速くならなかったからです。 (おそらく並列化で高速化を行うには、乗算の際に行うFFT（あるいはNTT?）の 並列化が必要なんだと思います。)

ここへ至る道程

そのようなわけでHaskellで円周率を求めるプログラムができたわけですが、 ここへ至る道程は実はそれほど容易なものではなかったのです。 今では当たり前に、期待したとおりに、期待した通りのパフォーマンスで このプログラムは動いてくれますが、 5年前、当時のGHCはたしかバージョン6.4だったと思いますが、 これはいくつかの問題を抱えていました。 そしてまともに計算を行うには、それらを解消する必要がありました。 それらの問題が1つずつ、処理系の方から改善されて、 今のGHCでは安心して計算が行えるようになったのです。 だいたいのことは冒頭に挙げた私のブログに書かれていることですが、 それらが今のGHCでは解消されてる、やったね！ということを少しご紹介したいと思います。

文字列への変換

GHCのInteger（多倍長整数型）は、内部的には2進表現で保持されています。 LSBから順に、延々と必要な分だけビット列が並んでいるというものです。 これは計算する上で速度的に有利ですが、 最終的に求めた値を10進で書きださなければなりませんので、 ここでその処理が発生することになります。 ところで、基数変換（ここでは2進→10進）は実は自明な処理ではありません。 Haskellでは show とやるだけで文字列になってしまいますが、 ことIntegerに関して言えば裏ですごい基数変換アルゴリズムが走ることになるのです。

基数変換とはどうやればいいのでしょうか？ 一番自明なのは、10で割り続けるというものでしょうか。 大きな数を10で割る。割った余りが一番下の桁です。 もう一回10で割る。あまりは2桁目です。 これを繰り返していくと、お望みの10進表現が得られるというわけです。 しかしこれは桁数に対して O(n2) の計算量が必要になります。 100万桁に対してはとてもじゃありませんが現実的な時間では動かないのです。

ところでGHC6.4の show の実装はまさにそのようになっていました。 値の計算自体は速く終わるのに、と悔しい思いをしながらも、 出力のためにKaratsuba基数変換を自前で実装して 文字列化しなければなりませんでした。

さて、現在のGHCのIntegerですが、 さすがにまともな実装に変わっていて、 割と高速に文字列化できるようになっています。 https://github.com/ghc/packages-base/blob/master/GHC/Show.lhs#L479

シフト

多倍長の浮動小数点数を扱う場合、 典型的なのは元の値を必要なだけシフトして大きな整数として扱う方法です。 例えば、小数点以下100桁計算したいのであれば、10100 を掛けて切り捨てた整数として保持します。 しかしこのようにすると掛け算の後に桁がずれるので、 ずれた分だけ、この場合は 10100 で割ってやる必要があります。

ところで、この割り算という操作がとてもコストが高いのです。 おそらくリーズナブルな実装は、ニュートン法などを用いて除数の逆数を求めて、 それを掛けるというものになるかと思いますが、これは掛け算とは 比べ物にならないほど時間がかかります。 なので、割り算は避けたいのです。 幸いなことにここでの割り算はシフトした数に限られます。 ならばこの数を割りやすい数、例えば、2x にすれば、 割り算は多倍長整数を右にシフトするだけになります。

そのかわり、求まる値が元の値に 2x を掛けたものになります。 これをそのまま出力してもよくわからないので、10進表記でもっともらしい数に変換する必要があります。 これは簡単で、求まった値に 10d （dは有効桁数）を掛けて 2x で割ればよいだけです。 もちろんここの割り算も右シフトで行うことができます。

ここで問題になるのが、シフトです。 GHC6.4での右シフトは、ああ、なんという事でしょう。 なんとも愚直に、2x で割る実装になっていたのです。 これでは何のためにシフトを使うようにしたのかわかりません！

そこで、当時高速にこれを処理するには、 Integerの内部表現に手をツッコミ、 #だらけのMagicHashを駆使して、 ごりごりメモリをいじくり、 無理やり整数をシフトする必要がありました。

さて、現在のGHCですが、これも改良されて、 gmpのシフト関数を呼ぶ実装になっています。 （https://github.com/ghc/packages-base/blob/master/Data/Bits.hs#L349） この改良は割りと最近だったでしょうか。 ともあれ、現在の我々は晴れて安心してIntegerをシフトできるようになったわけです。

そういうわけで

あれから5年たった今、GHC7.4.1では、特に変なHackを行うことなく、 Haskellで1億桁の計算ができるようになっていました！ やりました。 ようやく安心しておすすめできます。 みなさんもHaskellで多倍長計算、チャレンジしてみてください。