メルセンヌ・ツイスタの性質を理解したい。

この記事は前回の続きです。

今日のテーマ

今日のテーマは、メルセンヌ・ツイスタ(MT)の性質についてです。

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}

漸化式で表わされるということは、連続した生成された乱数をいくつか集めれば、その次に現れる数値が予測可能じゃないか！？！？！？？また逆に、今まで生成された乱数値も復元できるんじゃないか！？！？？？！

ワクワクしますね！

打倒Tempering

さて、さっそく……といいたいところですが、そういえば乱数値はTemperingとかいう操作をしてから出力していましたね。乱数を予測するには、コイツをどうにかしなければなりません。

前回はTempering行う行列 $T$ は正則だから逆行列が求まるよね？みたいな話をして終わりました。

以下の様な謎のビット演算が $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}

そういえば前回、何の脈絡もなく整数を各ビットで分けて行ベクトルとしていましたが、なんで行ベクトルを考えるのかというと、計算上都合がいいからです。例えば、XOR演算はベクトル同士の加算で表現できますし、ビットシフトは適当な正方行列との積で表現できます。

あっ！じゃあ↑の式も一つ一つのビット演算が行列として表せるじゃん！！！

ビットシフトの行列表現

これはすぐに思いつきそうです。単純に各要素をずらすだけなので、単位行列を列ベクトル分解して、それをずらしたものを掛ければよさそう。

例えば、8ビット整数について2ビットの右シフトを表す行列は以下になります。

S_8 = \left( \begin{array}{ccccc} 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{array} \right)

この記事中では、 $S_n$ をビットシフトを表す行列としますただし、 $n < 0$ のとき $|n|$ ビット左シフトを、 $n > 0$ のとき $n$ ビット右シフトを表すとします。

ANDの行列表現

これは対角成分がAND演算する整数の各ビットの値になった正方行列を考えれば良いです。例えば、8ビット整数について178 = 0b10110010とのANDを表す行列は以下のように書けます。

D_{178} = \left( \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{array} \right)

対角に10110010が現れています。なんでコレがAND演算を表すのかは、左から適当な行ベクトルを掛けて、手で計算してみればすぐにわかるはずです。

この記事中では、 $D_n$ を $n$ とのAND演算を表す行列とします。

XORの行列表現

A = 0b1001, B = 0b1001としたとき、A xor Bを考えてみます。

A	B	A xor B
1	0	1
0	1	1
0	0	0
1	1	0

各ビットについて見ると、1ビット同士の加算になっています。（桁があふれた分は無視です。）ということはどうやら、2つの整数のXOR演算は、それら整数を表す行ベクトルを単純に加算するだけで良さそうです。

Temperingの逆演算

さて、材料は揃いました。さっそく、Temperingの逆演算を表す行列を求めるとしましょう。

例のビット演算で書かれた式を行列で表してみます。

\begin{array}{l} \mathbf{y}_1 = \mathbf{x} + \mathbf{x} S_u = \mathbf{x} (I + S_u) \\ \mathbf{y}_2 = \mathbf{y}_1 + \mathbf{y}_1 S_{-s} D_b = \mathbf{y}_1 ( I + S_{-s} D_b ) \\ \mathbf{y}_3 = \mathbf{y}_2 + \mathbf{y}_2 S_{-t} D_c = \mathbf{y}_2 ( I + S_{-t} D_c ) \\ \mathbf{y}_4 = \mathbf{y}_3 + \mathbf{y}_3 S_l = \mathbf{y}_3 ( I + S_l ) \end{array}

で、式が4本もあると面倒ですし、1つにまとめてしまいましょう。

\mathbf{y} = \mathbf{x} (I + S_u) ( I + S_{-s} D_b ) ( I + S_{-t} D_c ) ( I + S_l )

ということで、ようやく $T$ の本性が分かりました。

T = (I + S_u) ( I + S_{-s} D_b ) ( I + S_{-t} D_c ) ( I + S_l )

で、コイツの逆行列 $T^{-1}$ を求めれば、それがTemperingの逆演算を表す行列です。

乱数を予想する

さて、Temperingの逆演算ができれば、話は早いですね！

前回作ったプログラムの MT内部状態 を表す配列を、先ほどの手順で復元した値N個で埋めてやれば、あとは前回説明した計算方法に従って次に出現する乱数を計算することが出来ます。やった！

乱数を復元する

じゃあ、今まで作られたであろう乱数を復元するには？

とりあえず内部状態 $\mathbf{x}$ をまず復元しなければならないワケですが、 MTの漸化式を見ると、 $x_k$ の上位ビットが $x_{k+n}$ の、下位ビットが $x_{k+n-1}$ の計算に使われてるなーって気が付くとおもいます。

じゃあ、その計算の逆をやれば、 $x_k$ が復元できるネ！っていうお話です。幸い、漸化式中に現れる演算はさっきやったXOR演算と、既に中身が分かっている行列 $A$ の乗算だけです！

Aはこんな行列なのでした。

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)

この定義を見ると、 $a_{w-1} \ne 0$ ならば明らかに正則ですねコレ。 MT19937では↑を満たしますし、逆行列も求まりますね！

元の論文には、高速に計算するためにこの形にするって書いてあるんですけど、コレ正則じゃなくても良いのかな？そうだとすると逆演算ができない？？？（知らない）

やってみよう！

早速やりましょう。とりあえず、前回作ったプログラムで生成した乱数をファイルに書き出して、それを入力して、その部分乱数列から全体を復元しましょう。

はい、やり方はさっきまでさんざん書いたとおりなので、サクッと実装します。今回もRubyです。

require "matrix"

# 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)

# ビット行列
class BitMat < Matrix
	# 行ベクトルを整数に
	def to_i
		# mod 2 しないとダメ
		self.row(0).to_a.map{|i| i.to_i.abs % 2 }.join.to_i(2)
	end
	# 整数を行ベクトルに
	def self.from_i(i)
		self[format("%." + W.to_s + "b", i).split("").map(&:to_i)]
	end

	# kだけビットシフトする演算を表す行列を生成 (k > 0 : 右シフト, k < 0 : 左シフト)
	def self.Shift(k)
		self[ *(1..W).map{|i| (1..W).map{|j| j == i+k ? 1 : 0 } } ]
	end
	# kとのAND演算を表す行列を生成
	def self.And(k)
		self.diagonal(*self.from_i(k).row(0).to_a)
	end

	# r行目を行ベクトルvで置き換え
	def []=(r, v)
		@rows[r] = v.row(0).to_a
	end
end

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

# 行列T
t = (
	(BitMat.I(W) + BitMat.Shift( U)) *
	(BitMat.I(W) + BitMat.Shift(-S) * BitMat.And(B)) *
	(BitMat.I(W) + BitMat.Shift(-T) * BitMat.And(C)) *
	(BitMat.I(W) + BitMat.Shift( L))
)
# Tの逆行列
t_inv = t ** -1

# 行列A
a = BitMat.Shift(1)
a[W-1] = BitMat.from_i(A)
# Aの逆行列
a_inv = a ** -1

# 乱数列を読み込む
input = []
while line = gets
	input.push($1.to_i) if /^\d+: (\d+)$/ =~ line
end

# E番目からN個だけの乱数を使う
E = input.length / 3

print "USE: " + E.to_s + " -> " + (E + N - 1).to_s + "\n";

# E番目～E+N-1番目の乱数を取得し、Temperingの逆演算をする
xr = input[E, N].map{|e| (BitMat.from_i(e) * t_inv).to_i }

# 内部状態
x = xr.dup

# E+N番目から順に乱数列を復元
(E + N).upto(input.length - 1) do |k|
	i = (k - E) % N

	# 乱数を計算
	z = x[i] & UPPER_MASK | x[(i + 1) % N] & LOWER_MASK
	x[i] = x[(i + M) % N] ^ (BitMat.from_i(z) * a).to_i

	# Temperingして一致するか確認
	y = (BitMat.from_i(x[i]) * t).to_i
	if y != input[k]
		abort "FAIL: " + k.to_s + "\n";
	end
end

print "RETRIEVED: " + (E + N).to_s + " -> ", input.length - 1,"\n";

# 内部状態
x = xr.dup

# E-1番目から順に0番目までの乱数列を復元
(E - 1).downto(0) do |k|
	i = (k - E) % N

	# z_i を復元
	z = (BitMat.from_i(x[i] ^ x[(i + M) % N]) * a_inv).to_i

	# z_{i-1} を復元
	zp = (BitMat.from_i(x[(i - 1 + N) % N] ^ x[(i - 1 + M) % N]) * a_inv).to_i

	# z_i, z_{i-1} から x[i] を復元
	x[i] = z & UPPER_MASK | zp & LOWER_MASK

	# Temperingして一致するか確認
	y = (BitMat.from_i(x[i]) * t).to_i
	if y != input[k]
		abort "FAIL: " + k.to_s + "\n";
	end
end

print "RETRIEVED: 0 -> " + (E - 1).to_s + "\n";

今回は、MTで乱数を生成する際にも行列を用いて計算してみました。なのでとても遅いですね。

ちなみに、わざわざ行列を使わなくてもTemperingの逆演算をビット演算で高速に行うことも出来ます。

MTはダメ？

わー！MTで作った乱数は予想・復元されちゃう！危ない！

というわけではなく、要は適材適所、その性質をよく理解して使うべき、ということでした。決して暗号用途に使っちゃダメですよ。（そういうCTF問題がどこかにありましたね……！）

次は xorshift を調べてみようかな？

6715.jp

メルセンヌ・ツイスタを倒す