疑似乱数発生法 M系列法のメモ
疑似乱数発生法 M系列法のメモ
前回は線形合同法のメモだったが,今回はもう1つの疑似乱数発生法として, メルセンヌ・ツイスタの基になっているM系列法のメモ.
M系列法
線形合同法とは違い1bitの乱数を出力する.
JIS規格 Z9031 5.4節にも記述されている.
一般式は以下になる.
$p$ を自然数,$c_{1},c_{2}, \cdots ,c_{p−1}$と$x_{0}, \cdots, x_{n+p}$を0又は1として, 漸化式
によって生成される 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_{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}にするには,$次の特性多項式(既約多項式)となるように選択する。既約というのは、それ以上因数分解できないという意味.
先程の例$p=3,c_2=1,c_1=0 $の特性多項式は以下になり因数分解できないので$2^{3}-1=7$になる.
既約にならない例
既約にならない例として, $c_{3}=1, c_{2}=1, c_{1}=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$になる.
参考リンク
JIS9031「乱数生成及びランダム化の手順 Procedure for random number generation and randomization」
https://www.nagoya-bunri.ac.jp//~t-ymzm/rand/rand_index.html