Ziggurat アルゴリズム

この前から、このブログでは、なんとなく特定の分布、特に正規分布の疑似乱数を生成する方法について書いているわけだが、
manieth.hatenablog.com
manieth.hatenablog.com

Ziggurat アルゴリズム

正規分布作成にZigguratアルゴリズムというのがあるらしい。
Ziggurat algorithm - Wikipedia, the free encyclopedia
で、こいつが速い。だが

However, since the ziggurat algorithm is more complex to implement it is best used when large quantities of random numbers are required.

って書いてあるんだけど、実装なんて一回してしまえばあとは使うだけなんで、複雑もクソもあるもんかって思うわけですよ。ところが、zigguratアルゴリズムのプログラムがなかなか見つからなかった。そしてようやく見つけたのが
http://people.sc.fsu.edu/~jburkardt/cpp_src/ziggurat/ziggurat.html.
でも、これ疑似乱数生成器KISSと一体になってるんですよ。で、こちらとしてはメルセンヌ・ツイスタ推しなんで、アルゴリズムを分離したい。
というわけで、Zigguratアルゴリズムを分離してC＋＋テンプレートライブラリにしました。これでSFMTやMTと組み合わせて、正規分布の疑似乱数を生成できる。でも、元のプログラムがそうだったので、単精度浮動小数点の正規分布乱数しか生成できません。倍精度に対応するには、ルックアップテーブルを変更するべきなのだろうと思うが、それにはちゃんと理解してやらないといけないので、とりあえず単精度版だけ。
github.com
このZigguratブランチにある。まだ、masterブランチではないので注意（mergeしたのでmasterブランチになった。2016-10-5）。単精度でいいのかという疑問が消えないのだが、SFMTと組み合わせると元のコードのKISS版よりも速いし、SFMTがKISSよりも乱数としての質で劣るとも思えないので、まあ元のコードを使っていた人は置き換えてもいいのかもと思う次第である。SFMTではメモリが大きすぎるという場合は、gitの中には入ってないけど、TinyMTとかXORSHIFT*とかの小さい疑似乱数生成器と組み合わせてもよいと思う。