SICPを読む(91) 3.1.2 - 乱数

Scheme

randの解説があるけど、randのソースが書いてなくてハマった。

(define random-init 12345)

(define rand
  (let ((x random-init))
       (lambda ()
         (set! x (rand-update x))
         x)))

(define (rand-update x)
  (modulo (+ (* 214013 x) 253011) 32767))

(rand) ; 10917

線形合同法なら、ax + bの余りを取る。数値は適当にやると、偏りがかなり出るので、数値に関してはCの仕様を使うといいと思う。

とりあえず進めそうだ。



他のアルゴリズムについて調べてみると、M系列やらメルセンヌ・ツイスタなんかがあるけど、案外複雑。


M系列について理解したところまでメモっておく。


pはメルセンヌ素数のM(p)を選ぶ。qの選定はイマイチ不明だけど、p > q。

M(n) = M(n-p) xor M(n-q)

  • とりあえずp個のテーブルに適当な値をセットし、xorを使って計算する。
  • 線形合同法とは違い、nは1個づつ進めるのが特徴。
  • 2^p-1回で全てのパターンを回るので、擬似乱数を作成することが出来ると。頭いい。
  • 初期値が全て0では無い限り、0以外のパターンは全て試すことになる。
  • テーブルp個に対して、2^p - 1という長い周期を持てるようになるというのは凄いかも。
  • サンプルはC言語による実用アルゴリズム辞典に載ってるけど、イマイチ理解できてない。後で写経しよう。

メルセンヌツイスタは、M系列を行列に発展させたバージョンらしい。


乱数も素数も深いな。