差分表示


*&date(Y-n-j[lL],2010/1/26); xorshift乱数のメモ

XorShift乱数は、[[こちら>http://ja.wikipedia.org/wiki/Xorshift]]とかのように短い関数で紹介されるけれど、たいていシード固定で書かれていて、シードを外部から与えたい場合はどうしたら、と悩んでしまう.
もちろん、ちょっとぐぐったらシード設定できるソース載せてるサイトあるし、
おおもとの[[paper>http://www.jstatsoft.org/v08/i14/paper]] に
all 0でなければよいようなことかいてあるようだけど、
どの程度適当でよいのか不安にもなる.
で、今さらながら[[2chの擬似乱数2>http://www.unkar.org/read/pc11.2ch.net/tech/1192628099]]
をみたら76に解説&seed設定付ソースがあった.
シードは極端な設定をすると不自然な部分列が出力されるからMTのを参考にしたとある.
ありがたく、これ流用させていただくことに...
と、よくよくみれば、[[ここ>http://meme.biology.tohoku.ac.jp/klabo-wiki/index.php?%A5%D7%A5%ED%A5%B0%A5%E9%A5%DF%A5%F3%A5%B0%B8%C0%B8%ECC%2B%2B]]の初期化とほぼ同じ(改良版)のよう.
//(試してないけど、この式みてたら線形合同法な乱数(たとえばs=s*1566083941+1とか)で初期化するのもありかなって気も)

その他、検索したときのサイトメモ.

[[mt-lite>http://mt-lite.sourceforge.net/doc/ja/index.html]]の[[実行速度>http://mt-lite.sourceforge.net/doc/ja/evaluation.speed.html]]に
xorshiftを含む各種乱数の速度比較あり. 

[[良い乱数・悪い乱数>http://www001.upp.so-net.ne.jp/isaku/rand.html]]に XorShiftと他の乱数との速度比較とかSSE2使ったバージョンとかがある. 同サイトの
[[こちら>http://www001.upp.so-net.ne.jp/isaku/rand2.html]]にもいろいろなコンパイラでの測定比較あり.

[[xorshift.pdf>http://www.iro.umontreal.ca/~lecuyer/myftp/papers/xorshift.pdf]]
[[上記サイト>http://www001.upp.so-net.ne.jp/isaku/rand.html]]で、"問題があると書かれている" とあったのでみた(翻訳ソフトでわからないなりに).
元のpaperより突っ込んだ性能テストをして、通過しないテストがいくつか(も)あるよう.
(XorShiftは3つのシフト数の有効な値の組み合わせがいくつか(も)あるけれど、
組み合わせによって性能にばらつきがあるよう)
あと改良版として256ビット版が載ってる.
(周期のビット数だけでなく、計算でのシフト数も増やして精度をあげている)


**** XorShift乱数の速度

[[擬似乱数2>http://www.unkar.org/read/pc11.2ch.net/tech/1192628099]]の 42 に
周期 32,64,96,128,160,192 ビット版のソースがあったので、
ちょっと速度比較してみた.

ソースは
[[コレ>http://www.6809.net/tenk/html/prog/xorshiftrand/XorShiftRand.h.html]]

[[コレ>http://www.6809.net/tenk/html/prog/xorshiftrand/MtRand.h.html]]

[[コレ>http://www.6809.net/tenk/html/prog/xorshiftrand/test.cpp.html]].
([[zip>http://www.6809.net/tenk/html/prog/xorshiftrand/XorShiftRand_test.zip]])

>96ビット版が要修正だったり64ビット版は出力も64ビットだったり だったので 結局 [[paper>http://www.jstatsoft.org/v08/i14/paper]] みて修正. あと、xorshift.pdfの256ビット版追加. (シード設定可のクラスにしたりでエンバク不安あり. 64,96,160版はA,B,Cの選択がこれでよいのか自信なし).&br();で paper の160のソース、>>Aになっているけど <<Aの誤記だと思われる. 他と比べてもだし、[[別のチェック>http://www.6809.net/tenk/html/prog/xorshiftrand/test2.cpp.html]]で結果がメタメタになった.

結果は以下. 環境は基本Phenom2 3GHz (玄箱HGは PowerPC 266MHz).
~
100000000回実行の平均. 単位:ナノ秒.
,                     ,  vc9&br();(32bit),  vc9&br();(64bit),mingw&br();g++3.4.5, mingw&br();g++ 4.4.0,    dmc8.51 ,     owc 1.8,   bcc 5.8.2,玄箱HG&br();g++4.1.2
,XorShift  32          ,       3.411,       2.662,       3.329,       2.997,      10.480,       3.355,       5.010,      23.100
,XorShift  64          ,       2.016,       2.167,       1.842,       1.996,       5.159,       3.000,       3.734,      34.700
,XorShift  96          ,       2.335,       4.345,       1.667,       2.041,       5.773,       3.023,       3.512,      38.500
,XorShift 128          ,       2.331,       3.963,       2.419,       1.995,       5.444,       3.230,       3.685,      46.200
,XorShift 160          ,       2.334,       4.206,       2.338,       2.026,       6.231,       3.002,       3.187,      46.200
,XorShift 192          ,       2.584,       3.874,       2.876,       2.336,       8.115,       3.157,       3.566,      46.200
,XorShift 256          ,       7.397,       6.264,       6.327,       6.491,      17.224,       6.184,       7.197,     134.800
,線形合同法の乱数      ,       1.509,       1.336,       1.664,       5.784,       3.581,       1.329,       3.154,      23.200
,clib rand             ,      18.559,      14.174,      15.080,      14.720,       6.534,       5.236,       3.231,     461.900
,MT Rand               ,       5.452,       5.679,       7.996,       5.869,      31.299,       9.289,       8.280,     127.200


下3つは比較用.
MT乱数([[mt19937ar-cok.c>http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c]])と
線形合同法を使った乱数とC標準ライブラリのrand()を使ったもの.

コンパイラごとに癖はあるけれど、大筋は似た傾向.

paperの例では1.8GHzで4~6ナノ秒とあったので、
3GHzの環境(vcやg++)で 2~4ナノ秒は、そんなものかもしれない.
~
なお、これら(の速い値)はインライン展開されている(はず).
関数呼び出しになっている場合は、(vcにて) 3~6ナノ秒台だった
(テストルーチンで乱数をファンクタにせずワンクッション関数を噛ましてたとき)

全体に周期32ビット版は他より遅くなっているのは、たぶん(アウトオブオーダーとかのレベルの)命令の並列化ができないためだろうか.

速度的には、
大筋ビット数が増えれば遅くなるけれど、32-192の範囲では
増える周期(質の改善)に比べ時間の増加はさほどでもないので、
どうせなら128でなく160や192あたりを使うのも手かもしれない.

256ビット版はプログラムが少々複雑になっているせいか
時間増加のよう. 128ビット版にくらべ2~3倍...
MTと対して変わらない速度(場合によっては負ける)なので出番はないだろう.

でinline展開で速いのはいいけれど、
ビット数が大きいほど(ささやかとはいえ)コードサイズが膨れやすくなるので、
乱数の質が高くなくていいなら、ビット数少ないのを選ぶのもよいかもしれない.~
でも線形合同法でいいなら、(inline展開されたら)そっちのほうが速いか.

サイズきになるならinlineしない形に書いて160,192あたりを...
もっとも、このソースの64,96,160(,192)は使ったパラメータがほんとにこれがいいのかわかっていないので
結局 128 を選んでおくのがベターな気はする.

~

あと、XorShiftの件とはずれるが、vcやg++でのcライブラリの rand()の速度が遅いのは、
マルチスレッド対応のオーバーヘッドのためだろうか.
インライン展開された線形合同法の乱数の10倍前後なのは
ペナルティが大きいとも、10数ナノ程度ならば大抵たいしたことないとも、どちらとも.
結局使用量/使用箇所しだい、だけれど、どっちにしろ精度のことを思えば、
とっとと XorShiftなりMTなり用意して用途別にインスタンスを管理したほうが、
楽になれると思う.

(参考サイトにある SIMD版のSFMTとかが使える環境なら, そっちのほうがいいだろうな)

~

**** 追記

乱数の精度の評価についてはさっぱりなんですが
(自分が使った処理において出目に不都合があるかどうかくらいなら、ともかく)、
[[同根のtest2>http://www.6809.net/tenk/html/prog/xorshiftrand/test2.cpp.html]]をしたときの結果はあまりよいという感じではなかった.
[[同根のtest2>http://www.6809.net/tenk/html/prog/xorshiftrand/test2.cpp.html]]をしたときの結果はあまりよいという感じではなかった.(この場合のテストとして適切かどうかもわかってないけど)

まあ、たいていの場合の乱数の利用頻度や精度を思うと
普通は MT を選んどくのが無難だと思う.

この程度の速度差(時間差)やプログラムサイズ差が
ネックになるような利用頻度(仕様)やターゲット環境なんて
ゲームでもまずないだろうし.~
どちらかというとライブラリ化でのオーバーヘッドのほうが
気になる可能性があるかも.
(この程度のサイズのルーチンなんだから、わざわざメーカー提供のライブラリ使わなくてもいいよね、とか)


----
#comment