読者です 読者をやめる 読者になる 読者になる

メルセンヌ・ツイスタの逆算に関する考察

数学 ブール代数 乱数 乱数調整

メルセンヌ・ツイスタを逆算する式を求める。
なお、 tempering 部分の逆関数は既にできているので省略する(→メルセンヌ・ツイスタのtemperingの逆関数に関する考察 - Plus Le Blog)。



32bitブール代数で考える。
以下、定数や変数はすべて整数とし、特に、添字などに使われている場合などの例外を除いて32bit非負整数とする。


※本記事の数式中で積記号が省略されている場合、算術積 \(\cdot\) ではなくすべて論理積 \(\cap\) である。
算術積を使う場合は、必ず \(x\cdot y\) のように省略せず明示するものとする。

※本記事では、「最下位ビット」を“第1ビット”,「最上位ビット」を“第32ビット”とする。
後に重要になるので念のため書いておくと、“第31ビット”は「最上位ビットの1つ下のビット」を指し、“下位30ビット”は「第1ビットから第30ビットまで」を指す。


まず、メルセンヌ・ツイスタを数式化してみる。
最初に準備として、メルセンヌ・ツイスタで使われる定数値を定義しておく。

(*1)
\begin{eqnarray*}
N &\overset{\mathrm{def}}{=}& 624\\
M &\overset{\mathrm{def}}{=}& 397\\
A &\overset{\mathrm{def}}{=}& \mathtt{0x9908b0df}\\
M_{u} &\overset{\mathrm{def}}{=}& \mathtt{0x80000000}\\
M_{l} &\overset{\mathrm{def}}{=}& \mathtt{0x7fffffff}
\end{eqnarray*}

※ \(M_{u},\,M_{l}\) はそれぞれ upper mask, lower mask の略であり、どちらも \(M\left(=397\right)\) とは一切関係ないので注意。
※定数 \(A\) も含め、これらの変数は元の論文の命名を参考にしている。
ちなみに元の論文では \(x\) の上位(1?)ビットと \(y\) の下位(31?)ビットの論理和を \(\left(x^{u}|y^{l}\right)\) のように表記していたが、分かりにくいので不採用。

さて、メルセンヌ・ツイスタを漸化式で書くと、次式のように表される。

(*2)
\begin{eqnarray*}
x_{i}&\overset{\mathrm{def}}{=}&
\begin{cases}
\text{初期値(任意の値)} & \left(0\le i< N\text{のとき}\right)\\
x_{i+M-N}\oplus\left(y_{i}>>1\right)\oplus\left(\left(y_{i}\cap 1\right)\cdot A\right) & \left(i\ge N\text{のとき}\right)\\
\end{cases}\\
y_{i}&\overset{\mathrm{def}}{=}&M_{u}x_{i-N}\cup M_{l}x_{i+1-N}
\end{eqnarray*}

この式(*2)を、逆算しやすい形に変形していこう。
まず、 \(y_{i}\) を消去するため、 \(y_{i}\) を変形してみる。
\(\alpha\cup\beta=\alpha\oplus\beta\oplus\alpha\beta\) より、

(*3)
\begin{eqnarray*}
y_{i}
&=&M_{u}x_{i-N}\cup M_{l}x_{i+1-N}\\
&=&M_{u}x_{i-N}\oplus M_{l}x_{i+1-N}\oplus M_{u}x_{i-N} M_{l}x_{i+1-N}\\
&=&M_{u}x_{i-N}\oplus M_{l}x_{i+1-N}\oplus x_{i-N}\cap x_{i+1-N}\cap\mathtt{0x80000000}\cap\mathtt{0x7fffffff}\\
&=&M_{u}x_{i-N}\oplus M_{l}x_{i+1-N}\oplus x_{i-N}\cap x_{i+1-N}\cap 0\\
&=&M_{u}x_{i-N}\oplus M_{l}x_{i+1-N}
\end{eqnarray*}

このとき \(y_{i}\cap 1\) について、 \(M_{u}\cap 1=\mathtt{0x80000000}\cap 1=0\) , \(M_{l}\cap 1=\mathtt{0x7fffffff}\cap 1=1\) なので、

(*4)
\begin{eqnarray*}
y_{i}\cap 1
&=&\left(M_{u}x_{i-N}\oplus M_{l}x_{i+1-N}\right)\cap 1\\
&=&M_{u}x_{i-N}\cap 1\oplus M_{l}x_{i+1-N}\cap 1\\
&=&x_{i+1-N}\cap 1
\end{eqnarray*}

よって(*2)の漸化式部分を書き換えると、式(*5)のようになる。

(*2)再掲
\[
x_{i}=x_{i+M-N}\oplus\left(y_{i}>>1\right)\oplus\left(\left(y_{i}\cap 1\right)\cdot A\right)\qquad\left(i\ge N\text{のとき}\right)
\]

(*5)
\(i\ge N\) なる \(i\) について、
\begin{eqnarray*}
x_{i}
&=&x_{i+M-N}\oplus\left(y_{i}>>1\right)\oplus\left(\left(y_{i}\cap 1\right)\cdot A\right)\\
&=&x_{i+M-N}\oplus\left(\left(M_{u}x_{i-N}\oplus M_{l}x_{i+1-N}\right)>>1\right)\oplus\left(\left(x_{i+1-N}\cap 1\right)\cdot A\right)\\
&=&x_{i+M-N}\oplus\left(\left(M_{u}x_{i-N}\right)>>1\right)\oplus\left(\left(M_{l}x_{i+1-N}\right)>>1\right)\oplus\left(\left(x_{i+1-N}\cap 1\right)\cdot A\right)
\end{eqnarray*}

さらに、 \(i\rightarrow i+N\) の変換を行うと、

(*6)
\(i\ge 0\) なる \(i\) について、
\[
x_{i+N}=x_{i+M}\oplus\left(\left(M_{u}x_{i}\right)>>1\right)\oplus\left(\left(M_{l}x_{i+1}\right)>>1\right)\oplus\left(\left(x_{i+1}\cap 1\right)\cdot A\right)
\]

よって右辺を左辺に移項することで、次式を得る。

(*7)
\[
x_{i+N}\oplus x_{i+M}\oplus\left(\left(M_{u}x_{i}\right)>>1\right)\oplus\left(\left(M_{l}x_{i+1}\right)>>1\right)\oplus\left(\left(x_{i+1}\cap 1\right)\cdot A\right)=0
\]

式(*7)を図で表すと次のようになる。
※画像はイメージです。

式(*7)の分解

さて今、より分かりやすくするために式(*7)を「最上位ビット(第32ビット)」「第31ビット」「下位30ビット」の3つに分解することを考える。

まず最上位ビットに注目する。
式(*7)において、 \(\left(M_{u}x_{i}\right)>>1\) および \(\left(M_{l}x_{i+1}\right)>>1\) の項については、右にビットシフトしているので最上位ビットは \(0\) である。

(*8)
\begin{eqnarray*}
\left(\left(\left(M_{u}x_{i}\right)>>1\right) \text{の最上位ビット}\right) &=& 0\\
\left(\left(\left(M_{l}x_{i+1}\right)>>1\right) \text{の最上位ビット}\right) &=& 0
\end{eqnarray*}

また、 \(\left(x_{i+1}\cap 1\right)\cdot A\) の項は、 \(x_{i+1}\) が偶数のとき \(0\) ,奇数のとき \(A\) になる。
\(A\left(=\mathtt{0x9908b0df}\right)\) の最上位ビットは \(1\) なので、 \(\left(x_{i+1}\cap 1\right)\cdot A\) の最上位ビットは \(x_{i+1}\) の最下位ビットに等しい。

(*9)
\[
\left(\left(\left(x_{i+1}\cap 1\right)\cdot A\right) \text{の最上位ビット}\right) = \left(x_{i+1} \text{の最下位ビット}\right)
\]

以上より、式(*7)の左辺の最上位ビットは、式(*10)のようになる。

(*7)再掲
\[
x_{i+N}\oplus x_{i+M}\oplus\left(\left(M_{u}x_{i}\right)>>1\right)\oplus\left(\left(M_{l}x_{i+1}\right)>>1\right)\oplus\left(\left(x_{i+1}\cap 1\right)\cdot A\right)=0
\]

(*10)
\begin{eqnarray*}
&& \left(\text{式(*7)の左辺の最上位ビット}\right)\\
&=& \left(\left(x_{i+N}\oplus x_{i+M}\oplus\left(\left(M_{u}x_{i}\right)>>1\right)\oplus\left(\left(M_{l}x_{i+1}\right)>>1\right)\oplus\left(\left(x_{i+1}\cap 1\right)\cdot A\right)\right)\text{の最上位ビット}\right)\\
&=& \left(\left(x_{i+N}\oplus x_{i+M}\right)\text{の最上位ビット}\right)\oplus\left(x_{i+1} \text{の最下位ビット}\right)
\end{eqnarray*}

これが式(*7)の右辺の最上位ビット、すなわち \(0\) に等しいので、

(*11)
\[
\left(\left(x_{i+N}\oplus x_{i+M}\right)\text{の最上位ビット}\right)\oplus\left(x_{i+1} \text{の最下位ビット}\right)=0
\]

適当に移項して変形すると、

(*12)
\begin{eqnarray*}
\left(x_{i+1} \text{の最下位ビット}\right)
&=&\left(\left(x_{i+N}\oplus x_{i+M}\right)\text{の最上位ビット}\right)\\
&=&\left(x_{i+N}\oplus x_{i+M}\right)>>31
\end{eqnarray*}


次に、式(*7)の第31ビットについて考える。
式(*7)において、 \(\left(M_{u}x_{i}\right)>>1\) の項については、 \(M_{u}=\mathtt{0x80000000}\) なので、

(*13)
\[
\left(\left(\left(M_{u}x_{i}\right)>>1\right) \text{の第31ビット}\right)=\left(x_{i} \text{の最上位ビット}\right)
\]

また、 \(\left(M_{l}x_{i+1}\right)>>1\) の項については、 \(M_{l}=\mathtt{0x7fffffff}\) なので、

(*14)
\[
\left(\left(\left(M_{l}x_{i+1}\right)>>1\right) \text{の第31ビット}\right)=0
\]

さらに、 \(\left(x_{i+1}\cap 1\right)\cdot A\) の項は、 \(x_{i+1}\) が偶数のとき \(0\) ,奇数のとき \(A\) になる。
\(A\left(=\mathtt{0x9908b0df}\right)\) の第31ビットは \(0\) なので、 \(\left(x_{i+1}\cap 1\right)\cdot A\) の第31ビットは \(x_{i+1}\) の偶奇に依らず常に \(0\) である。

(*15)
\[
\left(\left(\left(x_{i+1}\cap 1\right)\cdot A\right) \text{の第31ビット}\right)=0
\]

以上より、式(*7)の左辺の第31ビットは、式(*16)のようになる。

(*7)再掲
\[
x_{i+N}\oplus x_{i+M}\oplus\left(\left(M_{u}x_{i}\right)>>1\right)\oplus\left(\left(M_{l}x_{i+1}\right)>>1\right)\oplus\left(\left(x_{i+1}\cap 1\right)\cdot A\right)=0
\]

(*16)
\begin{eqnarray*}
&& \left(\text{式(*7)の左辺の第31ビット}\right)\\
&=& \left(\left(x_{i+N}\oplus x_{i+M}\oplus\left(\left(M_{u}x_{i}\right)>>1\right)\oplus\left(\left(M_{l}x_{i+1}\right)>>1\right)\oplus\left(\left(x_{i+1}\cap 1\right)\cdot A\right)\right)\text{の第31ビット}\right)\\
&=& \left(\left(x_{i+N}\oplus x_{i+M}\right)\text{の第31ビット}\right)\oplus\left(x_{i} \text{の最上位ビット}\right)
\end{eqnarray*}

これが式(*7)の右辺の第31ビット、すなわち \(0\) に等しいので、

(*17)
\[
\left(\left(x_{i+N}\oplus x_{i+M}\right)\text{の第31ビット}\right)\oplus\left(x_{i} \text{の最上位ビット}\right)=0
\]

移項すると、

(*18)
\[
\left(x_{i} \text{の最上位ビット}\right)=\left(\left(x_{i+N}\oplus x_{i+M}\right)\text{の第31ビット}\right)
\]

よって次式が成り立つ。

(*19)
\begin{eqnarray*}
\left(x_{i} \text{の最上位ビット}\right)<<31
&=&\left(\left(x_{i+N}\oplus x_{i+M}\right)\text{の第31ビット}\right)<<31\\
&=&M_{u}\left(\left(x_{i+N}\oplus x_{i+M}\right)<<1\right)
\end{eqnarray*}


最後に、式(*7)の下位30ビットについて考える。
式(*7)において、 \(\left(M_{u}x_{i}\right)>>1\) の項については、 \(M_{u}=\mathtt{0x80000000}\) なので、

(*20)
\[
\left(\left(\left(M_{u}x_{i}\right)>>1\right) \text{の下位30ビット}\right)=0
\]

また、 \(\left(M_{l}x_{i+1}\right)>>1\) の項については、 \(M_{l}=\mathtt{0x7fffffff}\) なので、

(*21)
\[
\left(\left(\left(M_{l}x_{i+1}\right)>>1\right) \text{の下位30ビット}\right)=\left(x_{i+1} \text{の第2~31ビット}\right)
\]

さらに、 \(\left(x_{i+1}\cap 1\right)\cdot A\) の項は、 \(x_{i+1}\cap 1\) が \(x_{i+1}\) の最下位ビットを表すことから、次のようになる。

(*22)
\[
\left(\left(\left(x_{i+1}\cap 1\right)\cdot A\right) \text{の下位30ビット}\right)=\left(\left(\left(x_{i+1} \text{の最下位ビット}\right)\cdot A\right) \text{の下位30ビット}\right)
\]

これに式(*12)を代入すると、式(*23)になる。

(*12)再掲
\[
\left(x_{i+1} \text{の最下位ビット}\right)=\left(x_{i+N}\oplus x_{i+M}\right)>>31
\]

(*23)
\[
\left(\left(\left(x_{i+1}\cap 1\right)\cdot A\right) \text{の下位30ビット}\right)=\left(\left(\left(\left(x_{i+N}\oplus x_{i+M}\right)>>31\right)\cdot A\right) \text{の下位30ビット}\right)
\]

以上より、式(*7)の左辺の下位30ビットは、式(*24)のようになる。

(*7)再掲
\[
x_{i+N}\oplus x_{i+M}\oplus\left(\left(M_{u}x_{i}\right)>>1\right)\oplus\left(\left(M_{l}x_{i+1}\right)>>1\right)\oplus\left(\left(x_{i+1}\cap 1\right)\cdot A\right)=0
\]

(*24)
\begin{eqnarray*}
&& \left(\text{式(*7)の左辺の下位30ビット}\right)\\
&=& \left(\left(x_{i+N}\oplus x_{i+M}\oplus\left(\left(M_{u}x_{i}\right)>>1\right)\oplus\left(\left(M_{l}x_{i+1}\right)>>1\right)\oplus\left(\left(x_{i+1}\cap 1\right)\cdot A\right)\right)\text{の下位30ビット}\right)\\
&=& \left(\left(x_{i+N}\oplus x_{i+M}\right)\text{の下位30ビット}\right)\oplus\left(x_{i+1} \text{の第2~31ビット}\right)\\
&&\;\oplus \left(\left(\left(\left(x_{i+N}\oplus x_{i+M}\right)>>31\right)\cdot A\right) \text{の下位30ビット}\right)\\
&=& \left(\left(x_{i+N}\oplus x_{i+M}\oplus\left(\left(\left(x_{i+N}\oplus x_{i+M}\right)>>31\right)\cdot A\right)\right)\text{の下位30ビット}\right)\\
&&\;\oplus\left(x_{i+1} \text{の第2~31ビット}\right)
\end{eqnarray*}

これが式(*7)の右辺の下位30ビット、すなわち \(0\) に等しいので、

(*25)
\begin{eqnarray*}
\left(\left(x_{i+N}\oplus x_{i+M}\oplus\left(\left(\left(x_{i+N}\oplus x_{i+M}\right)>>31\right)\cdot A\right)\right)\text{の下位30ビット}\right)\qquad&&\\
\oplus\left(x_{i+1} \text{の第2~31ビット}\right)&=&0
\end{eqnarray*}

移項すると、

(*26)
\begin{eqnarray*}
&&\left(x_{i+1} \text{の第2~31ビット}\right)\\
&&\qquad=\left(\left(x_{i+N}\oplus x_{i+M}\oplus\left(\left(\left(x_{i+N}\oplus x_{i+M}\right)>>31\right)\cdot A\right)\right)\text{の下位30ビット}\right)
\end{eqnarray*}

よって次式が成り立つ。

(*27)
\begin{eqnarray*}
&&\left(x_{i+1} \text{の第2~31ビット}\right)<<1\\
&&\qquad=\left(\left(x_{i+N}\oplus x_{i+M}\oplus\left(\left(\left(x_{i+N}\oplus x_{i+M}\right)>>31\right)\cdot A\right)\right)\text{の下位30ビット}\right)<<1\\
&&\qquad=M_{l}\left(\left(x_{i+N}\oplus x_{i+M}\oplus\left(\left(\left(x_{i+N}\oplus x_{i+M}\right)>>31\right)\cdot A\right)\right)<<1\right)
\end{eqnarray*}


以上、式(*12)(*19)(*27)をまとめると、以下のようになる。

(*28)
\(i\ge 0\) なる \(i\) について、

\begin{eqnarray*}
\left(x_{i} \text{の最上位ビット}\right)<<31&=&M_{u}\left(\left(x_{i+N}\oplus x_{i+M}\right)<<1\right)\\
\left(x_{i+1} \text{の第2~31ビット}\right)<<1&=&M_{l}\left(\left(x_{i+N}\oplus x_{i+M}\oplus\left(\left(\left(x_{i+N}\oplus x_{i+M}\right)>>31\right)\cdot A\right)\right)<<1\right)\\
\left(x_{i+1} \text{の最下位ビット}\right)&=&\left(x_{i+N}\oplus x_{i+M}\right)>>31
\end{eqnarray*}

メルセンヌ・ツイスタの逆算

さてこれでようやく準備が整ったので、ようやくメルセンヌ・ツイスタの逆算を考えることができる。
話を簡単にするためにまず、メルセンヌ・ツイスタのみによって生成された無限数列 \(\left\{a_{i}\right\}\;\left(-\infty\le i\le\infty\right)\) を考える。
この数列は、任意の \(i\) について式(*28)が成り立つものとする。

いま、この数列から \(N\) 個の連続する値 \(a_{I+1},\,a_{I+2},\,\dots,\,a_{I+N}\) が与えられたとする。
このとき、これらから \(a_{I}\) を求めてみよう。

条件より式(*28)は \(i=I-1,\,I\) のときも成り立つので、次式が成り立つ。

(*29)
\begin{eqnarray*}
\left(a_{I-1} \text{の最上位ビット}\right)<<31&=&M_{u}\left(\left(a_{I+N-1}\oplus a_{I+M-1}\right)<<1\right)\\
\left(a_{I} \text{の第2~31ビット}\right)<<1&=&M_{l}\left(\left(a_{I+N-1}\oplus a_{I+M-1}\oplus\left(\left(\left(a_{I+N-1}\oplus a_{I+M-1}\right)>>31\right)\cdot A\right)\right)<<1\right)\\
\left(a_{I} \text{の最下位ビット}\right)&=&\left(a_{I+N-1}\oplus a_{I+M-1}\right)>>31\\
\\
\left(a_{I} \text{の最上位ビット}\right)<<31&=&M_{u}\left(\left(a_{I+N}\oplus a_{I+M}\right)<<1\right)\\
\left(a_{I+1} \text{の第2~31ビット}\right)<<1&=&M_{l}\left(\left(a_{I+N}\oplus a_{I+M}\oplus\left(\left(\left(a_{I+N}\oplus a_{I+M}\right)>>31\right)\cdot A\right)\right)<<1\right)\\
\left(a_{I+1} \text{の最下位ビット}\right)&=&\left(a_{I+N}\oplus a_{I+M}\right)>>31
\end{eqnarray*}

この式の第1式は \(a_{I-1}\) を求める際に使う式なので、無視する。
また、式(*29)の第5,6式については、条件より \(\left\{a_{i}\right\}\) がメルセンヌ・ツイスタによって生成された数列ゆえに、この2式は常に成り立つ。
よって第5,6式も無視してよい。

【補足】
正確に言えば、式(*29)の第5,6式は「 \(a_{I+1},\,a_{I+2},\,\dots,\,a_{I+N}\) から \(a_{I}\) を逆算できる」ための必要十分条件である。

式(*29)の第2,3,4式から、 \(a_{I}\) が求められる。

(*30)
\begin{eqnarray*}
\left(a_{I} \text{の最上位ビット}\right)<<31&=&M_{u}\left(\left(a_{I+N}\oplus a_{I+M}\right)<<1\right)\\
\left(a_{I} \text{の第2~31ビット}\right)<<1&=&M_{l}\left(\left(a_{I+N-1}\oplus a_{I+M-1}\oplus\left(\left(\left(a_{I+N-1}\oplus a_{I+M-1}\right)>>31\right)\cdot A\right)\right)<<1\right)\\
\left(a_{I} \text{の最下位ビット}\right)&=&\left(a_{I+N-1}\oplus a_{I+M-1}\right)>>31
\end{eqnarray*}

よって \(a_{I}\) は次式で表される。

(*31)
\begin{eqnarray*}
a_{I}
&=&\left(\left(a_{I} \text{の最上位ビット}\right)<<31\right)\oplus\left(\left(a_{I} \text{の第2~31ビット}\right)<<1\right)\oplus\left(a_{I} \text{の最下位ビット}\right)\\
&=&M_{u}\left(\left(a_{I+N}\oplus a_{I+M}\right)<<1\right)\\
&\oplus&M_{l}\left(\left(a_{I+N-1}\oplus a_{I+M-1}\oplus\left(\left(\left(a_{I+N-1}\oplus a_{I+M-1}\right)>>31\right)\cdot A\right)\right)<<1\right)\\
&\oplus&\left(a_{I+N-1}\oplus a_{I+M-1}\right)>>31
\end{eqnarray*}


さて、式(*31)は \(\left\{a_{i}\right\}\) が無限数列の場合に成り立つ式だった。
もしも \(\left\{a_{i}\right\}\) が、初期値を \(a_{0},\,a_{1},\,\dots,\,a_{N-1}\) とする半無限数列だった場合は、式(*31)は \(I\ge 0\) かつ \(I-1\ge 0\) のとき、すなわち \(I\ge 1\) のときに成り立つ。
つまり、 \(N\) 個の初期値のうち \(a_{1},\,a_{2},\,\dots,\,a_{N-1}\) の\(N-1\) 個は求められるが、 \(a_{0}\) だけは求められない。
正確に言うと \(a_{0}\) は最上位ビットのみ求めることができて、次式で表される。

(*32)
\[
\left(a_{0} \text{の最上位ビット}\right)=\left(M_{l}\left(a_{N}\oplus a_{M}\right)\right)>>30
\]

まとめ

メルセンヌ・ツイスタの計算式

\(x_{i},\,x_{i+1},\,\dots,\,x_{i+N-1}\) の \(N\) 個の乱数値について、これらから次の乱数値 \(x_{i+N}\) を得る式は、次式で与えられる。

(*A)
\[
x_{i+N}=x_{i+M}\oplus\left(\left(M_{u}x_{i}\right)>>1\right)\oplus\left(\left(M_{l}x_{i+1}\right)>>1\right)\oplus\left(\left(x_{i+1}\cap 1\right)\cdot A\right)
\]

メルセンヌ・ツイスタの逆算式

\(x_{i+1},\,x_{i+2},\,\dots,\,x_{i+N}\) の \(N\) 個の乱数値について、これらから1つ前の乱数値 \(x_{i}\) を得る式は、次式で与えられる。

(*B)
\begin{eqnarray*}
x_{i}
&=&M_{u}\left(\left(x_{i+N}\oplus x_{i+M}\right)<<1\right)\\
&\oplus&M_{l}\left(\left(x_{i+N-1}\oplus x_{i+M-1}\oplus\left(\left(\left(x_{i+N-1}\oplus x_{i+M-1}\right)>>31\right)\cdot A\right)\right)<<1\right)\\
&\oplus&\left(x_{i+N-1}\oplus x_{i+M-1}\right)>>31
\end{eqnarray*}

逆算可能判定式

\(x_{i+1},\,x_{i+2},\,\dots,\,x_{i+N}\) の \(N\) 個の乱数値について、これらの値から式(*B)によって \(x_{i}\) が正しく逆算できるための必要十分条件は、次式で与えられる。

(*C)
\begin{eqnarray*}
\left(x_{i+1} \text{の第2~31ビット}\right)<<1&=&M_{l}\left(\left(x_{i+N}\oplus x_{i+M}\oplus\left(\left(\left(x_{i+N}\oplus x_{i+M}\right)>>31\right)\cdot A\right)\right)<<1\right)\\
\left(x_{i+1} \text{の最下位ビット}\right)&=&\left(x_{i+N}\oplus x_{i+M}\right)>>31
\end{eqnarray*}

この条件が満たされない場合、 \(x_{i}\) は最上位ビットしか求められない。
なお このとき、 \(x_{i}\) の最上位ビットは次式で与えられる。

(*D)
\[
\left(x_{i} \text{の最上位ビット}\right)=\left(M_{l}\left(x_{i+N}\oplus x_{i+M}\right)\right)>>30
\]

残りの下位31ビットは不定であり、すなわち自由度31である。
とすると解の存在条件となる式が31本あるはずであるが、それが式(*C)である。
※式(*C)は2本しかないように見えるが、第1式が30ビット分の情報量を持っているので、実際には31本の方程式(条件式)である。

実装してみた

以上を踏まえて、正逆双方向に計算可能なメルセンヌ・ツイスタを作ってみた。
速度より汎用性重視(キリッ なので高速化できる人は高速化してやってくださいな……

正逆双方向に計算できるメルセンヌ・ツイスタ - Plus Le Blog

ポケモンGen5でのメルセンヌ・ツイスタ

周知の通り[要出典]メルセンヌ・ツイスタはポケモンGen5でいわゆる「個体値乱数(列)」として使用されている。

ポケモンGen5におけるMTの初期値は、以下の式で与えられる。

(*33)
\[
x_{i}\overset{\mathrm{def}}{=}
\begin{cases}
\mathrm{seed} & \left(i=0 \text{のとき}\right)\\
\left(x_{i-1}\oplus\left(x_{i-1}>>30\right)\right)\cdot\mathtt{0x6c078965}+i & \left(1\le i< N\text{のとき}\right)\\
\end{cases}
\]

ポケモン乱数調整界においては、この「seed」すなわち \(x_{0}\) を得る手順の確立が求められている。
しかし先述の通り、MTの逆算によって確実に求められるのは \(x_{1}\) までであり、 \(x_{0}\) は最上位ビットしか求められない。
そのため \(x_{0}\) は別の方法で \(x_{1}\) から逆算する必要があり、それは実際に可能である。


式(*33)において、 \(i=1\) のとき、

(*34)
\[
x_{1}=\left(x_{0}\oplus\left(x_{0}>>30\right)\right)\cdot\mathtt{0x6c078965}+1
\]

右辺の \(+1\) を移項して、

(*35)
\[
x_{1}-1=\left(x_{0}\oplus\left(x_{0}>>30\right)\right)\cdot\mathtt{0x6c078965}
\]

よって \(Y\overset{\mathrm{def}}{=}x_{1}-1,\,X\overset{\mathrm{def}}{=}x_{0}\oplus\left(x_{0}>>30\right),\,c\overset{\mathrm{def}}{=}\mathtt{0x6c078965}\) とおくと、

(*36)
\[
Y=c\cdot X
\]

普通の数学なら \(X=\frac{Y}{c}\) で \(X\) を求められるのだが、この式の \(Y\) は上位ビットが切り捨てられているため除算によって求めることはできない。
しかしこの形の式は、 \(c\) が奇数のとき \(X\) を求めるアルゴリズムが存在する。
以下にそれを示す。

まず \(c\) が奇数なので、 \(X\) の第1ビット(最下位ビット) \(X\left[1\right]\) は、 \(Y\) の第1ビット \(Y\left[1\right]\) に等しい。

(*37)
\[
X\left[1\right]=Y\left[1\right]
\]

いま \(X'\overset{\mathrm{def}}{=}X-X\left[1\right]\) とすると、

(*38)
\begin{eqnarray*}
Y&=&c\cdot X\\
&=&c\cdot\left(X'+X\left[1\right]\right)\\
&=&c\cdot X'+c\cdot X\left[1\right]\\
&=&c\cdot X'+c\cdot Y\left[1\right]\\
\Longleftrightarrow Y-c\cdot Y\left[1\right]&=&c\cdot X'
\end{eqnarray*}

よって \(Y'\overset{\mathrm{def}}{=}Y-c\cdot Y\left[1\right]\) とすると、

(*39)
\[
Y'=c\cdot X'
\]

定義より \(X',\,Y'\) の第1ビットは \(0\) なので、式(*39)において、第1ビットから第2ビットへの繰り上がりは発生しない。
このことと \(c\) が奇数であることから、 \(X'\) の第2ビット \(X'\left[2\right]\) は、 \(Y'\) の第2ビット \(Y'\left[2\right]\) に等しい。

(*40)
\[
X'\left[2\right]=Y'\left[2\right]
\]

いま \(X''\overset{\mathrm{def}}{=}X'-\left(X'\left[2\right]<<1\right)\) とすると、

(*41)
\begin{eqnarray*}
Y'&=&c\cdot X'\\
&=&c\cdot\left(X''+\left(X'\left[2\right]<<1\right)\right)\\
&=&c\cdot X''+c\cdot\left(X'\left[2\right]<<1\right)\\
&=&c\cdot X''+c\cdot\left(Y'\left[2\right]<<1\right)\\
\Longleftrightarrow Y'-c\cdot\left(Y'\left[2\right]<<1\right)&=&c\cdot X''
\end{eqnarray*}

よって \(Y''\overset{\mathrm{def}}{=}Y'-c\cdot\left(Y'\left[2\right]<<1\right)\) とすると、

(*42)
\[
Y''=c\cdot X''
\]

定義より \(X'',\,Y''\) の第1,2ビットは \(0\) なので、式(*42)において、第2ビット以下から第3ビットへの繰り上がりは発生しない。
このことと \(c\) が奇数であることから、 \(X''\) の第3ビット \(X''\left[3\right]\) は、 \(Y''\) の第3ビット \(Y''\left[3\right]\) に等しい。

(*43)
\[
X''\left[3\right]=Y''\left[3\right]
\]

いま \(X'''\overset{\mathrm{def}}{=}X''-\left(X''\left[3\right]<<2\right)\) とすると、

(*44)
\begin{eqnarray*}
Y''&=&c\cdot X''\\
&=&c\cdot\left(X'''+\left(X''\left[3\right]<<2\right)\right)\\
&=&c\cdot X'''+c\cdot\left(X''\left[3\right]<<2\right)\\
&=&c\cdot X'''+c\cdot\left(Y''\left[3\right]<<2\right)\\
\Longleftrightarrow Y''-c\cdot\left(Y''\left[3\right]<<2\right)&=&c\cdot X'''
\end{eqnarray*}

よって \(Y'''\overset{\mathrm{def}}{=}Y''-c\cdot\left(Y''\left[3\right]<<2\right)\) とすると、

(*45)
\[
Y'''=c\cdot X'''
\]

定義より \(X''',\,Y'''\) の第1~3ビットは \(0\) である。

以下同様に繰り返し計算していけば左辺はいずれ \(0\) になるので、このアルゴリズムを終了することができる。
そしてそれまでに得られた \(X\left[1\right],\,X'\left[2\right],\,X''\left[3\right],\,\dots\) から \(X\) を計算することができる。

さてこの方法により \(X\) が求められたとする。
\(X\) は \(X=x_{0}\oplus\left(x_{0}>>30\right)\) だったので、 tempering の逆関数と同様にして、 \(x_{0}\) は次式で求められる。

(*46)
\[
x_{0}=X\oplus\left(X>>30\right)
\]

この \(x_{0}\) が、求めたかった seed である。

実装してみた

以上の流れをプログラムで書くと、次のようになる。

// C#
uint CalcSeed(uint x1) {
	const uint c = 0x6c078965u;
	const uint index = 1u;

	uint y = x1 - index;

	// 既知の c, y について、 y = c * x を満たす x を求める
	// (※ c が奇数なのでこの方法が使える)
	uint x = 0u;
	uint multiplier = c;
	for( uint mask=1u;mask!=0u;mask<<=1,multiplier<<=1 ) {
		if( ( y & mask ) != 0u ) {
			x |= mask;
			y -= multiplier;
		}
	}

	uint x0 = x ^ ( x >> 30 );
	return x0;
}