線形合同法で色々やってみる。

任意の未来の状態を計算できる乱数について調べてたら、線形合同法ならパラメータを計算し直すだけで何とかならないかと考えた。あと 2ch のログでもこういうのを発見した。(出典: http://unkar.org/r/tech/1192628099 [擬似乱数２ - プログラム技術@2ch掲示板])

140 ：デフォルトの名無しさん[sage]：2008/09/17(水) 19:03:50

XorShift ならちょっと考えればわかる。

線形合同法なら 2^{31}-1 だけ進めた x=Ax+C の A と C を求めれば可能。

メルセンヌツイスタやWELLなら

http://www.iro.umontreal.ca/~lecuyer/myftp/papers/jumpf2.pdf

の方法で周期よりも一つ少ないところへジャンプすれば不可能ではないはず。

というのも、訳あって SIMD 版と非 SIMD 版で全く同じ乱数を使用する必要があり、しかもある程度高速な方が好まれるという事情があるためだ。(あと本来の乱数性もそれほど必要ではない。) 線形合同法は、念の為に説明すると次の形の乱数生成法だ。(パラメータの制約などについては Wikipedia を参照してほしい。)

$X_{n+1} = A ￥cdot X_n + B ￥pmod{M}$

2 段重ね

では、これを二段重ねにすることを考えてみよう。剰余に関連する定理により、

$X_{n+2} = A ￥cdot X_{n+1} + B ￥pmod{M}$

$X_{n+2} = A(A ￥cdot X_n + B) + B ￥pmod{M}$

の両者は等価である。これをさらに展開すると面白いことになる。

$X_{n+2} = P_2 ￥cdot X_n + Q_2 ￥pmod{M}$

$P_2 = A^2$

$Q_2 = AB + B$

となる。つまり、線形合同法のパラメータを変更するだけで、2 段飛ばしの乱数を生成できるようになるのだ。

n 段重ね

さらに一般化しよう。 $k ￥geq 1$ について、具体的な計算は省略するが、次のようになる。

$X_{n+k} = P_k ￥cdot X_n + Q_k ￥pmod{M}$

$P_k = A^k ￥pmod{M}$

$Q_k = B ￥sum^{k-1}_{m=0} A^m ￥pmod{M}$

2^n 段重ね, k+l 段重ね

このままだと $Q_k$ を計算するコストが結構ある。これをカバーするには、二段重ねを応用して $2^n$ 段重ねをやって、これを最終的に合成すれば比較的高速に生成できる。2 つの n 段重ねの系列を合成 (k+l 段重ね) するためには、次の計算をすれば良い。

$X_{n+k+l} = P_{k+l} ￥cdot X_n + Q_{k+l} ￥pmod{M}$

$P_{k+l} = P_k P_l ￥pmod{M}$

$Q_{k+l} = P_k Q_l + Q_k ￥pmod{M}$

これを使ってできるのが、上にある 2ch ログで答えがあったような、「線形合同法による乱数列の巻き戻し」だ。周期が f の乱数列の場合、f-1 回後の乱数を生成するパラメータを生成すれば、乱数列の巻き戻しが可能になる。(一般的には周期が $2^{32}$ や $2^{64}$ となるから、その周期の一つ前を計算すれば良いわけだ。)