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

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

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

疑似乱数発生法 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$になる.