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

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

はじめに

↓これの逆関数の3行目の定数値がどこから出てきたのか気になったので調べてみた。


準備

\(N\) bitのブール代数を考える。以下で扱う文字(定数・変数)はすべて整数とする。
定数 \(a\left(\ge 0\right),\,n\) および変数 \(x\left(\ge 0\right)\) について、関数 \(S,\,g,\,f\) を次のように定義する。

(*1)
\begin{eqnarray*}
S\left(x\right)&\overset{\mathrm{def}}{=}&
\begin{cases}
x >> n & \left(n>0 \text{のとき}\right)\\
x & \left(n=0 \text{のとき}\right)\\
x << \left|n\right| & \left(n<0 \text{のとき}\right)
\end{cases}\\
g\left(x\right)&\overset{\mathrm{def}}{=}&a\,S\left(x\right)\\
f\left(x\right)&\overset{\mathrm{def}}{=}&x\oplus g(x)\\
&=&x\oplus a\,S\left(x\right)
\end{eqnarray*}

(※省略された積記号は 算術積\(\cdot\) ではなくすべて論理積 \(\cap\) である。算術積を使う場合は、省略せずに常に \(x\cdot y\) のように積記号を明示するものとする。以降も同様。)
(※関数 \(S\) は、\(0\) や負数の場合を考慮した右ビットシフトである。)


このとき、関数 \(f,\,g,\,S\) には以下の関係が成り立つ。
簡単に導けるので証明は省略する(意訳:めんどい)。

ビットシフトの基本的な性質

論理積,論理和,排他的論理和のビットシフトは、先にビットシフトしたもの同士の論理積,論理和,排他的論理和に等しい。

(*2)
\begin{eqnarray*}
S\left(x\cap y\right) &=&S\left(x\right)\cap S\left(y\right)\\
S\left(x\cup y\right) &=&S\left(x\right)\cup S\left(y\right)\\
S\left(x\oplus y\right) &=&S\left(x\right)\oplus S\left(y\right)
\end{eqnarray*}

ビットシフト関数 \(S\) の合成

\(k\left(\ge 0\right)\) のとき、次式が成り立つ。

(*3)
\[
S^{k}\left(x\right)=
\begin{cases}
x >> \left(k\cdot n\right) & \left(k\cdot n>0 \text{のとき}\right)\\
x & \left(k\cdot n=0 \text{のとき}\right)\\
x << \left|k\cdot n\right| & \left(k\cdot n<0 \text{のとき}\right)
\end{cases}
\]

なお式(*3)は、 \(k=0\) のとき \(S^{k}\left(x\right)\) が恒等関数になることをも表している。

(*4)
\[
S^{0}\left(x\right)=x
\]

関数 \(g\) の合成

\(k\left(\ge 1\right)\) のとき、次式が成り立つ。

(*5)
\[
g^{k}\left(x\right)=S^{k}\left(x\right)\bigcap_{i=0}^{k-1}S^{i}\left(a\right)
\]

特に \(a=\bar{0}\) のときは \(g\left(x\right)=S\left(x\right)\) なので、次式が成り立つ。
(※ \(\bar{x}\) は \(x\) の補数を表す。よって \(\bar{0}\) は全ビットが \(1\) である数であり、例えば \(N=32\) の場合 \(\bar{0}=\mathtt{0xffffffff}\) である。)

(*6)
\[
g^{k}\left(x\right)=S^{k}\left(x\right)
\]

※式(*5)のように \(\bigcap\) を使って書くと分かりづらいが、 \(k=0\) のときは \(g^{0}\left(x\right)=x\) であり、恒等関数である。

(*7)
\[
g^{0}\left(x\right)=x
\]

for文感覚で「 \(\alpha>\beta\) のとき \(\displaystyle\bigcap_{i=\alpha}^{\beta}\mathrm{hoge}\overset{\mathrm{def}}{=}\bar{0}\) 」のように読めば、式(*5)は \(k=0\) のときも成り立つ(Q.なんで \(\bar{0}\) ? A.論理積 \(\cap\) の単位元だから。)。

関数 \(f\) の逆関数 その1

\(n=0\) のとき、 \(f\left(x\right)=x\oplus a\,S\left(x\right)=x\oplus a\,x=\bar{a}\,x\) である。
これは( \(a=0\) の場合を除いて)単射でないので、 \(f\) の逆関数は存在しない。
(※ \(a=0\) の場合 \(f\) は恒等関数 \(f\left(x\right)=x\) なので、逆関数は \(f^{-1}\left(x\right)=x\) である。)

【補足】 \(x\oplus a\,x=\bar{a}\,x\) について
以下の、2つの法則と1つの定義を利用している。

\begin{eqnarray*}
x\oplus a\,x
&=&\bar{0}\,x\oplus a\,x\\
&=&\left(\bar{0}\oplus a\right)\,x\\
&=&\bar{a}\,x
\end{eqnarray*}


\(n\ne 0\) のとき、 \(f\) の逆関数は次式で表される。

(*8)
\begin{eqnarray*}
f^{-1}\left(x\right)
&=&\bigoplus_{k=0}^{t}g^{k}\left(x\right)\\
&=&\bigoplus_{k=0}^{t}S^{k}\left(x\right)\bigcap_{i=0}^{k-1}S^{i}\left(a\right)
\end{eqnarray*}


ただし、 \(t\overset{\mathrm{def}}{=}\max\left(\left\lceil \frac{N}{\left|n\right|}\right\rceil-1,0\right)\) である。

特に \(a=\bar{0}\) のときは、(*6)より次式が成り立つ。

(*9)
\[
f^{-1}\left(x\right)=\bigoplus_{k=0}^{t}S^{k}\left(x\right)
\]

(*6)再掲
\(a=\bar{0}\) のとき、次式が成り立つ。
\[
g^{k}\left(x\right)=S^{k}\left(x\right)
\]


逆関数を求める例題を解いてみよう。

(*10)例題1
【問】 \(N=32\) のとき、次の関数の逆関数を求めよ。
\begin{eqnarray*}
f_{1}\left(x\right)&=&x\oplus \left(x>>11\right)\\
f_{2}\left(x\right)&=&x\oplus \left(x<<7\right)\cap \mathtt{0x9d2c5680}\\
f_{3}\left(x\right)&=&x\oplus \left(x<<15\right)\cap \mathtt{0xefc60000}\\
f_{4}\left(x\right)&=&x\oplus \left(x>>18\right)
\end{eqnarray*}


【解】
順番が違うが、 \(f_{4}\left(x\right)=x\oplus \left(x>>18\right)\) の逆関数から求める。
これは \(a=\bar{0}=\mathtt{0xffffffff}\) の場合であり、1番簡単だからである。

\(n=18\) なので、\(t=\left\lceil \frac{N}{\left|n\right|}\right\rceil-1=\left\lceil \frac{32}{\left|18\right|}\right\rceil-1=1\) より、

(*11) \(f_{4}\) の逆関数
\begin{eqnarray*}
f_{4}^{-1}\left(x\right)
&=&\bigoplus_{k=0}^{1}S^{k}\left(x\right)\\
&=&S^{0}\left(x\right)\oplus S^{1}\left(x\right)\\
&=&x\oplus S\left(x\right)\\
&=&x\oplus \left(x>>18\right)
\end{eqnarray*}


次に、 \(f_{1}\left(x\right)=x\oplus \left(x>>11\right)\) の逆関数を求める。
これも \(a=\bar{0}=\mathtt{0xffffffff}\) の場合なので簡単。

\(n=11\) なので、\(t=\left\lceil \frac{N}{\left|n\right|}\right\rceil-1=\left\lceil \frac{32}{\left|11\right|}\right\rceil-1=2\) より、

(*12) \(f_{1}\) の逆関数
\begin{eqnarray*}
f_{1}^{-1}\left(x\right)
&=&\bigoplus_{k=0}^{2}S^{k}\left(x\right)\\
&=&S^{0}\left(x\right)\oplus S^{1}\left(x\right)\oplus S^{2}\left(x\right)\\
&=&x\oplus S\left(x\right)\oplus S^{2}\left(x\right)\\
&=&x\oplus \left(x>>11\right)\oplus \left(x>>22\right)
\end{eqnarray*}

式(*11)と式(*12)を見比べて気づいたかもしれないが、 \(t\) が大きくなるほど計算量が増えて面倒になる。
つまり、元の関数のビットシフト幅 \(\left|n\right|\) が大きいほうが、逆関数を求めるのは簡単になる。


次に、ビット幅が大きい \(f_{3}\left(x\right)=x\oplus \left(x<<15\right)\cap \mathtt{0xefc60000}\) の逆関数を求める。
左ビットシフトなので \(n\) は負になることに注意。

\(n=-15\) なので、\(t=\left\lceil \frac{N}{\left|n\right|}\right\rceil-1=\left\lceil \frac{32}{\left|-15\right|}\right\rceil-1=2\) より、

(*13) \(f_{3}\) の逆関数
\begin{eqnarray*}
f_{3}^{-1}\left(x\right)
&=&\bigoplus_{k=0}^{2}S^{k}\left(x\right)\bigcap_{i=0}^{k-1}S^{i}\left(\mathtt{0xefc60000}\right)\\
&=&S^{0}\left(x\right)\oplus S^{0}\left(\mathtt{0xefc60000}\right)S^{1}\left(x\right)\oplus S^{0}\left(\mathtt{0xefc60000}\right)S^{1}\left(\mathtt{0xefc60000}\right)S^{2}\left(x\right)\\
&=&x\oplus\mathtt{0xefc60000}\cap S\left(x\right)\oplus\mathtt{0xefc60000}\cap\mathtt{0x00000000}\cap S^{2}\left(x\right)\\
&=&x\oplus\mathtt{0xefc60000}\cap S\left(x\right)\\
&=&x\oplus\left(x<<15\right)\cap\mathtt{0xefc60000}
\end{eqnarray*}


最後に、 \(f_{2}\left(x\right)=x\oplus \left(x<<7\right)\cap \mathtt{0x9d2c5680}\) の逆関数を求める。

\(n=-7\) なので、\(t=\left\lceil \frac{N}{\left|n\right|}\right\rceil-1=\left\lceil \frac{32}{\left|-7\right|}\right\rceil-1=4\) より、

(*14) \(f_{2}\) の逆関数
\begin{eqnarray*}
f_{2}^{-1}\left(x\right)
&=&\bigoplus_{k=0}^{4}S^{k}\left(x\right)\bigcap_{i=0}^{k-1}S^{i}\left(\mathtt{0x9d2c5680}\right)\\
&=&S^{0}\left(x\right)\\
&\oplus&S^{0}\left(\mathtt{0x9d2c5680}\right)S^{1}\left(x\right)\\
&\oplus&S^{0}\left(\mathtt{0x9d2c5680}\right)S^{1}\left(\mathtt{0x9d2c5680}\right)S^{2}\left(x\right)\\
&\oplus&S^{0}\left(\mathtt{0x9d2c5680}\right)S^{1}\left(\mathtt{0x9d2c5680}\right)S^{2}\left(\mathtt{0x9d2c5680}\right)S^{3}\left(x\right)\\
&\oplus&S^{0}\left(\mathtt{0x9d2c5680}\right)S^{1}\left(\mathtt{0x9d2c5680}\right)S^{2}\left(\mathtt{0x9d2c5680}\right)S^{3}\left(\mathtt{0x9d2c5680}\right)S^{4}\left(x\right)\\
&=&x\\
&\oplus&\mathtt{0x9d2c5680}\cap S\left(x\right)\\
&\oplus&\mathtt{0x9d2c5680}\cap\mathtt{0x962b4000}\cap S^{2}\left(x\right)\\
&\oplus&\mathtt{0x9d2c5680}\cap\mathtt{0x962b4000}\cap\mathtt{0x15a00000}\cap S^{3}\left(x\right)\\
&\oplus&\mathtt{0x9d2c5680}\cap\mathtt{0x962b4000}\cap\mathtt{0x15a00000}\cap\mathtt{0xd0000000}\cap S^{4}\left(x\right)\\
&=&x\\
&\oplus&\mathtt{0x9d2c5680}\cap S\left(x\right)\\
&\oplus&\mathtt{0x94284000}\cap S^{2}\left(x\right)\\
&\oplus&\mathtt{0x14200000}\cap S^{3}\left(x\right)\\
&\oplus&\mathtt{0x10000000}\cap S^{4}\left(x\right)\\
&=&x\\
&\oplus&\left(x<<7\right)\cap\mathtt{0x9d2c5680}\\
&\oplus&\left(x<<14\right)\cap\mathtt{0x94284000}\\
&\oplus&\left(x<<21\right)\cap\mathtt{0x14200000}\\
&\oplus&\left(x<<28\right)\cap\mathtt{0x10000000}
\end{eqnarray*}


以上より式(*11)~(*14)をまとめると、次のようになる。

(*10)【問】例題1(再掲)
\begin{eqnarray*}
f_{1}\left(x\right)&=&x\oplus \left(x>>11\right)\\
f_{2}\left(x\right)&=&x\oplus \left(x<<7\right)\cap \mathtt{0x9d2c5680}\\
f_{3}\left(x\right)&=&x\oplus \left(x<<15\right)\cap \mathtt{0xefc60000}\\
f_{4}\left(x\right)&=&x\oplus \left(x>>18\right)
\end{eqnarray*}

(*15)【解】例題1
\begin{eqnarray*}
f_{1}^{-1}\left(x\right)&=&x\oplus \left(x>>11\right)\oplus \left(x>>22\right)\\
f_{2}^{-1}\left(x\right)
&=&x\oplus\left(x<<7\right)\cap\mathtt{0x9d2c5680}\oplus\left(x<<14\right)\cap\mathtt{0x94284000}\\
&&\oplus\left(x<<21\right)\cap\mathtt{0x14200000}\oplus\left(x<<28\right)\cap\mathtt{0x10000000}\\
f_{3}^{-1}\left(x\right)&=&x\oplus\left(x<<15\right)\cap\mathtt{0xefc60000}\\
f_{4}^{-1}\left(x\right)&=&x\oplus \left(x>>18\right)
\end{eqnarray*}

関数 \(f\) の合成

\(k\left(\ge 1\right)\) のとき、次式が成り立つ。

(*16)
\begin{eqnarray*}
f^{k}\left(x\right)
&=&\bigoplus_{i=0}^{k}\left(\left({}_{k}C_{i}\bmod{2}\right)\cdot g^{i}\left(x\right)\right)\\
&=&\bigoplus_{i=0}^{k}\left(\left({}_{k}C_{i}\bmod{2}\right)\cdot S^{i}\left(x\right)\bigcap_{j=0}^{i-1}S^{j}\left(a\right)\right)
\end{eqnarray*}

特に \(k=2^{m}\;\left(0\le m\in\mathbb{Z}\right)\) の場合、次のようになる。

(*17)
\begin{eqnarray*}
f^{k}\left(x\right)
&=&x\oplus g^{k}\left(x\right)\\
&=&x\oplus S^{k}\left(x\right)\bigcap_{j=0}^{k-1}S^{j}\left(a\right)
\end{eqnarray*}

関数 \(f\) の逆関数 その2

\(n=0\) のとき \(f\) は逆関数を持たないので、 \(n\ne 0\) の場合を考える。
まず、数列 \(\left\{x_{i}\right\}\) を次のように定義する。

(*18)
\[
\left\{
\begin{array}{lcll}
x_{0} & \overset{\mathrm{def}}{=} & x & \\
x_{i+1} & \overset{\mathrm{def}}{=} & f^{2^{i}}\left(x_{i}\right) & \\
& = & \displaystyle x_{i}\oplus S^{2^{i}}\left(x_{i}\right)\bigcap_{j=0}^{2^{i}-1}S^{j}\left(a\right) & \left(i\ge 0 \text{のとき}\right)
\end{array}
\right.
\]

このとき、 \(f\) の逆関数は次式で表せる。

(*19)
\[
f^{-1}\left(x\right)=x_{t}
\]


ただし、 \(t\overset{\mathrm{def}}{=}\max\left(\left\lceil\log_{2}\frac{N}{\left|n\right|}\right\rceil,0\right)\) である。


逆関数の求め方その1は逆関数を1つの式で表現できるが、その2は漸化式でしか表現できない。
しかし繰り返し計算の回数 \(t\) を比較すると、その1が \(t=\max\left(\left\lceil \frac{N}{\left|n\right|}\right\rceil-1,0\right)\) だったのに対し、その2は \(t=\max\left(\left\lceil\log_{2}\frac{N}{\left|n\right|}\right\rceil,0\right)\) である。
これらは、ビットシフト幅 \(\left|n\right|\) を小さくしていったとき、その2の \(t\) よりその1の \(t\) の方が増加が早いということを示している。
言い換えると、その2はその1より計算量が少ないということだ。

【追記】

(*16)再掲
\[
f^{k}\left(x\right)=\bigoplus_{i=0}^{k}\left(\left({}_{k}C_{i}\bmod{2}\right)\cdot g^{i}\left(x\right)\right)
\]

式(*16)において、 \(k=2^{m}-1\;\left(0\le m\in\mathbb{Z}\right)\) のとき \({}_{k}C_{i}\bmod{2}=1\;\left(0\le i\le k\right)\) なので、次式が成り立つ。

\[
f^{k}\left(x\right)=\bigoplus_{i=0}^{k}g^{i}\left(x\right)
\]

これは式(*8)の右辺と一致する。

(*8)再掲
\[
f^{-1}\left(x\right)=\bigoplus_{k=0}^{t}g^{k}\left(x\right)
\]


ただし、 \(t=\max\left(\left\lceil \frac{N}{\left|n\right|}\right\rceil-1,0\right)\) である。

このことと、式(*8)の \(t\) について \(i>t\) を満たすとき \(g^{i}\left(x\right)=0\) となることから、 \(f\) の逆関数は次式で表される。

\[
f^{-1}\left(x\right)=f^{2^{m}-1}\left(x\right)
\]


ただし、 \(m\) は、 \(2^{m}-1\ge\max\left(\left\lceil \frac{N}{\left|n\right|}\right\rceil-1,0\right)\) を満たす最小の \(m\) である。
言い換えると、 \(m=\max\left(\left\lceil\log_{2}\frac{N}{\left|n\right|}\right\rceil,0\right)\) である。


……「その2は漸化式でしか表現できない(キリッ」とか言ったけど、普通に1本の式で書けました。
というか \(f\) どうしの合成だから普通に代入できるの見落としてました。
もっとも、計算量減らす工夫としては漸化式の形にする必要があったので無駄ではなかった……ということにしておく。


実際に確認してみよう。
例題1の中でビットシフト幅が1番小さい、関数 \(f_{2}\) の逆関数を求めてみる。
その1と異なり、その2の方法は逆関数を1つの関数の形で書けないので、「逆関数を求める漸化式」を求めることになる。

(*20)例題2
【問】 \(N=32\) のとき、次の関数の逆関数を求めるプログラムを示せ。
\[
f_{2}\left(x\right)=x\oplus \left(x<<7\right)\cap \mathtt{0x9d2c5680}
\]


【解】
\(n=-7\) なので、 \(t=\max\left(\left\lceil\log_{2}\frac{32}{\left\lceil -7\right\rceil}\right\rceil,0\right)=3\) より、

(*21)
\begin{eqnarray*}
x_{0} &=&x\\
x_{1} &=&f_{2}\left(x_{0}\right)\\
&=&x_{0}\oplus S^{0}\left(\mathtt{0x9d2c5680}\right)S^{1}\left(x_{0}\right)\\
&=&x_{0}\oplus\mathtt{0x9d2c5680}\cap S\left(x_{0}\right)\\
&=&x_{0}\oplus\left(x_{0}<<7\right)\cap\mathtt{0x9d2c5680}\\
x_{2} &=&f_{2}^{2}\left(x_{1}\right)\\
&=&x_{1}\oplus S^{0}\left(\mathtt{0x9d2c5680}\right)S^{1}\left(\mathtt{0x9d2c5680}\right)S^{2}\left(x_{1}\right)\\
&=&x_{1}\oplus\mathtt{0x9d2c5680}\cap\mathtt{0x962b4000}\cap S^{2}\left(x_{1}\right)\\
&=&x_{1}\oplus\left(x_{1}<<14\right)\cap\mathtt{0x94284000}\\
x_{3} &=&f_{2}^{4}\left(x_{2}\right)\\
&=&x_{2}\oplus S^{0}\left(\mathtt{0x9d2c5680}\right)S^{1}\left(\mathtt{0x9d2c5680}\right)S^{2}\left(\mathtt{0x9d2c5680}\right)S^{3}\left(\mathtt{0x9d2c5680}\right)S^{4}\left(x_{2}\right)\\
&=&x_{2}\oplus\mathtt{0x9d2c5680}\cap\mathtt{0x962b4000}\cap\mathtt{0x15a00000}\cap\mathtt{0xd0000000}\cap S^{4}\left(x_{2}\right)\\
&=&x_{2}\oplus\left(x_{2}<<28\right)\cap\mathtt{0x10000000}\\
f_{2}^{-1}\left(x\right)&=&x_{3}
\end{eqnarray*}

まとめると、

(*22)
\begin{eqnarray*}
x_{0} &=&x\\
x_{1} &=&x_{0}\oplus\left(x_{0}<<7\right)\cap\mathtt{0x9d2c5680}\\
x_{2} &=&x_{1}\oplus\left(x_{1}<<14\right)\cap\mathtt{0x94284000}\\
x_{3} &=&x_{2}\oplus\left(x_{2}<<28\right)\cap\mathtt{0x10000000}\\
f_{2}^{-1}\left(x\right)&=&x_{3}
\end{eqnarray*}

よって次のようなプログラムで逆関数を求めることができる。

// C#
// その2の方法で求めた、 f_2(x) の逆関数
uint f2_inv(uint x) {
	x ^= ( x <<  7 ) & 0x9d2c5680u;
	x ^= ( x << 14 ) & 0x94284000u;
	x ^= ( x << 28 ) & 0x10000000u;
	return x;
}

(*15)参考:その1で求めた逆関数(再掲)
\begin{eqnarray*}
f_{2}^{-1}\left(x\right)
&=&x\\
&\oplus&\left(x<<7\right)\cap\mathtt{0x9d2c5680}\\
&\oplus&\left(x<<14\right)\cap\mathtt{0x94284000}\\
&\oplus&\left(x<<21\right)\cap\mathtt{0x14200000}\\
&\oplus&\left(x<<28\right)\cap\mathtt{0x10000000}
\end{eqnarray*}

メルセンヌ・ツイスタの tempering の逆関数

以上を踏まえると、メルセンヌ・ツイスタの tempering 部分の逆関数は次のように書ける。

// C#
// メルセンヌ・ツイスタの tempering 部分
uint tempering(uint x) {
	x ^= ( x >> 11 );		// f_1(x)
	x ^= ( x <<  7 ) & 0x9d2c5680u;	// f_2(x)
	x ^= ( x << 15 ) & 0xefc60000u;	// f_3(x)
	x ^= ( x >> 18 );		// f_4(x)
	return x;
}
// その2の方法で作った、 tempering の逆関数
uint tempering_inv_plus(uint x) {
	x ^= ( x >> 18 );		// f_4(x) の逆関数
	x ^= ( x << 15 ) & 0xefc60000u;	// f_3(x) の逆関数
	x ^= ( x <<  7 ) & 0x9d2c5680u;	// f_2(x) の逆関数 ここから
	x ^= ( x << 14 ) & 0x94284000u;	// 
	x ^= ( x << 28 ) & 0x10000000u;	// ここまで
	x ^= ( x >> 11 ) ^ ( x >> 22 );	// f_1(x) の逆関数

	// ※その2の方法で f_1(x) の逆関数を求めた場合は、次のようになる。
//	x ^= ( x >> 11 );		// f_1(x) の逆関数 ここから
//	x ^= ( x >> 22 );		// ここまで
	return x;
}
// 参考:冒頭のoupoさんの提示した逆関数(その1の方法で作った、 tempering の逆関数)
uint tempering_inv_oupo(uint x) {
	x ^=   ( x >> 18 );			// f_4(x) の逆関数
	x ^=   ( x << 15 ) & 0xefc60000u;	// f_3(x) の逆関数
	x ^= ( ( x <<  7 ) & 0x9d2c5680u )	// f_2(x) の逆関数 ここから
	   ^ ( ( x << 14 ) & 0x94284000u )	// 
	   ^ ( ( x << 21 ) & 0x14200000u )	// 
	   ^ ( ( x << 28 ) & 0x10000000u );	// ここまで
	x ^=   ( x >> 11 ) ^ ( x >> 22 );	// f_1(x) の逆関数
	return x;
}

oupoさんの式と少し違うけど、uint型の全範囲の入力に対して計算結果が同じになることは確認済み。
ちなみに、uint型の全範囲の入力に対して逆関数を適用した値を求めるだけのプログラムの実行時間を計測したところ、 tempering_inv_plus は20936ミリ秒、 tempering_inv_oupo は22698ミリ秒であり、有意差は認められなかった。
ビッシシフト幅 \(\left|n\right|\) がもっと小さくないと、その2がその1より劇的に有利になるわけではないようだ。

【追記】有意差あった。というかその2よりその1の方が速かった。詳しくは後述の追記2を参照。

おわりに

いずれも厳密な証明はしてないので注意。
もしかしたら \(a,\,n,\,N\) の組み合わせによって例外があるかもしれないので、必要に応じて各自確認してから使用することを推奨します。

ちなみにoupoさんは行列を使って冒頭の逆関数を求めたようだが、普通に(ブール)代数的に解けた。
……図解の1つでも描いた方がいいのかなぁ。

【追記】おまけ

// C#
// メルセンヌ・ツイスタの tempering 部分
uint tempering(uint x) {
	x ^= ( x >> 11 );
	x ^= ( x <<  7 ) & 0x9d2c5680u;
	x ^= ( x << 15 ) & 0xefc60000u;
	x ^= ( x >> 18 );
	return x;
}
// tempering の逆関数
uint tempering_inv(uint x) {
	x ^= ( x >> 18 );
	x ^= ( x << 15 ) & 0xefc60000u;
	x ^= ( x <<  7 ) & 0x9d2c5680u;
	x ^= ( x << 14 ) & 0x94284000u;
	x ^= ( x << 28 ) & 0x10000000u;
	x ^= ( x >> 11 );
	x ^= ( x >> 22 );
	return x;
}

tempering の逆関数が上記の tempering_inv であることは既に示したとおり。
では、 tempering_inv の逆関数は本当に元の tempering に一致するのだろうか?
確認してみた。


例題1で使った関数 \(f_{1}\left(x\right)\) ~ \(f_{4}\left(x\right)\) について、 \(g_{i}\left(x\right)\overset{\mathrm{def}}{=}x\oplus f_{i}\left(x\right)\) と定義する。

(*10)再掲
\begin{eqnarray*}
f_{1}\left(x\right)&=&x\oplus \left(x>>11\right)\\
f_{2}\left(x\right)&=&x\oplus \left(x<<7\right)\cap \mathtt{0x9d2c5680}\\
f_{3}\left(x\right)&=&x\oplus \left(x<<15\right)\cap \mathtt{0xefc60000}\\
f_{4}\left(x\right)&=&x\oplus \left(x>>18\right)
\end{eqnarray*}

(*23)
\begin{eqnarray*}
g_{1}\left(x\right)&=&\left(x>>11\right)\\
g_{2}\left(x\right)&=&\left(x<<7\right)\cap \mathtt{0x9d2c5680}\\
g_{3}\left(x\right)&=&\left(x<<15\right)\cap \mathtt{0xefc60000}\\
g_{4}\left(x\right)&=&\left(x>>18\right)
\end{eqnarray*}

式(*10)(*23)を使うと、 tempering , tempering_inv は次のように書ける(長いのでそれぞれ \(T,\,T_{inv}\) と表記する)。

(*24)tempering
\[
\begin{array}{rclclcl}
x_{0}&=&x\\
x_{1}&=&f_{1}\left(x_{0}\right)&=&x_{0}\oplus g_{1}\left(x_{0}\right)&=&x_{0}\oplus \left(x_{0}>>11\right)\\
x_{2}&=&f_{2}\left(x_{1}\right)&=&x_{1}\oplus g_{2}\left(x_{1}\right)&=&x_{1}\oplus \left(x_{1}<<7\right)\cap \mathtt{0x9d2c5680}\\
x_{3}&=&f_{3}\left(x_{2}\right)&=&x_{2}\oplus g_{3}\left(x_{2}\right)&=&x_{2}\oplus \left(x_{2}<<15\right)\cap \mathtt{0xefc60000}\\
x_{4}&=&f_{4}\left(x_{3}\right)&=&x_{3}\oplus g_{4}\left(x_{3}\right)&=&x_{3}\oplus \left(x_{3}>>18\right)\\
T\left(x\right)&=&x_{4}
\end{array}
\]

(*25)tempering_inv
\[
\begin{array}{rclclcl}
y_{0}&=&y\\
y_{1}&=&f_{4}\left(y_{0}\right)&=&y_{0}\oplus g_{4}\left(y_{0}\right)&=&y_{0}\oplus\left(y_{0}>>18\right)\\
y_{2}&=&f_{3}\left(y_{1}\right)&=&y_{1}\oplus g_{3}\left(y_{1}\right)&=&y_{1}\oplus\left(y_{1}<<15\right)\cap\mathtt{0xefc60000}\\
y_{3}&=&f_{2}\left(y_{2}\right)&=&y_{2}\oplus g_{2}\left(y_{2}\right)&=&y_{2}\oplus\left(y_{2}<<7\right)\cap\mathtt{0x9d2c5680}\\
y_{4}&=&f_{2}^{2}\left(y_{3}\right)&=&y_{3}\oplus g_{2}^{2}\left(y_{3}\right)&=&y_{3}\oplus\left(y_{3}<<14\right)\cap\mathtt{0x94284000}\\
y_{5}&=&f_{2}^{4}\left(y_{4}\right)&=&y_{4}\oplus g_{2}^{4}\left(y_{4}\right)&=&y_{4}\oplus\left(y_{4}<<28\right)\cap\mathtt{0x10000000}\\
y_{6}&=&f_{1}\left(y_{5}\right)&=&y_{5}\oplus g_{1}\left(y_{5}\right)&=&y_{5}\oplus\left(y_{5}>>11\right)\\
y_{7}&=&f_{1}^{2}\left(y_{6}\right)&=&y_{6}\oplus g_{1}^{2}\left(y_{6}\right)&=&y_{6}\oplus\left(y_{6}>>22\right)\\
T_{inv}\left(y\right)&=&y_{7}
\end{array}
\]

(※式(*25)では、式(*17)を使ったことに注意。)

(*17)再掲
\(k=2^{m}\;\left(0\le m\in\mathbb{Z}\right)\) のとき、次式が成り立つ。
\[
f^{k}\left(x\right)=x\oplus g^{k}\left(x\right)
\]


さていま、その1の方法を使って tempering_inv の逆関数を求めてみる。
例題とたいして変わらないので計算過程を省略すると、次のようになる。

(*26)
\[
\begin{array}{rclcl}
z_{0}&=&z\\
z_{1}&=&y_{7}^{-1}\left(z_{0}\right)&=&z_{0}\oplus g_{1}^{2}\left(z_{0}\right)\\
z_{2}&=&y_{6}^{-1}\left(z_{1}\right)&=&z_{1}\oplus g_{1}\left(z_{1}\right)\oplus g_{1}^{2}\left(z_{1}\right)\\
z_{3}&=&y_{5}^{-1}\left(z_{2}\right)&=&z_{2}\oplus g_{2}^{4}\left(z_{2}\right)\\
z_{4}&=&y_{4}^{-1}\left(z_{3}\right)&=&z_{3}\oplus g_{2}^{2}\left(z_{3}\right)\oplus g_{2}^{4}\left(z_{3}\right)\\
z_{5}&=&y_{3}^{-1}\left(z_{4}\right)&=&z_{4}\oplus g_{2}\left(z_{4}\right)\oplus g_{2}^{2}\left(z_{4}\right)\oplus g_{2}^{3}\left(z_{4}\right)\oplus g_{2}^{4}\left(z_{4}\right)\\
z_{6}&=&y_{2}^{-1}\left(z_{5}\right)&=&z_{5}\oplus g_{3}\left(z_{5}\right)\quad\left(\,=f_{3}\left(z_{5}\right)\,\right)\\
z_{7}&=&y_{1}^{-1}\left(z_{6}\right)&=&z_{6}\oplus g_{4}\left(z_{6}\right)\quad\left(\,=f_{4}\left(z_{6}\right)\,\right)\\
T_{inv}^{-1}\left(z\right)&=&z_{7}
\end{array}
\]

(※式中で \(y_{i}\) を関数のように扱っているが、その意味はおおよそ見当がつくと思うので厳密な定義は省略する。訳:面倒。)

ここで \(z_{2}\) に \(z_{1}\) を代入すると、

(*27)
\begin{eqnarray*}
z_{1}&=&z_{0}\oplus g_{1}^{2}\left(z_{0}\right)\\
z_{2} &=&z_{1}\oplus g_{1}\left(z_{1}\right)\oplus g_{1}^{2}\left(z_{1}\right)\\
&=&\left(z_{0}\oplus g_{1}^{2}\left(z_{0}\right)\right)\oplus g_{1}\left(z_{0}\oplus g_{1}^{2}\left(z_{0}\right)\right)\oplus g_{1}^{2}\left(z_{0}\oplus g_{1}^{2}\left(z_{0}\right)\right)\\
&=&z_{0}\oplus g_{1}^{2}\left(z_{0}\right)\oplus g_{1}\left(z_{0}\right)\oplus g_{1}^{3}\left(z_{0}\right)\oplus g_{1}^{2}\left(z_{0}\right)\oplus g_{1}^{4}\left(z_{0}\right)\\
&=&z_{0}\oplus g_{1}\left(z_{0}\right)\oplus g_{1}^{3}\left(z_{0}\right)\oplus g_{1}^{4}\left(z_{0}\right)\\
&=&z_{0}\oplus g_{1}\left(z_{0}\right)\\
&=&f_{1}\left(z_{0}\right)
\end{eqnarray*}

同様に、 \(z_{4}\) に \(z_{3}\) を代入したものを \(z_{5}\) に代入すると、

(*28)
\begin{eqnarray*}
z_{3} &=&z_{2}\oplus g_{2}^{4}\left(z_{2}\right)\\
z_{4} &=&z_{3}\oplus g_{2}^{2}\left(z_{3}\right)\oplus g_{2}^{4}\left(z_{3}\right)\\
&=&\left(z_{2}\oplus g_{2}^{4}\left(z_{2}\right)\right)\oplus g_{2}^{2}\left(z_{2}\oplus g_{2}^{4}\left(z_{2}\right)\right)\oplus g_{2}^{4}\left(z_{2}\oplus g_{2}^{4}\left(z_{2}\right)\right)\\
&=&z_{2}\oplus g_{2}^{4}\left(z_{2}\right)\oplus g_{2}^{2}\left(z_{2}\right)\oplus g_{2}^{6}\left(z_{2}\right)\oplus g_{2}^{4}\left(z_{2}\right)\oplus g_{2}^{8}\left(z_{2}\right)\\
&=&z_{2}\oplus g_{2}^{2}\left(z_{2}\right)\oplus g_{2}^{6}\left(z_{2}\right)\oplus g_{2}^{8}\left(z_{2}\right)\\
&=&z_{2}\oplus g_{2}^{2}\left(z_{2}\right)\\
z_{5} &=&z_{4}\oplus g_{2}\left(z_{4}\right)\oplus g_{2}^{2}\left(z_{4}\right)\oplus g_{2}^{3}\left(z_{4}\right)\oplus g_{2}^{4}\left(z_{4}\right)\\
&=&\left(z_{2}\oplus g_{2}^{2}\left(z_{2}\right)\right)\oplus g_{2}\left(z_{2}\oplus g_{2}^{2}\left(z_{2}\right)\right)\oplus g_{2}^{2}\left(z_{2}\oplus g_{2}^{2}\left(z_{2}\right)\right)\\
&&\;\oplus g_{2}^{3}\left(z_{2}\oplus g_{2}^{2}\left(z_{2}\right)\right)\oplus g_{2}^{4}\left(z_{2}\oplus g_{2}^{2}\left(z_{2}\right)\right)\\
&=&z_{2}\oplus g_{2}^{2}\left(z_{2}\right)\oplus g_{2}\left(z_{2}\right)\oplus g_{2}^{3}\left(z_{2}\right)\oplus g_{2}^{2}\left(z_{2}\right)\oplus g_{2}^{4}\left(z_{2}\right)\\
&&\;\oplus g_{2}^{3}\left(z_{2}\right)\oplus g_{2}^{5}\left(z_{2}\right)\oplus g_{2}^{4}\left(z_{2}\right)\oplus g_{2}^{6}\left(z_{2}\right)\\
&=&z_{2}\oplus g_{2}\left(z_{2}\right)\oplus g_{2}^{5}\left(z_{2}\right)\oplus g_{2}^{6}\left(z_{2}\right)\\
&=&z_{2}\oplus g_{2}\left(z_{2}\right)\\
&=&f_{2}\left(z_{2}\right)
\end{eqnarray*}


式(*26)~(*28)をまとめると、次のようになる。

(*29)
\begin{eqnarray*}
z_{0}&=&z\\
z_{2}&=&f_{1}\left(z_{0}\right)\\
z_{5}&=&f_{2}\left(z_{2}\right)\\
z_{6}&=&f_{3}\left(z_{5}\right)\\
z_{7}&=&f_{4}\left(z_{6}\right)\\
T_{inv}^{-1}\left(z\right)&=&z_{7}
\end{eqnarray*}

これは tempering の式(*24)と一致しているので、 tempering_inv の逆関数は \(T_{inv}=T\) であることが示された。

(*24)tempering(再掲)
\begin{eqnarray*}
x_{0}&=&x\\
x_{1}&=&f_{1}\left(x_{0}\right)\\
x_{2}&=&f_{2}\left(x_{1}\right)\\
x_{3}&=&f_{3}\left(x_{2}\right)\\
x_{4}&=&f_{4}\left(x_{3}\right)\\
T\left(x\right)&=&x_{4}
\end{eqnarray*}

【追記】おまけ2

// C#
static void Main(string[] args) {
	int t1 = 0;
	int t2 = 0;
	int max = 10;
	uint x = 0u;

	t1 = Environment.TickCount;
	for( int i=0;i<max;i++ ) {
		do {
			// f_2(x) その1
			x ^=  ( ( x <<  7 ) & 0x9d2c5680u )
				^ ( ( x << 14 ) & 0x94284000u )
				^ ( ( x << 21 ) & 0x14200000u )
				^ ( ( x << 28 ) & 0x10000000u );
		} while( ++x != 0u );
	}
	t2 = Environment.TickCount;
	Console.WriteLine("その1 : {0}ミリ秒",t2-t1);

	t1 = Environment.TickCount;
	for( int i=0;i<max;i++ ) {
		do {
			// f_2(x) その2
			x ^= ( x <<  7 ) & 0x9d2c5680u;
			x ^= ( x << 14 ) & 0x94284000u;
			x ^= ( x << 28 ) & 0x10000000u;
		} while( ++x != 0u );
	}
	t2 = Environment.TickCount;
	Console.WriteLine("その2 : {0}ミリ秒",t2-t1);

	//対照実験
	t1 = Environment.TickCount;
	for( int i=0;i<max;i++ ) {
		do {
			// NOP
		} while( ++x != 0u );
	}
	t2 = Environment.TickCount;
	Console.WriteLine("NOP    : {0}ミリ秒",t2-t1);

	return;
}

上のプログラムでその1とその2の計算速度を測ったところ、実行結果は次のようになった。

その1 : 142116ミリ秒
その2 : 166843ミリ秒
NOP    : 25756ミリ秒
続行するには何かキーを押してください . . .

\(\frac{166843}{142116}=1.174\) すなわち約17%の速度差。
ループ部分の25756ミリ秒を双方から除いて考えるなら、\(\frac{166843-25756}{142116-25756}=\frac{141087}{116360}=1.213\) すなわち約21%の速度差となる。
この差を有意差でないとみなすのはさすがに無理があるだろう。

\(f_{2}\) の逆関数の計算で使う演算の回数はそれぞれ、

これを見ると、「ビットシフト1回+論理積1回+排他的論理和1回」より「代入2回」の方が時間がかかっていることが分かる。
代入は重い処理なのだろうか?

ビットシフト,論理積排他的論理和に各1クロック、代入に3クロック必要と仮定すると、その1が計15クロック、その2が計18クロック必要なので、 \(\frac{18}{15}=1.2\) となり前述の速度比 1.213 とほぼ同じ値になる。
よって各演算の必要クロック数は上の値であると推測できる。

この予想が正しいとすれば、ビットシフト幅 \(\left|n\right|\) が大きい場合は代入の3クロックがボトルネックになるので、その1の方法を使うのがよいということになる。
一方 \(\left|n\right|\) が小さい場合、その2の方法ならビットシフト,論理積排他的論理和の回数を対数的回数に抑えられるので、その2の方法を使うのがよさそうだ。
\(\left|n\right|\) が中くらいの値の場合? 知らんがな(´・ω・`)

【追記】おまけ3

参考画像描いたけど疲れたので解説は後日(書くとは言ってない)。

【追記】おまけ4

通常のメルセンヌ・ツイスタでは \(n\) や \(a\) が定数なので \(\displaystyle\bigcap_{i=0}^{k-1}S^{i}\left(a\right)\) の値をあらかじめ求めておける。
パラメータを動的生成するメルセンヌ・ツイスタ(※pdf)では今回の実装は使えないが、代わりにXorshift の逆関数に関する考察/正逆双方向に計算できる Xorshift - Plus Le Blogのようにwhileループを使って実装すればいい。

// C#
private readonly int _u;	// 11
private readonly int _s;	//  7
private readonly int _t;	// 15
private readonly int _l;	// 18
private readonly uint _b;	// 0x9d2c5680u
private readonly uint _c;	// 0xefc60000u

private static uint f_R(uint x,int n) {
	x ^= ( x >> n );
	return x;
}
private static uint f_L(uint x,int n,uint a) {
	x ^= ( x << n ) & a;
	return x;
}
private static uint f_R_inv(uint x,int n) {
	while( n < N ) {
		x ^= ( x >> n );
		n <<= 1;
	}
	return x;
}
private static uint f_L_inv(uint x,int n,uint a) {
	while( n < N && a != 0u ) {
		x ^= ( x << n ) & a;
		a &= ( a << n );
		n <<= 1;
	}
	return x;
}

// f_R, f_L を使ったメルセンヌ・ツイスタの tempering 部分の実装例
uint tempering(uint x) {
	x = f_R(x,_u);		// f_1(x)
	x = f_L(x,_s,_b);	// f_2(x)
	x = f_L(x,_t,_c);	// f_3(x)
	x = f_R(x,_l);		// f_4(x)
	return x;
}

// f_R, f_L を使った tempering の逆関数の実装例
uint tempering_inv_plus(uint x) {
	x = f_R_inv(x,_l);		// f_4(x) の逆関数
	x = f_L_inv(x,_t,_c);	// f_3(x) の逆関数
	x = f_L_inv(x,_s,_b);	// f_2(x) の逆関数
	x = f_R_inv(x,_u);		// f_1(x) の逆関数
	return x;
}