M系列乱数を作ってみた。
M系列乱数を作ってみた。
M系列乱数は、
R(n) = R(n-p) xor R(n-q)
の漸化式で作れる。
と言っても最初が必要なので、既存のアルゴリズムを使用して、乱数テーブルを作っておく必要がある。
nを１個ずつ進めながらxorで計算していく。nの余りを取っていけば、テーブルを増やす必要はないのでヴェクタを使って破壊的に更新していく事にする。
動作が確認しやすいように、2番目に小さいメルセンヌ素数M(3)=2^3 - 1=7を使用して、擬似乱数を作ってみます。テーブル数pに対して、2^p-1という周期が出来ることを確認したいと思います。
Gaucheでxorする時はlogxorかな。処理系によってビット演算の仕様が違うみたい。randomも無いので、srfi-27か組み込みライブラリが必要。ということで、MzScheme専用です。
(define rand (let* ((p 3) (q 1) (i 0) (l (list->vector (list-tabulate p (lambda (x) (random 2))))) (shift (lambda (x) (modulo x p)))) (lambda () (set! i (modulo (+ i 1) p)) (vector-set! l i (bitwise-xor (vector-ref l (shift (- i 3))) (vector-ref l (shift (- i 1))))) (list (vector-ref l i) i l))))
テスト。
(rand) ; (0 1 #(0 0 1)) <- ここから (rand) ; (1 2 #(0 0 1)) (rand) ; (1 0 #(1 0 1)) (rand) ; (1 1 #(1 1 1)) (rand) ; (0 2 #(1 1 0)) (rand) ; (1 0 #(1 1 0)) (rand) ; (0 1 #(1 0 0)) <- ここまでで一周 (rand) ; (0 2 #(1 0 0)) (rand) ; (1 0 #(1 0 0)) (rand) ; (1 1 #(1 1 0)) (rand) ; (1 2 #(1 1 1)) (rand) ; (0 0 #(0 1 1))
はっきり周期があるってのがおもろいな。
(random p)にして、0,1,2,3でランダムにしてみる。
(rand) ; (3 1 #(2 3 2)) (rand) ; (1 2 #(2 3 1)) (rand) ; (3 0 #(3 3 1)) (rand) ; (0 1 #(3 0 1)) (rand) ; (1 2 #(3 0 1)) (rand) ; (2 0 #(2 0 1)) (rand) ; (2 1 #(2 2 1)) (rand) ; (3 2 #(2 2 3)) (rand) ; (1 0 #(1 2 3)) (rand) ; (3 1 #(1 3 3)) ...
うっは、擬似乱数だ。
7回の周期があることが確認出来ました。1周するとインデックスが1個ズレるというのも面白い。
2^521 - 1
ちょっとまともにしてみる。
p=521,q=32を使う事にしてっと、
(define rand (let* ((rand-max 2147483647) (p 521) (q 32) (i 0) (l (list->vector (list-tabulate p (lambda (x) (random rand-max))))) (shift (lambda (x) (modulo x p)))) (lambda x (set! i (modulo (+ i 1) p)) (vector-set! l i (bitwise-xor (vector-ref l (shift (- i 3))) (vector-ref l (shift (- i 1))))) ((if (null? x) (lambda (r) (* 1.0 (/ r rand-max))) (lambda (r) (quotient r (quotient rand-max (car x))))) (vector-ref l i)))))
MzSchemeと同じように、(rand)すると、小数で、(rand 数)で整数という感じ。
(rand 100) ; 4 (rand 100) ; 83 (rand 100) ; 40 (rand 100) ; 39 (rand 100) ; 70
周期2^521 - 1がどれくらいなのか、計算してみます。
(- (expt 2 521) 1) ; 6864797660130609714981900799081393217269435300143305409394463459185543183397656 ; 052122559640661454554977296311391480858037121987999716643812574028291115057151
充分ランダムなサイコロが作れることがわかった。
もうちょい改造
循環してるんだから、循環リストで表現出来るよね。
しかも、n = pなので、R(n)とR(n-p)は、テーブル上では同じ位置にあるので、
(define rand (let* ((rand-max 2147483647) (p 521) (q 32) (l (apply circular-list (list-tabulate p (lambda (x) (random rand-max))))) (lq (drop l (- p q)))) (lambda x (set! l (cdr l)) (set! lq (cdr lq)) (set-car! l (bitwise-xor (car l) (car lq))) (if (null? x) (* 1.0 (/ (car l) rand-max)) (quotient (car l) (quotient rand-max (car x)))))))
余りを取る必要が無くなった。
だいぶスッキリ。
まぁ、フィボナッチみたいなもんだなと思って、図にしたら
全パターンが網羅できる様子を捉えた。
一番右の奴をずらして考えてみた。
A B C
V V
C AB BC
V V
BC ABC AC
V V
AC A B見事に解ける様子がわかる。
Bについて見てみると、
B = ABC xor AC = C xor AB xor AB xor BC = C xor A xor B xor A xor B xor B xor C = B
Bだけ1個多くなる。途中の分岐でうまくBが消えるんだな。
フィボナッチ恐るべし。
全パターンの数は二項定理の横列の和-1になることはわかったけど、どういう意味なんだろう不思議だ。
お、C言語による最新アルゴリズム辞典の有限体(ガロア体)の所に解説あり。mod 2の世界を構築し、521次の原始多項式を解いたと。
(x^4 + x^3 + 1) + (x^4 + x^2 + 1) = x^3 + x^2
アレ・・・計算ミスじゃないようだ。1 + 1 = 0。mod 2の世界なのだ・・・mod 2の世界で、集合{1,x,x^2...}をx^512 + x^32 + 1で割った余りになると。わ・わからん。
お、えぇもん発見。
1と0だけで割り算したり、引き算したり・・・筆算がありえなすぎる。
1000 / 1011 = 1 余り 11
圧倒的にありえない世界観。わかったけど、わからん（謎
問題点
うまい初期テーブルを出さないと、最大〜最小の値を取る事が出来ないので偏りが出てしまう。どのようにして最大、最小の値ををリストの中に組み入れるかと言う問題が残る。
が、参考書籍として挙げられている、乱数の本は、Amazon.co.jp： 乱数 (UP応用数学選書): 伏見 正則: 本1989年発行で、絶版。
お、思いついた。0とrand-maxを用意しておいて1〜rand-max-1までをp - 2個生成する。必ず最大と最小が取れるし、xorしてしまうので、ランダムになる。気になるならランダムなところに挿入すればいい。
