はしくれエンジニアもどきのメモ

情報系技術・哲学・デザインなどの勉強メモ・備忘録です。

トップ > 乱数 > 疑似乱数発生法 M系列法のメモ

疑似乱数発生法 M系列法のメモ

乱数 math

疑似乱数発生法 M系列法のメモ

前回は線形合同法のメモだったが，今回はもう1つの疑似乱数発生法として， メルセンヌ・ツイスタの基になっているM系列法のメモ．

cartman0.hatenablog.com

M系列法

線形合同法とは違い1bitの乱数を出力する．

JIS規格 Z9031 5.4節にも記述されている．

一般式は以下になる．

$p$ を自然数，$c_{1}，c_{2}, \cdots ，c_{p−1}$と$x_{0}, \cdots, x_{n+p}$を0又は1として， 漸化式

x_{n+p} = c_{p-1}x_{n+p-1} \oplus c_{p-2}x_{n+p-2} \oplus \cdots \oplus c_{1}x_{n+1} \oplus x_{n}

によって生成される 0,1 の系列 $x_{1}，x_{2}, \cdots$を考える.

$x_{n＋N} ＝ x_{n}$が全ての $n$ について成り立つ最小の正の整数 $N$ を， この系列の周期という． 周期が $2^{p}−1$ であるとき，この系列のことを M系列という。

例

例： $p=3, c_{2} = 1, c_{1}=0$とすると，漸化式は

x_{n+3} = x_{n+2} + x_{n}

初期シードを$(x_{2}, x_{1}, x_{0})=(0, 0, 1)$ とすると

n=3 2 1 0
1(=$0 \oplus 1$) 0 0 1
n= 4 3 2 1
1(= $1 \oplus 0$) 1 0 0
n=5 4 3 2
1(= $1\oplus 0$) 1 1 0
n=6 5 4 3
0(= $1\oplus 1$) 1 1 1
n=7 6 5 4
1(= $0\oplus 1$) 0 1 1
n=8 7 6 5
0(= $1\oplus 1$) 1 0 1
n=9 8 7 6
0(= $0\oplus 0$) 0 1 0
n=10 9 8 7
1(= $0\oplus 1$) 0 0 1

周期は7．3bitの値が全部出現するので$2^{3}-1 = 7$となる．

かんたんな実装

seed = 0b001
print("seed:", bin(seed))

x = seed
for i in range(1,7+1):
    x_n = x & 0b1
    x_n2 = (x >> 2) & 0b1
    x_n3 = x_n ^ x_n2
    x = (x_n3 << 2) | (x >> 1)
    print(bin(x))

seed: 0b1
0b100
0b110
0b111
0b11
0b101
0b10
0b1

パラメータ$p, (c_{1},c_{2},\cdots,c_{p−1})$ の選び方

周期を$2^{p-1}にするには，$次の特性多項式（既約多項式）となるように選択する。既約というのは、それ以上因数分解できないという意味．

t^{p} + c_{p-1}t^{p-1} + \cdots + c_{1}t + 1

先程の例$p=3,c_2=1,c_1=0 $の特性多項式は以下になり因数分解できないので$2^{3}-1=7$になる．

t^3+t^2+1

既約にならない例

既約にならない例として， $c_{3}=1, c_{2}=1, c_{1}=1$を試す．

特性方程式は，

t^3+t^2+t+1=(t+1)(t^2+1)


n=3 2 1 0
1(=$0 \oplus 0 \oplus 1$) 0 0 1
n=4 3 2 1
1(=$1 \oplus 0 \oplus 0$) 1 0 0
n=5 4 3 2
0(=$1 \oplus 1 \oplus 0$) 1 1 0
n=6 5 4 3
0(=$0\oplus 1 \oplus 1$) 0 1 1
n=7 6 5 4
1(=$0\oplus 0 \oplus 1$) 0 0 1

周期は$4\neq 2^{3}-1$になる．