メニエスの名のもとに

プログラミング関係を中心としたぐだぐだブログ

MTToolBox バージョンアップ＆インターフェイス変更

MTToolBox

人の迷惑顧みず、やってきましたバージョンアップ。というわけで、MTToolBoxの以前から（私の頭の中で）気になっていたところを修正しました。その結果として、以前のバージョンとの互換性がなくなりました。まあ、あれですよ。どっかでやらねばならない修正ということです。つまりは、最初の設計が悪かった。ごめんなさい。
まあ、そうは言っても、ユーザーなんて片手くらいしかいないでしょ。
互換性がなくなったので、思い切ってバージョンを1.0にしました。なんとなく1.0でゴールした気になって、更新が止まりそうな気もします。
github.com

Ziggurat アルゴリズム

C++ 疑似乱数

この前から、このブログでは、なんとなく特定の分布、特に正規分布の疑似乱数を生成する方法について書いているわけだが、
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*とかの小さい疑似乱数生成器と組み合わせてもよいと思う。

dSFMTを元にした一様分布と正規分布をちょっと修正した

uniformIntFromDoubleクラスのテンプレートから返却値の型タイプを削除して、返却値はint32_t固定にした。いや、だってint64_tとか指定されても対応してないのですよ。それに、int8_t とか指定されても特別速い処理が出来るわけでもないし。
あとは、outputとspeedにカウントがなどが指定できるようした。
そして、sfmt-dist-0.1.tar.gzファイルを追加した。これがないと autotoolsから実行しないとならなくなるかも。
いや、これでだいたいいいや。
あとは要望があって、対応する暇と気力がある時まで凍結。
manieth.hatenablog.com

dSFMTを元にした一様分布と正規分布

なんかね、コンパイルする度に速さが変わるんで、訳がわからない。
とりあえず、他の環境でもコンパイル出来るようにしたので、興味がある人はどれが速いのか試してみて欲しい。
github.com

一般的なautotoolsなので

./configure
make
make check

で作成とチェックをして、実行ファイルはインストールされないので、

cd src
./speed
speed [-n|-u] [-m|-M|-S|-d] [seed]
-n     normal distribution
-u     uniform distribution
-m     mersennetwister19937
-M     MersenneTwister19937
-S     SFMT19937
-d     dSFMT19937
seed   seed

と使用法が表示されるので

./speed -u -d

などとやると100000000回乱数を生成してその時間を表示する。しまった、回数も指定できるようにするべきだった。で、速度が早くても分布はどうなのよって話になるのでその時は output を使うと同じような指定で乱数を1000個出力するので、出力ファイルを R かなんかで検定して納得できたら使ってみるというのはどうでしょうか。

で、speed でうちの環境で測ってみると、前回とは大違い。なんでこう毎回違う結果になるのか
uniform_int の出力

std::mt19937 dSFMT
970ms 167ms

normal_distributionの出力

std::mt19937 dSFMT
2258ms 1044ms

前回は
manieth.hatenablog.com

もし、dSFMTなり、SFMTなりが速くて使いたいなら

make install

で、ヘッダファイルとライブラリがインストールされる。デフォルトだと/usr/localに入ってしまうので、./configureするときに--prefixでインストールディレクトリを指定した方がいいかも。

-O2 と -O3 は大違い

C++ プログラム

Oだけに。
いやあ、これまでg++ でコンパイルする時に -O3 を付けてたんだけど、-O2 とはそんなに違わないと思ってたんですよ。で、GitHubにソースを置くようになってからautotoolsを使うようになって、Makefile.am でやはり -O3 と指定して訳です。でも、ちょっとしたテストをするときのために、Makefile.private というのを作って、automakeの枠外で好き勝手にmakeしたりしていたのですが、今回いろいろやっていたら実行速度が大幅に違うのである。なぜ実行速度が違うのかわけが分からなくてしばらく悩んでいたのですが、autotools でmakeのメッセージをよく見ると、-O3 と指定したあとに、-O2 の指定が出てくるのである。ぐぐってみると、statckoverflow サイトに回答があった。
configure.acに

AC_INIT
...
: ${CFLAGS=""}
: ${CXXFLAGS=""}
AC_PROG_CXX
AC_PROG_CC
...

と書かなければならないらしい。
いや-O2と-O3で何が違うのかよくわからないが、どうもC++のテンプレートの最適化が違うようである。randomのuniform_int_distributionが大幅に違うのだ。そうそうg++のバージョンは5.4.0である。
uniform_int_distributionは優秀でした。
std::mt19937での10^8個の乱数生成

32bit int uniform int
318ms 381ms

sfmt19937をstd::randomのインターフェイスに合わせて10^8個の乱数生成

32bit int uniform int
114ms 168ms

ほんとにこんなに速いのか、計測を間違っているんじゃないかと思うくらい。

O2に戻してみると
std::mt19937での10^8個の乱数生成

32bit int uniform int
435ms 469ms

sfmt19937をstd::randomのインターフェイスに合わせて10^8個の乱数生成

32bit int uniform int
112ms 165ms

あれぇ、そんなに変わらないじゃないか。おかしいなぁ。

メルセンヌ・ツイスターの小ネタ：自動ベクトル化

C++ 疑似乱数

なんとなく時間が半端に余ったので、１年位前に気づいたような気がするが、今はどうなっているか知らない小ネタをひとつ。

MersenneTwister をやや今風（あるいは俺流）に書くとこんな感じなんだけど

#define MATA UINT32_C(0x9908b0df)

#if defined(NOT_USE_MAG_ARRAY)
    const uint32_t mata32 = MATA;
#else
    const uint32_t mag[2] = {0, MATA};
#endif

   static inline uint32_t recursion32(uint32_t x, uint32_t y, uint32_t z) {
        x = (x & UPPER_MASK) | (y & LOWER_MASK);
#if defined(NOT_USE_MAG_ARRAY)
        return z ^ (x >> 1)
            ^ static_cast<uint32_t>((-static_cast<int32_t>(x & 1)) & mata32);
#else
        return z ^ (x >> 1) ^ mag[x & 1];
#endif
    }

    void fillState32(uint32_t * state)
    {
        int i = 0;

        for (; i + pos < size; i++) {
            state[i] = recursion32(state[i], state[i + 1], state[i + pos]);
        }
        for (; i < size - 1; i++) {
            state[i] = recursion32(state[i], state[i + 1],
                                   state[i + pos - size]);
        }
        state[size - 1] = recursion32(state[size - 1], state[0],
                                      state[pos - 1]);
    }

ヘッダーインクルードとか出力関数とかいろいろ省略してある。
これでNOT_USE_MAG_ARRAY を定義すると、配列参照ではなくて、ちょっとトリッキーなことをしてパラメータを入れるようになる。やってることは同じで、最下位ビットが1ならパラメータをXORする。

NOT_USE_MAG_ARRAYを定義して、clang++ で -O2 でコンパイルすると自動ベクトル化される。g++ では自動ベクトル化されなかった。-03だったかも。バージョンとかいろいろ適当。

一様乱数, それとは別に単精度浮動小数点数乱数

疑似乱数

コメントに答えるには、コメント欄が狭いかなと思ったが、書いてみたらそれほど長くもなかったような。

一様乱数

C++ のstd::uniform_int_distributionが遅いという話。ある範囲の数を一様に出力する乱数が欲しい。これは簡単なようでいてそれほど簡単ではない。
メルセンヌ・ツイスタのページのQandAコーナーよく聞かれる質問

サイコロみたいに、[0..N-1]の一様乱数を生成したい
N があまり大きくなくて(だいたい < 2^16) かつ、あまり精度を必要としなければ、 [0,1)区間版(倍精度）を使用して、出力をN倍して切り捨てれば十分です。 (注意:必ず半開区間バージョン[0,1) を使うこと。閉区間バージョンでこれをやると、約４２億回に１回の割合で 1.0(16進整数で言うと0xFFFFFFFF) が出現し、値として N が返されるというたちの悪いバグが起こります。 和田維作さんの指摘。2002/Jan.8th added)
WARNING: この方法は、Nが２の累乗でない限り、分布に誤差が生じます。 誤差のオーダーはN×2^-32になります。 より優れた、より安全な方法は、 整数乱数を生成し、N<=2nとなる最小のnをとって、二進で上位n桁 を見て、N以上の乱数は捨てることです。

関係ないけど、上付き文字を指定しても下付き文字に見えるんだよな。このページで下付き文字に見えるとしても上付き文字だと思ってください。
C++ のstd::uniform_int_distributionはここでいうところの「より優れた、より安全な方法」を使っていると思われる。
この生成法ならば、（可能な限り）正確な一様乱数を生成することができるが、ifによる条件判断が入る。それどころか、運が悪ければ、１個生成するのに、複数回条件判断で乱数を捨てなければならない。
しかし、サイコロ(N=6)のように疑似乱数生成器の出力ビット数に対して、Nが十分に小さければ、誤差は小さいので最初の方法で生成しても良いかもしれない。そこはユーザーの判断となるだろう。つまり、Nが小さく速度が欲しければstd::uniform_int_distributionは使わないという選択もあるということである。
また上記のQandAで「誤差のオーダーはN×2^-32」と言っているのは、32ビットMTから倍精度に変換したものから作っているから、53ビット精度のdouble を使用すれば、誤差のオーダーはN×2^-53になる。
分布の誤差を小さくする方法としては、64ビット整数から80ビット拡張倍精度浮動小数点数にして、簡易法を使えばN×2^-64の誤差に出来る。あるいは、４倍精度を使う。
どれが速いかは、例によって全部やって比較しないと分からない。

単精度、半精度浮動小数点数

32bit整数から作ればよいのだが、dSFMTのように浮動小数点数形式で作ればより速いと考えられる。しかし、ここで一つ問題がある。それは擬似乱数生成器の統計的検定をするために、TestU01というパッケージを使っていることである。いや、このパッケージの中のBigCrushというテストが非常に多くの統計的検定をやってくれるので大変便利なのである。だが、BigCrushは基本的に32ビット用である。64bitの疑似乱数を検定する時は、上位32ビットだけの出力、下位32ビットだけの出力、上位下位交互の出力というように分けて検定している。
しかし、単精度をダイレクトに作ると、23ビットの生成、半精度だと10ビットの生成である。これをそのままBigCrushに掛けると、まず検定を通らない。BigCrushを使わずに個々の統計的検定関数を呼べばよいのだが、パラメータが多くて相当に大変（だと予想される）。さらに、適当なパラメータを設定して検定を通ったとしても、そのパラメータが正しいのか、意味のある検定をしているのかという疑問が出てくる。まあ、論文を投稿すれば査読者に指摘されるだろう。
（実をいえば、BigCrushで使っているパラメータは本当に適切なのかという疑問もある。だが、エライ教授が監修して作ったパッケージでみんなが使っているので、信頼されている）
32ビット整数から作った場合に比べて本当に速いのかというのも、作ってみないと分からない。
この辺は、学部生なり院生なりが卒論、修論として取り組むと面白いかも知れない。