メルセンヌ・ツイスタについて。

メルセンヌ・ツイスタ(MT)は擬似乱数列を作るアルゴリズムの一つで、他の手法と比べると欠点が少なくて高品質な擬似乱数列を高速に作れるんだって。スゴイ！プログラムをかじった人なら、多分聞いたことがあるんじゃないかと思います。

日本初のアルゴリズムだし、日本語文献あるかな？って思ったんですけど、良い物が見つからなかった（かなしい）ので、元の論文を読みながらRubyでMTを実装して理解を深めたいと思います。

なんか強いアルゴリズムだ！って聞くとめっちゃ複雑なんじゃないかって思いがちですけど、MTはとっても†単純†です。

MTの定義

$w$ ビットの整数からなる乱数列を生成する場合を考えます。この時、整数を各ビットで分けて $w$ 次元行ベクトルとして考えることとします。

すると、MTによって生成される乱数列は以下の線形漸化式によって表されます。

\begin{array}{c} \mathbf{x}_{k+n} \, = \, \mathbf{x}_{k+m} \oplus ( \mathbf{x}^u_k \: | \: \mathbf{x}^l_{k+1} ) A \\ (k = 0,1, \cdots) \end{array}

この式に登場する $n,m$ は定数で、 $1 \le m \le n$ を満たします。

$( \mathbf{x}^u_k \: | \: \mathbf{x}^l_{k+1} )$ は、 $\mathbf{x}_k$ の上位 $w-r$ ビットと $\mathbf{x}_{k+1}$ の下位 $r$ ビットを連結した行ベクトルを表しています。この $r$ も定数で、 $0 \le r \le w-1$ を満たします。

$A$ は以下のように定義された $w \times w$ 正方行列です。

A = \left( \begin{array}{ccccc} 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & \ddots & 0 \\ 0 & 0 & 0 & 0 & 1 \\ a_{w-1} & a_{w-2} & \cdots & \cdots & a_0 \end{array} \right)

これによって、 $\mathbf{x} A$ を以下のように高速に計算することができます。

\begin{array}{c} \mathbf{x} A = \begin{cases} (\mathbf{x} \gg 1) & (x_0 = 0) \\ (\mathbf{x} \gg 1) \oplus \mathbf{a} & (x_0 = 1) \end{cases} \\ \\ \mathbf{x} = (x_{w-1}, x_{w-2}, \cdots, x_0) \\ \mathbf{a} = (a_{w-1}, a_{w-2}, \cdots, a_0) \end{array}

こうして、 $w$ ビットの整数、もとい $w$ 次元行ベクトルがたくさんできるわけですが、これをそのまま乱数として出力するのはなんだかマズいらしくて、値を程よく分布させるため、出力する行ベクトルに $w \times w$ の適当な正則行列 $T$ を右から掛けます。

$T$ を右から掛ける事に相当する演算として、実際には以下の様な演算を行います。

\begin{array}{l} \mathbf{y}_1 = \mathbf{x} \oplus (\mathbf{x} \gg u) \\ \mathbf{y}_2 = \mathbf{y}_1 \oplus ((\mathbf{y}_1 \ll s) \: \& \: \mathbf{b}) \\ \mathbf{y}_3 = \mathbf{y}_2 \oplus ((\mathbf{y}_2 \ll t) \: \& \: \mathbf{c}) \\ \mathbf{y}_4 = \mathbf{y}_3 \oplus (\mathbf{y}_3 \gg l) \end{array}

$u,s,t,l$ は定数で、&はビットAND演算を表し、 $\mathbf{b}, \mathbf{c}$ は適当な行ベクトルです。

こうして得られた $\mathbf{y_4}$ を出力します。（この操作を Tempering と言うそうです。）

で、途中にいろいろ定数やらなんやらが登場したんですが、これらを

\begin{array}{l} w = 32 \\ n = 624 \\ m = 397 \\ r = 31 \\ u = 11 \\ s = 7 \\ t = 15 \\ l = 18 \\ \mathbf{a} = \mathbf{0x9908B0DF} \\ \mathbf{b} = \mathbf{0x9D2C5680} \\ \mathbf{c} = \mathbf{0xEFC60000} \end{array}

とすると、周期が $2^{19937}-1$ とめちゃくちゃ長い、かの有名なMT19937になります。

MTを実装する

なんだか、は？ってカンジですね！

整理すると、

長さnの配列を用意して、適当な値で埋める（これがシード値になります）
i番目の乱数を得る
- Step.1 z = x[i] & 0b111..1000..0 | x[(i+1)%n] & 0b000..0111..1 を計算する
- Step.2 x[i] = x[(i+m)%n] ^ (z >> 1) ^ (z & 1 == 0 ? 0 : a) を計算する
- Step.3 Temperingした値を返す

これだけです。シンプルだ！

Step.1は $( \mathbf{x}^u_k \: | \: \mathbf{x}^l_{k+1} )$ を求めることに相当します。 Step.2は最初の漸化式を適用することに相当し、XORで繋がれた後ろの2項は、行列 $A$ を掛けることに相当します。

自分みたいなプログラマな人間は、たぶんソースコードを見れば一発で理解できるんじゃないかと思います。ということで、書いてみたのが以下。

# MT19937
W = 32
N = 624
M = 397
R = 31
U = 11
S = 7
T = 15
L = 18
A = 0x9908B0DF
B = 0x9D2C5680
C = 0xEFC60000

# ビットマスク用
WHOLE_MASK = ("1" * W).to_i(2)
UPPER_MASK = ("1" * (W - R) + "0" * R).to_i(2)
LOWER_MASK = ("0" * (W - R) + "1" * R).to_i(2)

# MT乱数列
class MT
	# seedを受け取って初期化
	def initialize(seed)
		# MT内部状態
		@i = 0
		@x = [seed & WHOLE_MASK]

		# 初期化 (mt19937ar.cに準拠)
		1.upto(N-1) do |i|
			@x[i] = (1812433253 * (@x[i-1] ^ (@x[i-1] >> 30)) + i) & WHOLE_MASK
		end
	end

	# MTで乱数を生成
	def next
		# Step.1
		z = @x[@i] & UPPER_MASK | @x[(@i + 1) % N] & LOWER_MASK

		# Step.2
		@x[@i] = @x[(@i + M) % N] ^ (z >> 1) ^ (z & 1 == 0 ? 0 : A)

		# Step.3
		y = @x[@i]
		y = y ^ (y >> U)
		y = y ^ ((y << S) & B)
		y = y ^ ((y << T) & C)
		y = y ^ (y >> L)

		# カウンタを変更して、生成した乱数を返す
		@i = (@i + 1) % N
		return y
	end
end

########################################

# 使ってみる
mt = MT.new(20150919) # ← シード値
2048.times do |i|
	print i, ": ", mt.next, "\n"
end

（あんまり自信がないですが、MTの考案者が書いたC++プログラムと出力が一致したので、多分大丈夫。）

こうやってみると、本当に単純ですね！それでいて性能がたいへんよろしいのですから、スゴイものです。

ちなみに↑はもともとPerlで書いてたんですが、諸事情によりRubyで書き直しました。以下はPerl版。

発展

とりあえずMTがこういうものだって、わかった気になれたわけです。個人的には、 Tempering のトコロが面白いと思っていて、ココをもうちょっと掘り下げてみたいと思います。

Temperingは正則行列 $T$ を右から掛ける演算ですが、この $T$ が実際にはどんな行列なのかには、触れませんでした。

でも、正則行列ってことは逆行列が存在して、Temperingの逆演算もできて……？みたいなお話です。（なんか楽しそうな気がしません？）

ということで、次回に続く！！！

6715.jp

メルセンヌ・ツイスタをわかった気になる

MTの定義

MTを実装する

発展