TinyMTの更新を行列で書いてみようのコーナー(前編)
はじめに
この記事はPokémon RNG Advent Calendar 2017 21日目の記事です。
21日の枠が空いていたので、怖い人に脅されて急遽昔やったネタを引っ張り出してきました。
これ。
twitter.comTinyMTの更新関数(遷移関数)を行列で表してみました(MAT等の値はポケモンSMで使われている値)。 状態空間(status[0]~[3])を1本の横ベクトルで表すと、更新関数f(x)はこの行列Aを用いて「f(x) = xA」と表せます。 見やすいよう1のセルを黒くしました。 pic.twitter.com/e32oMSDOkb
— 夜綱 (@sub_827) 2017年3月30日
さっさと書き上げるはずだったのですが、なんと計算したメモがロクに読めたもんじゃないという事態が発生し、一から計算し直す羽目に。
その結果クリスマスはルーズリーフと格闘して過ごしました。どうにも計画性が無さすぎて困る。
というわけで、一気に書き上げる気力が残っていないので、前後編に分かれます。
本題
去年のAdvent Calenderではおだんぽよ氏が行列を用いてTinyMTの逆算処理を求めていましたが、 この記事ではその行列を具体的に求めていきたいと思います。
更新処理の確認
まずはTinyMTの更新処理のソースコードを見てみましょう。
inline static void tinymt32_next_state(tinymt32_t * random) {
uint32_t x;
uint32_t y;
y = random->status[3];
x = (random->status[0] & TINYMT32_MASK)
^ random->status[1]
^ random->status[2];
x ^= (x << TINYMT32_SH0);
y ^= (y >> TINYMT32_SH0) ^ x;
random->status[0] = random->status[1];
random->status[1] = random->status[2];
random->status[2] = x ^ (y << TINYMT32_SH1);
random->status[3] = y;
random->status[1] ^= -((int32_t)(y & 1)) & random->mat1;
random->status[2] ^= -((int32_t)(y & 1)) & random->mat2;
}
https://github.com/MersenneTwister-Lab/TinyMTより
これを、random->status[n] = s[n]と置いたり、定数に具体的な値を入れたりして、もうちょっと読みやすく書き直すと以下のようになります。
y = s[3];
x = (s[0] & 0x7FFFFFFF) ^ s[1] ^ s[2];
x ^= (x << 1);
y ^= (y >> 1) ^ x;
s'[0] = s[1];
s'[1] = s[2];
s'[2] = x ^ (y << 10);
s'[3] = y;
if(y&1==1){
s'[1] ^= mat1 // 0x8f7011eeだけど却ってややこしくなるからそのまま
s'[2] ^= mat2 // 0xfc78ff1fだけど同上
}
基本事項の確認
この前の記事と同様に、整数を2進数表記して上のベクトルとして扱っていきます。
(最終的には内部状態を128次ベクトルとして次の状態へ移す行列を考えますが、
各statusごとに分割して処理が書かれているので、一度各statusごとの式に落とし込んでから128次ベクトルにまとめ上げていきます。)
このあたりはoupo先生の記事とかMTを作った松本眞氏の講義ノートとかを参考にしたので詳しくはそちらで。
○各statusやx,y等の32bit整数を32次横ベクトルとみなし、それぞれ,
等と表します。
○これにより、各処理は以下のように表すことができます。
・排他的論理和xor(^)
32bit整数a, bの排他的論理和を取る処理a^bは、に相当します。
・ビットシフト(<< n, >> n)
単位行列の対角上にある1を、左nビットシフトは左に、右nビットシフトは右にnずらした行列との積で表せます。
例) 11 << 2 (%16)
$$
\begin{align}
\begin{bmatrix}
1 & 0 & 1 & 1
\end{bmatrix}
\begin{bmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
\end{bmatrix}=
\begin{bmatrix}
1 & 1 & 0 & 0
\end{bmatrix}
\end{align}
$$
・論理積and(&)
a&bは、で表せます。言い換えると、aと対角上にbの各成分を並べた行列との積です。
例) 11 & 6
$$
\begin{align}
\begin{bmatrix}
1 & 0 & 1 & 1
\end{bmatrix}
\begin{bmatrix}
0 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 \\
\end{bmatrix}=
\begin{bmatrix}
0 & 0 & 1 & 0
\end{bmatrix}
\end{align}
$$
・また、if(y&1==1) a^=bは、32列目にの各要素並べた行列
を用いて、
と表せます。
例) if(y&1) a^=11を表す行列
$$
A =
\begin{align}
\begin{bmatrix}
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
1 & 0 & 1 & 1 \\
\end{bmatrix}
\end{align}
$$
さて、準備が整ったところで計算していきましょう。
1行目～4行目
・1行目 そのまんまです。
$$
{\boldsymbol y} = {\boldsymbol s}_3
$$
・2行目 に
& 0x7FFFFFFFに相当する行列を掛けます。
(の下位31bitを取るので、
と置きます。)
$$
\begin{align}
{\boldsymbol x}
& = {\boldsymbol s}_0
\begin{bmatrix}
0 & & & \\
& 1 & & \\
& & \ddots & \\
& & & 1
\end{bmatrix} +
{\boldsymbol s}_1 +
{\boldsymbol s}_2 \\
& = {\boldsymbol s}_0 U_{31} + {\boldsymbol s}_1 + {\boldsymbol s}_2
\end{align}
$$
・3行目 左辺と右辺で文字を変えるべきなのでしょうが、=を左辺への代入だと思ってください。。。
(以下、左nビットシフトを表す行列を、右nビットシフトを表す行列を
と置きます。)
$$
\begin{align}
{\boldsymbol x} &= {\boldsymbol x} + {\boldsymbol x} L_{1} \\
&= {\boldsymbol x} (I + L_{1}) \\
&= {\boldsymbol s}_{0}U_{31} (I+L_{1}) + {\boldsymbol s}_{1}(I+L_{1}) + {\boldsymbol s}_{2} (I+L_{1})
\end{align}
$$
・4行目
$$
\begin{align}
{\boldsymbol y} &= {\boldsymbol y} + {\boldsymbol y} R_{1} + {\boldsymbol x} \\
&= {\boldsymbol y} (I + R_{1}) + {\boldsymbol x} \\
&= {\boldsymbol s}_{0}U_{31} (I+L_{1}) + {\boldsymbol s}_{1}(I+L_{1}) + {\boldsymbol s}_{2} (I+L_{1}) + {\boldsymbol s}_{3} (I + R_{1})
\end{align}
$$
これでが出ました。
次はこれらをstatusへ代入していきます。
5行目以降
更新後のstatusをとします。
・5, 6行目 ただ代入するだけです。
$$
\begin{align}
{\boldsymbol s}'_{0} = {\boldsymbol s}_{1} \\
{\boldsymbol s}'_{1} = {\boldsymbol s}_{2}
\end{align}
$$
・7行目 ちょっと長ったらしくなります。見やすくするために,
と置きます。
$$
\begin{align}
{\boldsymbol s}'_{2} &= {\boldsymbol x} + {\boldsymbol y}L_{10} \\
&= {\boldsymbol s}_{0}U_{31} (I+L_{1}) + {\boldsymbol s}_{1}(I+L_{1}) + {\boldsymbol s}_{2} (I+L_{1}) \\
&\quad + ({\boldsymbol s}_{0}U_{31} (I+L_{1}) + {\boldsymbol s}_{1}(I+L_{1}) + {\boldsymbol s}_{2} (I+L_{1}) + {\boldsymbol s}_{3} (I + R_{1}))L_{10} \\
&= {\boldsymbol s}_{0}U_{31}L'_{1}L'_{10} + {\boldsymbol s}_{1}L'_{1}L'_{10} + {\boldsymbol s}_{2}L'_{1}L'_{10} + {\boldsymbol s}_{3}R'_{1}L_{10}
\end{align}
$$
・8行目
を代入するだけです。
$$
\begin{align}
{\boldsymbol s}'_{3} = {\boldsymbol s}_{0}U_{31} L'_{1} + {\boldsymbol s}_{1}L'_{1} + {\boldsymbol s}_{2} L'_{1} + {\boldsymbol s}_{3} R'_{1}
\end{align}
$$
・
if内の処理
if(y&1) s[1] ^= mat1を表す行列を、
if(y&1) s[2] ^= mat2を表す行列をと置きます。
を作用させると最下位bitは必ず0になるので、
(よって
)が成り立つことに注意します。
$$
\begin{align}
{\boldsymbol s}'_{1} &= {s}'_{1} + {\boldsymbol y}C_{mat1} \\
&= {\boldsymbol s}_{2} + ({\boldsymbol s}_{0}U_{31} L'_{1} + {\boldsymbol s}_{1}L'_{1} + {\boldsymbol s}_{2} L'_{1} + {\boldsymbol s}_{3} R'_{1})C_{mat1} \\
&= {\boldsymbol s}_{0}U_{31} L'_{1}C_{mat1} + {\boldsymbol s}_{1}L'_{1}C_{mat1} + {\boldsymbol s}_{2} (I+L'_{1}C_{mat1}) + {\boldsymbol s}_{3} R'_{1}C_{mat1} \\
&= {\boldsymbol s}_{0}U_{31} C_{mat1} + {\boldsymbol s}_{1}C_{mat1} + {\boldsymbol s}_{2} (I+C_{mat1}) + {\boldsymbol s}_{3} (I + R_{1})C_{mat1}
\end{align}
$$
$$
\begin{align}
{\boldsymbol s}'_{2} &= {\boldsymbol s}'_{2} + {\boldsymbol y}C_{mat2} \\
&= ({\boldsymbol s}_{0}U_{31}L'_{1}L'_{10} + {\boldsymbol s}_{1}L'_{1}L'_{10} + {\boldsymbol s}_{2}L'_{1}L'_{10} + {\boldsymbol s}_{3}R'_{1}L'_{10}) +({\boldsymbol s}_{0}U_{31} L'_{1} + {\boldsymbol s}_{1}L'_{1} + {\boldsymbol s}_{2} L'_{1} + {\boldsymbol s}_{3} R'_{1})C_{mat2} \\
&= {\boldsymbol s}_{0}U_{31}(L'_{1}L'_{10} + C_{mat2}) + {\boldsymbol s}_{1}(L'_{1}L'_{10} + C_{mat2}) + {\boldsymbol s}_{2}(L'_{1}L'_{10} + C_{mat2}) + {\boldsymbol s}_{3}(R'_{1}L_{10} + R'_{1}C_{mat2})
\end{align}
$$
整理してみる。
何とか処理を式に落とし込むことができたので、行列にしやすいように整理してみます。
$$
{\boldsymbol s}'_{0} = {\boldsymbol s}_{0} {\boldsymbol O} + {\boldsymbol s}_{1} + {\boldsymbol s}_{2}{\boldsymbol O}+ {\boldsymbol s}_{3}{\boldsymbol O} \\
{\boldsymbol s}'_{1} = {\boldsymbol s}_{0}U_{31} C_{mat1} + {\boldsymbol s}_{1}C_{mat1} + {\boldsymbol s}_{2} (I+C_{mat1}) + {\boldsymbol s}_{3} (I + R_{1})C_{mat1} \\
{\boldsymbol s}'_{2} = {\boldsymbol s}_{0}U_{31}(L'_{1}L'_{10} + C_{mat2}) + {\boldsymbol s}_{1}(L'_{1}L'_{10} + C_{mat2}) + {\boldsymbol s}_{2}(L'_{1}L'_{10} + C_{mat2}) + {\boldsymbol s}_{3}(R'_{1}L_{10} + R'_{1}C_{mat2}) \\
{\boldsymbol s}'_{3} = {\boldsymbol s}_{0}U_{31} L'_{1} + {\boldsymbol s}_{1}L'_{1} + {\boldsymbol s}_{2} L'_{1} + {\boldsymbol s}_{3} R'_{1}
$$
さて、これで各と
を並べたベクトルをそれぞれ作り、内部状態の更新を行列で表現すると
$$
\begin{align}
\begin{bmatrix}
{\boldsymbol s}'_{0} & {\boldsymbol s}'_{1} & {\boldsymbol s}'_{2} & {\boldsymbol s}'_{3}
\end{bmatrix} =
\begin{bmatrix}
{\boldsymbol s}_{0} & {\boldsymbol s}_{1} & {\boldsymbol s}_{2} & {\boldsymbol s}_{3}
\end{bmatrix}
\begin{bmatrix}
{\boldsymbol O} & U_{31}C_{mat1} & U_{31}(L'_{1}L'_{10}+ C_{mat2}) & U_{31} L'_{1} \\
I & C_{mat1} & L'_{1}L'_{10} + C_{mat2} & L'_{1} \\
{\boldsymbol O} & I+C_{mat1} & L'_{1}L'_{10} + C_{mat2} & L'_{1} \\
{\boldsymbol O} & (I + R_{1})C_{mat1} & R'_{1}L_{10} + R'_{1}C_{mat2} & R'_{1}
\end{bmatrix}
\end{align}
$$
となることがわかりました。
あとはこの行列を具体的に{0, 1}で表すだけですね！
後編に続く！！！ →TinyMTの更新を行列で書いてみようのコーナー(後編) - yatsuna_827’s blog
