検証方法
検証方法としては、単純に、内部状態の各 bit に対し、その bit だけが立った乱数生成器を用意し、
その乱数生成器を状態遷移させていき、その内部状態の変化を観察していきます。
これは、乱数生成アルゴリズムが xor と shift のみで作られている以上、
内部状態の特定 1 bit のみが違う二つの乱数系列の差を比較する場合、
単純に内部状態が全て 0 の乱数系列（これは何回状態遷移しても変化しない）と、特定 1 bit のみが立った状態を初期状態とした乱数系列を比較すればそれでよい、ということです*1。
問題は、どのように観察するか、ですが、ここでは単純に、立っている bit の総数で評価します。
細かい bit の並びは置いといて、立っている bit の総数が、全状態量の bit 数の半数に近ければ、
Xorshift のアルゴリズムの特徴から、とりあえず無関係だと言っていいだろう、ということです。
で、この「とりあえず無関係といっていい」度合いの評価を、今回はエントロピー関数
H(X) = -p*log2(p) - (1-p)*log2(1-p)
を使って行ないます。 p としては大さっぱに、内部状態からランダムに選んだ bit が 1 である確率、つまり
p = (立っているbit数) / 128
を用いることとし、この確率から得られた H(X) の値を、考えられる 128 通り全ての系列で平均し、
その推移を観察することで、どの程度 読み飛ばせば「無関係だ」と言っても構わないレベルになるかを調べます。
検証コード
では、実際に測定しましょう。
今回検証に使ったコードは、以下の通りです：
// xor128.hpp （何度も使いそうなのでヘッダに分割） #ifndef INCLUDED_XOR128_HPP_ #define INCLUDED_XOR128_HPP_ #include <cassert> #include <boost/cstdint.hpp> #include <bitset> // ヘッダに分割したので一応名前空間に入れておこう namespace etude { struct Xor128 { typedef boost::uint32_t result_type; typedef boost::uint32_t int_type; // 乱数を構築する // 初期化をどうすべきかはまだ決めてないので、これのみ提供。 Xor128( int_type x_, int_type y_, int_type z_, int_type w_ ) : x(x_), y(y_), z(z_), w(w_) { assert( x != 0 || y != 0 || z != 0 || w != 0 ); } // 乱数生成 result_type operator()() { int_type const t = x ^ ( x << 11 ); x = y; y = z; z = w; w = ( w ^ ( w >> 19 ) ) ^ ( t ^ ( t >> 8 ) ); return w; } // 内部状態を取得する typedef std::bitset<128> state_type; state_type get_state() const { state_type state = x; state <<= 32; state |= state_type(y); state <<= 32; state |= state_type(z); state <<= 32; state |= state_type(w); return state; } // state_type からの構築 Xor128( state_type state ) { state_type const mask = int_type(-1); w = ( mask & state ).to_ulong(); state >>= 32; z = ( mask & state ).to_ulong(); state >>= 32; y = ( mask & state ).to_ulong(); state >>= 32; x = ( mask & state ).to_ulong(); } private: int_type x, y, z, w; }; } #endif // #ifndef INCLUDED_XOR128_HPP_
// 状態遷移をモデル化 #include "xor128.hpp" // エントロピー関数 #include <cmath> double calc_entropy( double p ) { static double const log2_ = std::log(2.0); double const p_ = 1 - p; return ( -p * std::log(p) - p_ * std::log(p_) ) / log2_; } // x の内部状態から適当に 1 bit 選んだとき、 // { 0, 0, 0, 0 } の内部状態を持つ乱数生成器の該当 bit （つまり 0 ）と // 異なっている確率のエントロピーを計算する。 double calc_entropy( etude::Xor128 const& x ) { std::size_t const all_bits = 128; // get_state() は std::bitset を返すので、 count で 1 な bit を数えられる std::size_t const high_bits = x.get_state().count(); if( high_bits == 0 || high_bits == all_bits ){ return 0; } double const p = double(high_bits) / all_bits; return calc_entropy(p); } #include <iostream> #include <vector> #include <boost/foreach.hpp> int main() { // 内部状態の bit 数 std::size_t const state_bits = 128; // 内部状態のリスト // vec[n] は、初期状態において下から n bit 目のみが立っている乱数生成器 // アルゴリズムが xor と shift のみなので、差異は { 0, 0, 0, 0 } と比較すればいい std::vector<etude::Xor128> vec; vec.reserve( state_bits ); // 初期化 { // 内部状態を表す bitset を用意。まず下 1 bit だけを立てる etude::Xor128::state_type state = 1; // state を 0…01 から順にシフトしていく for( std::size_t i = 0; i < state_bits; ++i ) { assert( state.count() == 1 ); vec.push_back( state ); assert( vec.back().get_state() == state ); state <<= 1; } } // 実際のループ std::size_t const steps = 100; for( std::size_t i = 0; i < steps; ++i ) { double sum = 0; double sum_sq = 0; // 全ての乱数生成器について統計＆状態遷移 BOOST_FOREACH( etude::Xor128& x, vec ) { double const entropy = calc_entropy(x); sum += entropy; sum_sq += entropy*entropy; // 乱数生成器の状態を次に進める x(); } // 統計結果を出力 double const mean = sum / vec.size(); double const moment2 = sum_sq / vec.size(); double const variance = moment2 - mean * mean; double const stddev = std::sqrt( variance ); std::cout << i << ", " << mean << ", " << stddev << std::endl; } return 0; }
結果
実行結果： http://ideone.com/AusHb
グラフ（ 60 step まで）： http://twitpic.com/2sus2u
考察
グラフを見る限り、内部状態 1 bit の差が全体に波及するまでには、
平均して 40 回くらい読み飛ばせば十分であることが読み取れます。
指摘にあった「40回前後で 128 bit 全ての bit が立ち得ることが分かった」というのは、恐らくこのことでしょう*2。
これは前回の日記の見積りとは大きく異なります。
その理由として考えられるのは、
「内部状態に n bit の差があった場合、乱数を 1 回生成すると、その内部状態の差異は r*n bit になる」
という乱暴な仮定が、僕の思っていた以上に荒唐無稽だった、ということです。
この仮定を持ち出すにあたり、僕は
「内部状態の差異がお互いに打ち消しあう可能性は、そこまで大きくないだろう」
と考えていましたが、実際には、少し実際の乱数系列を生成してみるとわかる通り、
殆どの bit の変化は xor によって起こされるものである以上、次の乱数生成で容易に打ち消されるため、
シフト幅が一定であることも相まって、かなり高い頻度で内部状態の変化は打ち消しあってしまうのです。
なので、この打ち消し合い分を 何らかの方法で補正しないことには、まともな見積りにはならないようです。
まとめ
というわけで、数学に詳しい人は、この辺の補正を適当に考えてみるといいと思います。
いや、僕も思いついたら書きますが、それよりは、
今回の測定手法を使って、今度は 32 bit の種から初期状態を生成したとき、
どのように初期状態を生成すれば読み飛ばし回数を減らせるか、というのを試行錯誤したいので。
その際、今回は単純に平均したわけですが、
実用上は種の下位 bit の方が上位 bit に比べて変化の度合いが大きいと考えられる以上、
種の下位 bit の変化に起因する変化に重みを付けた評価を行ってみようかな、とか思ってます。
もうちょっとだけ続くんじゃぞい。
