前回の日記で Xorshift アルゴリズムについて書いた時、最後に
「初期化された後にしばらく乱数列を読み飛ばせば偏りは簡単に回避できる」
と書きましたが、では偏りを回避するには、何回くらい読み飛ばせばよいのか、少し考えてみました。
で、それを考えるにあたって、まず「偏りがない」とはどういうことかを考える必要がありますが、今回は
「乱数の種（≒内部状態）が僅かにしか違わなくても、生成される乱数系列が大きく変わる」
場合に「偏りがない」とすることにします。
言い換えるなら、
「乱数の種を 1 bit でも変化させれば、それによって大きく異なる乱数系列を得られる」
ような時に、乱数の種による偏りがない、と表現することにし、
そのようにするには 何回 乱数列を読み飛ばせばいいかを考えてみる、ということです。
また、話を単純にするため、今回は単純に「内部状態」のみを扱い、
与えられた種から内部状態を作る、という作業は特に考慮しないことにします。
つまり、内部状態が僅かに異なる二つの乱数生成系列があったとき、
それぞれの乱数系列を何回進めれば 内部状態の僅かな差異が増幅されて
お互いが「殆ど無関係」と言っても構わないようになるか、ということを考えることにするのです。
では、それを前提に、具体的に考えてみましょう。
といっても、そのあたりを厳密に計算するのは いろいろと面倒です。
なので、そのあたりは数学の得意な方に任せるとして、
今回は適度にいいかげんな「まぁ大体そう言っていいでしょ？」というレベルの見積りをしてみます。
具体的には、まず、割と簡単に計算できる、
「内部状態の 1 bit の差異が、乱数生成 1 回で、平均して何 bit に影響するか？」
というのを求めます。
んで、その得られた値を仮に r と置き、ここで ものすごく乱暴な仮定をします。
それは、
「内部状態に n bit の差があった場合、乱数を 1 回生成すると、その内部状態の差異は r*n bit になる」
という仮定です。かなりデタラメなトンデモ仮定ですが、大雑把な見積もりなので許してください。
で、その仮定を使うと、乱数生成アルゴリズムの内部状態量を N bit としたとき、
「最初の x bit の差異が乱数生成アルゴリズムの内部状態 N bit 全てに影響するまでに必要な、平均の回数 n(x, N) 」
を、 n(x, N) = log(N/x) / log(r) で計算することができるので、
それをもって解を得られたことにしよう、ということです。乱暴ですがまぁ悪く無いと思います、たぶん。
では実際に計算してみましょう。
今回は 128bit の Xorshift アルゴリズムのなので、内部状態 N は 128 です。 x は順当に 1 bit とします。
んで、問題の「内部状態 1 bit の差異が、乱数生成 1 回で（平均）何 bit に影響するか」ですが、
これはアルゴリズムを睨みながら考えると、
int_t const t = x ^ ( x << 11 ); x = y; y = z; z = w; w = ( w ^ ( w >> 19 ) ) ^ ( t ^ ( t >> 8 ) );
- y, z に差異があった場合には、差異は移動するだけで 1 bit のまま（確率 1/2 ）
- w に差異があった場合には、z に移動される他、次の w にも影響する（確率 1/4 ）
- x に差異があった場合には、生成される t が影響され、それにより次の w に影響する（確率 1/4 ）
ということが分かります。
このうち最初のケースはすぐ分かるので、 w と x に差異があった場合のことを考てみましょう。
まず w に 1 bit の差異がある場合、その影響は次の w を求める式
w = ( w ^ ( w >> 19 ) ) ^ ( t ^ ( t >> 8 ) );
から、 w のビット列のどの部分に差異があるかによって、二通りに分かれます。
w の下位 19 bitの何処かに差異がある場合には、ビットシフトで xor の片方の差異は消されてしまうので、
それが原因で w にもたらされる差異は 1 bit のみとなり、それと z に影響する 1 bit で、合計 2 bit の差異が生じます。
そうでない場合には w ^ ( w >> 19 ) によって 2 bit 分（と z に影響する 1 bit）の差異が生じるので、
平均して 2 * 19/32 + 3 * ( 1 - 19/32 ) = 2.40625 bit の差異がある、と見積もれます。
次に x に 1 bit の差異があった場合ですが、そのことに起因する差異は、
t = x ^ ( x << 11 ) と w = ( w ^ ( w >> 19 ) ) ^ ( t ^ ( t >> 8 ) ) より、
- x の上位 11 bit のいずれかに差異がある場合と
- x の下位 8 bit のいずれかに差異がある場合
- それ以外の場合
の 3 通りに分けて考えることができます。
まず x の上位 11 bit に差異があった場合ですが、これは t = x ^ ( x << 11 ) では増幅されず、
次の t ^ ( t >> 8 ) ではそれが倍に増幅されるので、差異は 2 bit になります。
次に x の下位 8 bit に差異があった場合は、まず t = x ^ ( x << 11 ) で増幅され 2 bit になりますが、
次の t ^ ( t >> 8 ) の右辺で下位 8 bit にある差異が消されてしまうため、差異は 3 bit です。
最後に上記のどちらでもない場合には、 t = x ^ ( x << 11 ) で差異が 2 bit に増幅され、
それが t ^ ( t >> 8 ) で更に倍になりますから、 4 bit の差異が生じることになります。
これらを平均すると、 x のいずれかの bit に 1 bit の差異があった場合、
それによって生じる内部状態の差異は、平均して
2 * 11/32 + 3 * 8/32 + 4 * ( 1 - 11/32 - 8/32 ) = 3.0625 bit と見積もることができます。
後は、これらから得られた値を平均すれば、
1/2 * 1 + 1/4 * 2.40625 + 1/4 * 3.0625 = 1.8671875 bit
の差異が、 1 bit 分の内部状態の差異に対して生じる、と見積もれることになります。
後はこの値を n(x, N) = log(N/x) / log(r) に代入することで、 n ≒ 7.77 という結果が得られます。
つまり、これは平均の場合なので目安でしかないですが、
xor128 の場合は、ざっと 8 回ほど乱数を読み飛ばせば、初期状態に起因する偏りはそれなりに取り除ける、
という見積りを得ることができたのです。めでたしめでたし。
まぁ実際には間違ってるんですが。
検証結果はこちら： http://d.hatena.ne.jp/gintenlabo/20100928/1285702435
どこが間違っているのか、どの程度間違ってるのか、
まだ検証結果を読んでないなら、考えてみると面白いかもしれませんね。
まぁ実際にはこれは平均でしかないので、万全を期すためには 最悪の場合を見積る必要がありそうですが、
その場合（要するに z の上位 11 bit に差異がある場合）については、読者への宿題ということにします。
最悪の場合は種から内部状態を作る段階で取り除けるので、ここで計算する必要もないでしょうし。
