SFMTの256,512bit化への道その２
浮動小数点数版については前回の記事を参照してください。
整数版
整数版はすでにAVX2版が存在するけど、AVX512F版は存在しない。なので作ってもいいかもしれない。
実は、dSFMTで使用されているようなlungという変数を使った短いループを作る方法の方がゼロ過剰状態の回復などで優れている。なので浮動小数点版のSFMTではその方法を使っているのだが、整数版も同じようにしたらどうだろうという意見は、松本眞先生から当時出ていた。しかし、ちょっと試したところうまく行かなかったのであった。
これをもう一度試して、256bit, 512bit用の整数版SFMTを作るのも悪くはないだろう。
そんな感じで適当に考えたのがこの図である。
の部分は32bit単位の置換(permutation)を表している。
パラメータp1の値によって、32bit単位のローテーションの回数を指定する。
512bit版を作るならこんな風になるだろう
w256.hpp
ちょっとコードを書こうかと思ったが、途中ではてな記法に変更することはできなかったようだ。あとまわし。
SFMTの256,512bit化への道
結構遠い道のりかも知れない
疑似乱数生成には以下の三つが重要なのである。
- 周期が最大であること
- 乱数性が良いこと
- 速度が速いこと
面倒くさいね。
dSFMTの256, 512bit化
SFMTには整数を生成するいわゆるSFMTと浮動小数点数生成に特化したdSFMTというのがある。そして、整数用のSFMTは256bit版がある。
これはSFMTの配布にも取り入れられている。
この256ビット版SFMTは128ビット版と同じ乱数列を生成するものである。
それじゃあ、まあ、dSFMTの256bit版でも考えてみるか。ってな話から始まるのである。
で、SFMTの発想と同じで、アルゴリズムというか漸化式を変えてしまえば楽なんじゃないの、って考えるわけです。別に生成する乱数列が違ってもいいじゃないか。それはそれで周期と乱数性が保証されていれば。で、dSFMTの論文に乗っている図はこんなの
IEEE754倍精度浮動小数点数のビットパターンをF2線形（数学の人はアフィンという）演算で作るというのがミソなわけだ。これで作ると出来た乱数rは 1.0 <= r < 2.0 となるので、1.0を引いて出力するのである。
左側の箱が積んであるところが配列のつもりである。ここは１回生成する毎にインデックスを増やしてアクセスする内容をずらす。
のの部分は、64ビット毎の左シフトで、シフト量はパラメータになっている。instrinsics で言うなら_mm_slli_epi64。シフト量やマスク値、配列からの参照位置のパラメータセットは、最大周期になるように選ばれるので、勝手に変えていいものではない。つまり、128ビットのdSFMTのものを256ビット用に使うことは出来ない。
の部分はF2ベクトルとしての加算。またはビットごとの排他的論理和。instrinsicsでいうならずばり、_mm_xor_si128である。１回で3個足す命令はないけど。
の部分は、64ビット毎の右シフトで、シフト量は12ビット固定になっている。上位12ビットをゼロにすると同時に、上位ビットを下位ビットに(次の足し算で)混ぜ込むというスグレ技である。instrinsics で言うなら_mm_srli_epi64。
の部分は、128ビットのマスクをかける部分である。ここもマスクの値はパラメータになっている。instrinsics は_mm_and_si128である。
の部分は32bit単位の置換(permutation)を表している。
128ビットレジスタの内容を32ビット単位でぐるっと回しているだけ。rotateともいう。
instrinsicsは_mm_shuffle_epi32であり、このshuffle指定に使う値は 0x1b である。
これが32ビット単位でぐるっと回す値になる。
この部分が結構重要なんだけど、乱数性を良くするのとゼロ過剰状態からの回復を早くする上で役立っているはずなのである。ここはたとえば128ビットレジスタを使って全ビット使う。（IEEE754形式ではない）
赤い矢印が１回の生成の最後のステップを表す。SIMD生成なので1回で128ビット分生成する。浮動小数点数として2個を1度に作るわけだ。
256ビット版dSFMTの（想像）図
256ビット版にするなら上の図をこんな風に変えればいいんじゃないの。
ほとんど変わってないように見えるけど、全部の場所で256ビットを扱うように変わるってこと。
この部分は128bit版と同じ32ビットのpermuteなのだが、
instrinsicsはshuffleではなくて、_mm256_permutevar8x32_epi32となる。
ずっとshuffleという名前の命令を探していて見つからなかったが、そのまんまのpermuteだった。
512ビット版dSFMTの（想像）図
512bit版を作るならこんな風になるだろう。
この図を見ると、1度に2個作っていたのが1度に8個作るようになるんだから、多少のオーバーヘッドがあっても絶対に速いと思うかもしれない。しかし、それが素人の浅はかさってやつである。1度に8個作っても、使う時は1個ずつ使うのである。8個のうち何個目まで使ったかを調べるためにif文を1個入れるだけで、すべての高速化を無効にするだけの破壊力があるのだ。で、まとめて376個作るとかやるのだが、376個作るのなんてループ回数にすると47回ですよ。47回のループではループオーバーヘッドが大きいじゃないですか。つまり、SIMDのサイズが大きくなるとループのオーバーヘッドの割合が大きくなるのであう。かと言って一度に作る量を多くするとキャッシュに入らなくなる。まあ、乱数だけはキャッシュに入ったとしても、実際に使う時はユーザーのデータもキャッシュを使うわけで、ベンチマークだけ速いが、実際は遅いということにもなりかねない。なんてことは、実際に速度を測るところまで進まないとわからないのである。
それはともかく、このままじゃ絵に描いた餅なんだけどね。図を描いただけで動作するなら苦労はしないわけですよ。これでうまくいくのかどうか、最大周期を満たすパラメータを求めないとならないし、そもそもパラメータが一つも求まらない、つまりこの図の方式ではダメだということもある。そして、パラメータが求まったら乱数性が十分高いことを確認しなければならない。リリースするためには、最低周期を保証するための周期保証ベクターを計算して求めなければならない。
まあ、そういったこと全部をやってくれる（はずの）ライブラリMTToolBoxって便利なものがある。私が作ったものだけど。拙者の作ったMTToolBoxとてよいものでござるよ。ゲヒヒヒヒ。
更に、速度を計測しなければならないが、それにはちゃんとinstrinsicsを使って高速に実装する必要がある。これが結構大変というか、俺はすぐにミスをするので、速かったと思ったら出力が違っていたりする。そして、ここで遅かったら、最初から考え直さなければならない。なんと心折れる作業であることか。
まあ、学者だったらさらにその後、論文にしたりしなければならない訳だが、その過程が一番苦しいのかも知れない。今回はやらない。というか今後はやらない。というか、既にXSaddからやってない。
流行遅れ
どうも最近の流行は、crush free とかなんとか言うらしいですよ。よく知らんけど。TestU01という疑似乱数のテストをするライブラリがあって、その中にあらかじめ組み込まれた一連のテストがあるんですが、BigCrushテストというのが様々な統計的検定なんかを疑似乱数生成器に対して実行してくれる。このBigCrushテストを全部パスするというのが最近の疑似乱数開発のトレンドらしい。
しかし、上の図の方式で作られたdSFMTはBigCrushを通らないんですよね。必ず落ちるテストがある。だいたい Linear Complexity と Matrix Rank というテストに落ちる。それはdSFMTがF2線形疑似欄数生成器だからで、出力関数を非F2線形にすればパスするのです。
でも、でも、dSFMTはIEEE754倍精度浮動小数点数の内部ビット形式を直接生成するというトリッキーなことをしている関係上、出力関数を非F2線形にすることが簡単には出来ないのです。これを自業自得という。
実のところ統計的検定というのは、実際の疑似乱数の出力の一部だけ見ているので、初期化によって通ったり通らなかったりすることもある不完全なものである。統計的検定を通ったからといって、まったく問題ないと言い切ることは出来ないのだ。それに対して、メルセンヌ・ツイスタやdSFMTで計算している均等分布次元は理論的に求められた均一性なので、堅実な保証である。と言えないないこともない。
まあ、流行なんてものはそのうち一周回って戻ってくることもあるからね。
アンドロイドアプリ「ナンプレ練習」がとつぜん大量にダウンロードされて、広告がクリックされ、大金が転がり込んでくればなぁ。
朝から飲んだくれてアニメ見てられるのに。
ナンプレ：自動生成なのにうっかり綺麗な形になってしまった。
このブログでは、ナンプレ問題の自動生成とかをやっていた訳です。で、自動生成だからパズル雑誌に掲載されるような図柄の美しいパズルではなくて、難しいパズルを作ろうとしていたのでした。そのために、9x9の通常版ナンプレでは対称性をちょっと捨てて１個だけなら対称でなくてもいいという条件でパズルを生成していました。
そうは言っても、必ず対称でなくするというようにしていた訳ではないので、時には対称なパズルが生成されてしまうこともあるのです。そして、美しい形でなくてもいいやと思っていても偶然、そこそこ美しい形になることもあるのです。
そういう判定をしていないので、たまたまそうなってしまったのです。このパズルの生成された形が美しいかどうかの判定を機械学習ででもやらせれば、もう完全にパズル雑誌に掲載されるレベルのパズルが自動生成できると思います。が、美しさにそんなにこだわりがないので私はやる気がない。
もし、そういうことをやりたい人がいれば、GitHub - MSaito/NumPl: Number Place Puzzle Generatorを使ってやってみてください。（おっと、push してなかった。git push っと。いや、流石に自分のアプリ（のXYZ-Wing対応版）出す前に公開するのは人が良すぎると思って公開していなかったのだ。まあ、今日アプリを公開したから）
というわけで、デバッグ中にみつけた対称形のパズルである。
ナンプレ練習 広告付き無料版 - Google Play の Android アプリの XYZ-Wing #0003
XYZ-Wingの解き方：
Ⅳ－５ ＸＹ-Ｗｉｎｇ - 南碁空 ナンプレ（数独）攻略の広場
と英語だけど
dSFMTdc サンプル
dSMFT のパラメータを生成するプログラム
MTToolBox のサンプルにdSFMTのパラメータ生成プログラムを追加しました。
いまのところ、私が作った疑似乱数生成器の中である意味一番出来がいいのがdSFMTである。まあ、IEEE754形式で生成というのも、lungを使うというのも松本眞先生の提案だが、高速化のポイントになっている初期化でIEEE754フォーマットで使う定数をセットして、生成時には定数を参照しないというのは私の工夫なのである。
これで、過去に作った疑似乱数生成器は全部MTToolBoxのサンプルに入りました。疑似乱数生成器を作成する上で、私の知っているノウハウは公開したということになり、肩の荷が降りた気がします。github.com
dSFMTの論文にミスがありました。
均等分布次元計算プログラムのバグです。
dSFMTの場合、均等分布次元は初期値によって変動する可能性があり、その変動の中でも最悪でもこれだけはあるという値を論文に掲載していました。しかし、正しい値はそれよりも悪い値になります。なお、論文の結論自体は変わりません。
ミスがあることはSFMTのミスの時に分かっていたんだけど、正しい値を計算するのに無駄に手間取ってしまった。９月は政治的状況にうんざりしてやけ酒飲んでいたせいもある。
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/articles.html#dsfmt
SFMTの論文にミスがありました。
均等分布次元計算プログラムのバグです。
SFMTの場合、均等分布次元は初期値によって変動する可能性があり、その変動の中でも最悪でもこれだけはあるという値を論文に掲載していました。しかし、正しい値はそれよりも悪い値になります。なお、論文の結論自体は変わりません。
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/articles.html#sfmt
sfmtdc サンプル
技術的負債
MSaito/MTToolBox · GitHub
MTToolBoxにSFMTDCのサンプルを追加した。
昔の自分のプログラムを直すのが、こんなに大変だとは思わなかった。
でも、俺にしか出来ないという程ではなくても、他の人がやったら無駄に時間を使うと思うので、優秀な人材の時間を浪費しないためにも俺がやるしかないのであった。なんとか動くが修正部分のドキュメントは対応されていない。いつかドキュメントも修正できたらいいな。
なお、ついでなので、SFMTのシフト量を固定にするsfmtdcfixedも作成しました。これで、コンパイルしなおさなくても複数のパラメータに対応することが出来るSFMTが作れるはずです。ただし、固定されたシフト量はいい加減なので、修正して使ってください。シフト量を決める方法としては、sfmtdcでシフト量を固定しないでパラメータを生成し、その中で均等分布次元のよいもののシフト量を選んでそれに固定してしまうのがよいと思います。mexpごとに適切なシフト量にした方が良いのかもしれないし、mexpには依存しないのかも知れないが、いまちょっとその実験をしている余裕はありません。誰かやって、卒論とか修論にするといいと思います。
このサンプルを元にして頑張れば、256ビットや512ビットのSIMDに対応した疑似乱数生成器を作ることも出来るはずですが、保証は出来ない。どこかで128ビットにしか対応していない部分があるかも知れない。
なお、SFMTのみで、dSFMTのdynamic creatorのサンプルはまだ作れていないのです。
世界に混沌を
昔、松本眞先生のゼミのときに、MTToolBoxという疑似乱数開発ツールを作って公開する。そうしたら、卒論とか修論で疑似乱数を作る人が続出して、新しい疑似乱数生成法が雨後の筍のようににょきにょきと現れて、どれを使ったらいいやら分からなくなり、（疑似乱数の）世界が混沌に満たされるようになる。そうしたいと言ったのですが、全然そうなる気配はありません。一方、カオス系を使った乱数などという混沌は後を絶たないようですが。
今季アニメ
もう８月も下旬ですよ。「Charlotte」と「うしおととら」「わかば＊がーる」「境界のRINNE」「ミス・モノクローム」くらいですかねぇ。バンダイチャンネル見放題の範囲内では。「ガッチャマンクラウズ」も見ているけど、ちょっと判断に苦しむ。
来週はＳＦ大会
だからどうってわけでもないけど。米子に行くのでその場で私に現金をくれるとかいう人はメールでもください。物凄く楽で物凄く儲かる仕事でもいいです。
