メニエスの名のもとに

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

SFMTの256,512bit化への道

結構遠い道のりかも知れない

疑似乱数生成には以下の三つが重要なのである。

  • 周期が最大であること
  • 乱数性が良いこと
  • 速度が速いこと

面倒くさいね。

dSFMTの256, 512bit化

SFMTには整数を生成するいわゆるSFMTと浮動小数点数生成に特化したdSFMTというのがある。そして、整数用のSFMTは256bit版がある。

github.com

これはSFMTの配布にも取り入れられている。

github.com

この256ビット版SFMTは128ビット版と同じ乱数列を生成するものである。


それじゃあ、まあ、dSFMTの256bit版でも考えてみるか。ってな話から始まるのである。

で、SFMTの発想と同じで、アルゴリズムというか漸化式を変えてしまえば楽なんじゃないの、って考えるわけです。別に生成する乱数列が違ってもいいじゃないか。それはそれで周期と乱数性が保証されていれば。で、dSFMTの論文に乗っている図はこんなの

f:id:manieth:20160630133723p:plain

IEEE754倍精度浮動小数点数のビットパターンをF2線形(数学の人はアフィンという)演算で作るというのがミソなわけだ。これで作ると出来た乱数rは 1.0 <= r < 2.0 となるので、1.0を引いて出力するのである。

左側の箱が積んであるところが配列のつもりである。ここは1回生成する毎にインデックスを増やしてアクセスする内容をずらす。

f:id:manieth:20160630133721p:plainのの部分は、64ビット毎の左シフトで、シフト量はパラメータになっている。instrinsics で言うなら_mm_slli_epi64。シフト量やマスク値、配列からの参照位置のパラメータセットは、最大周期になるように選ばれるので、勝手に変えていいものではない。つまり、128ビットのdSFMTのものを256ビット用に使うことは出来ない。

f:id:manieth:20160630133720p:plainの部分はF2ベクトルとしての加算。またはビットごとの排他的論理和。instrinsicsでいうならずばり、_mm_xor_si128である。1回で3個足す命令はないけど。

f:id:manieth:20160630133718p:plainの部分は、64ビット毎の右シフトで、シフト量は12ビット固定になっている。上位12ビットをゼロにすると同時に、上位ビットを下位ビットに(次の足し算で)混ぜ込むというスグレ技である。instrinsics で言うなら_mm_srli_epi64。

f:id:manieth:20160630133717p:plainの部分は、128ビットのマスクをかける部分である。ここもマスクの値はパラメータになっている。instrinsics は_mm_and_si128である。

f:id:manieth:20160630133719p:plainの部分は32bit単位の置換(permutation)を表している。
128ビットレジスタの内容を32ビット単位でぐるっと回しているだけ。rotateともいう。
instrinsicsは_mm_shuffle_epi32であり、このshuffle指定に使う値は 0x1b である。
これが32ビット単位でぐるっと回す値になる。

f:id:manieth:20160630133716p:plainこの部分が結構重要なんだけど、乱数性を良くするのとゼロ過剰状態からの回復を早くする上で役立っているはずなのである。ここはたとえば128ビットレジスタを使って全ビット使う。(IEEE754形式ではない)

赤い矢印が1回の生成の最後のステップを表す。SIMD生成なので1回で128ビット分生成する。浮動小数点数として2個を1度に作るわけだ。

256ビット版dSFMTの(想像)図

256ビット版にするなら上の図をこんな風に変えればいいんじゃないの。

f:id:manieth:20160630133722p:plain

ほとんど変わってないように見えるけど、全部の場所で256ビットを扱うように変わるってこと。

f:id:manieth:20160630133715p:plainこの部分は128bit版と同じ32ビットのpermuteなのだが、
instrinsicsはshuffleではなくて、_mm256_permutevar8x32_epi32となる。
ずっとshuffleという名前の命令を探していて見つからなかったが、そのまんまのpermuteだった。

512ビット版dSFMTの(想像)図

512bit版を作るならこんな風になるだろう。

f:id:manieth:20160630133724p:plain

この図を見ると、1度に2個作っていたのが1度に8個作るようになるんだから、多少のオーバーヘッドがあっても絶対に速いと思うかもしれない。しかし、それが素人の浅はかさってやつである。1度に8個作っても、使う時は1個ずつ使うのである。8個のうち何個目まで使ったかを調べるためにif文を1個入れるだけで、すべての高速化を無効にするだけの破壊力があるのだ。で、まとめて376個作るとかやるのだが、376個作るのなんてループ回数にすると47回ですよ。47回のループではループオーバーヘッドが大きいじゃないですか。つまり、SIMDのサイズが大きくなるとループのオーバーヘッドの割合が大きくなるのであう。かと言って一度に作る量を多くするとキャッシュに入らなくなる。まあ、乱数だけはキャッシュに入ったとしても、実際に使う時はユーザーのデータもキャッシュを使うわけで、ベンチマークだけ速いが、実際は遅いということにもなりかねない。なんてことは、実際に速度を測るところまで進まないとわからないのである。

それはともかく、このままじゃ絵に描いた餅なんだけどね。図を描いただけで動作するなら苦労はしないわけですよ。これでうまくいくのかどうか、最大周期を満たすパラメータを求めないとならないし、そもそもパラメータが一つも求まらない、つまりこの図の方式ではダメだということもある。そして、パラメータが求まったら乱数性が十分高いことを確認しなければならない。リリースするためには、最低周期を保証するための周期保証ベクターを計算して求めなければならない。

まあ、そういったこと全部をやってくれる(はずの)ライブラリMTToolBoxって便利なものがある。私が作ったものだけど。拙者の作ったMTToolBoxとてよいものでござるよ。ゲヒヒヒヒ。

github.com

更に、速度を計測しなければならないが、それにはちゃんとinstrinsicsを使って高速に実装する必要がある。これが結構大変というか、俺はすぐにミスをするので、速かったと思ったら出力が違っていたりする。そして、ここで遅かったら、最初から考え直さなければならない。なんと心折れる作業であることか。

まあ、学者だったらさらにその後、論文にしたりしなければならない訳だが、その過程が一番苦しいのかも知れない。今回はやらない。というか今後はやらない。というか、既にXSaddからやってない。

流行遅れ

どうも最近の流行は、crush free とかなんとか言うらしいですよ。よく知らんけど。TestU01という疑似乱数のテストをするライブラリがあって、その中にあらかじめ組み込まれた一連のテストがあるんですが、BigCrushテストというのが様々な統計的検定なんかを疑似乱数生成器に対して実行してくれる。このBigCrushテストを全部パスするというのが最近の疑似乱数開発のトレンドらしい。

しかし、上の図の方式で作られたdSFMTはBigCrushを通らないんですよね。必ず落ちるテストがある。だいたい Linear Complexity と Matrix Rank というテストに落ちる。それはdSFMTがF2線形疑似欄数生成器だからで、出力関数を非F2線形にすればパスするのです。

でも、でも、dSFMTはIEEE754倍精度浮動小数点数の内部ビット形式を直接生成するというトリッキーなことをしている関係上、出力関数を非F2線形にすることが簡単には出来ないのです。これを自業自得という。

実のところ統計的検定というのは、実際の疑似乱数の出力の一部だけ見ているので、初期化によって通ったり通らなかったりすることもある不完全なものである。統計的検定を通ったからといって、まったく問題ないと言い切ることは出来ないのだ。それに対して、メルセンヌ・ツイスタやdSFMTで計算している均等分布次元は理論的に求められた均一性なので、堅実な保証である。と言えないないこともない。

まあ、流行なんてものはそのうち一周回って戻ってくることもあるからね。

アンドロイドアプリ「ナンプレ練習」がとつぜん大量にダウンロードされて、広告がクリックされ、大金が転がり込んでくればなぁ。

朝から飲んだくれてアニメ見てられるのに。