人気ブログランキング | 話題のタグを見る

System.Random (13) 線形合同法

System.Random 解読シリーズ

擬似乱数列をコンピュータで発生させる時、よく使われるのが線形合同法である。簡単に言うと、現在の値の次にくる乱数を計算で発生させる方法だ。現在の値の次の値を、現在の値を関数の引数にすることによって次にくる値を計算する。これを次々と行っていけば乱数の数列を作ることができる。現在の値を次の値を計算するための引数にするので、最初の値(seed) が同じならば、同じ乱数の数列を発生させることができる。

線形合同法の関数は次のような簡単な関数だ。

Xn + 1 = (A * Xn + B) mod M

計算の過程で seed と同じ値が生成されれば、それからの数列はそれまでの数列が繰り返されることになるから、線形合同法の出力は周期的になる。周期の長さや出現する値の分布を調節するためには次のようなルールがある。

1. BとMが互いに素である。
2. A-1が、Mの持つ全ての素因数で割りきれる。
3. Mが4の倍数である場合は、A-1も4の倍数である。

理論だけを調べてもわかりにくいので ghci で実験してみた。

Prelude> let next x a b m = (x * a + b) `mod` m
Prelude> next 8 3 5 13
3
Prelude> next 3 3 5 13
1
Prelude> next 1 3 5 13
8

A = 3, B = 5, M = 13 で最初の値 (seed) を8にしたら、3回繰り返すと8に戻った。したがってこの数列は 8, 3, 1 を繰り返すはずだ。同じ操作を繰り返すのは面倒なのでプログラムしてみた。

Prelude> let nexts x a b m = x : nexts (next x a b m) a b m
Prelude> take 10 $ nexts 8 3 5 13
[8,3,1,8,3,1,8,3,1,8]

A, B, M の値を工夫すると周期をMと同じにして、0〜 M-1 の数を全て出現させることができる。

A = 13, B = 5, M = 24 の場合を実験してみた。

Prelude> take 30 $ nexts 0 13 5 24
[0,5,22,3,20,1,18,23,16,21,14,19,12,17,10,15,8,13,6,11,4,9,2,7,0,5,22,3,20,1]

わかりにくいので整列させてみた。

Prelude> Data.List.sort $ take 24 $ nexts 0 13 5 24
[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]

System.Random の stdNext のソースは次のようになる。線形合同法そのままではなく、2つの乱数列を作り、その差を利用しているようだが理論的なことはわからなかった。次の乱数の発生のためには現在の乱数ではなくランダムジェネレータである整数のペアが使われている。

stdNext :: StdGen -> (Int, StdGen)
-- Returns values in the range stdRange
stdNext (StdGen s1 s2) = (fromIntegral z', StdGen s1'' s2'')
        where        z'   = if z < 1 then z + 2147483562 else z
                z    = s1'' - s2''

                k    = s1 `quot` 53668
                s1'  = 40014 * (s1 - k * 53668) - k * 12211
                s1'' = if s1' < 0 then s1' + 2147483563 else s1'
    
                k'   = s2 `quot` 52774
                s2'  = 40692 * (s2 - k' * 52774) - k' * 3791
                s2'' = if s2' < 0 then s2' + 2147483399 else s2'
by tnomura9 | 2015-03-01 18:40 | Haskell | Comments(0)
<< 吉田松陰とソクラテス グラフと排中律 >>