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

はじめに

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

temperingの逆関数できた pic.twitter.com/rssTofDVgG
— oupo (隠居中) (@_oupo_) 2014年10月22日

準備

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

（＊１）
\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\) には以下の関係が成り立つ。
簡単に導けるので証明は省略する（意訳：めんどい）。

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

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

（＊２）
\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)\) のとき、次式が成り立つ。

（＊３）
\[
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}
\]

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

（＊４）
\[
S^{0}\left(x\right)=x
\]

関数 \(g\) の合成

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

（＊５）
\[
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}\) である。）

（＊６）
\[
g^{k}\left(x\right)=S^{k}\left(x\right)
\]

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

（＊７）
\[
g^{0}\left(x\right)=x
\]

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

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

\(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\) について
以下の、２つの法則と１つの定義を利用している。

法則：論理積と排他的論理和の間には分配法則が成り立つ

定義： \(N\) bit整数の排他的論理和とは、各bitごとの排他的論理和である

法則： \(1\) bitのブール代数において、 \(\bar{0}\oplus x=1\oplus x=\bar{x}\) が成り立つ

\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\) の逆関数は次式で表される。

（＊８）
\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}\) のときは、（＊６）より次式が成り立つ。

（＊９）
\[
f^{-1}\left(x\right)=\bigoplus_{k=0}^{t}S^{k}\left(x\right)
\]

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

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

（＊10）例題１
【問】 \(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}\) の場合であり、１番簡単だからである。

\(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）【問】例題１（再掲）
\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）【解】例題１
\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\) の逆関数　その２

\(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)\) である。

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

【追記】

（＊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)
\]

これは式（＊８）の右辺と一致する。

（＊８）再掲
\[
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)\) である。

このことと、式（＊８）の \(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)\) である。

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

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

（＊20）例題２
【問】 \(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#
// その２の方法で求めた、 f_2(x) の逆関数
uint f2_inv(uint x) {
	x ^= ( x <<  7 ) & 0x9d2c5680u;
	x ^= ( x << 14 ) & 0x94284000u;
	x ^= ( x << 28 ) & 0x10000000u;
	return x;
}

（＊15）参考：その１で求めた逆関数（再掲）
\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;
}
// その２の方法で作った、 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) の逆関数

	// ※その２の方法で f_1(x) の逆関数を求めた場合は、次のようになる。
//	x ^= ( x >> 11 );		// f_1(x) の逆関数　ここから
//	x ^= ( x >> 22 );		// ここまで
	return x;
}
// 参考：冒頭のoupoさんの提示した逆関数（その１の方法で作った、 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|\) がもっと小さくないと、その２がその１より劇的に有利になるわけではないようだ。
【追記】有意差あった。というかその２よりその１の方が速かった。詳しくは後述の追記２を参照。

おわりに

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

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

【追記】おまけ

// 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 に一致するのだろうか？
確認してみた。

例題１で使った関数 \(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)
\]

さていま、その１の方法を使って 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*}

【追記】おまけ２

// 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) その１
			x ^=  ( ( x <<  7 ) & 0x9d2c5680u )
				^ ( ( x << 14 ) & 0x94284000u )
				^ ( ( x << 21 ) & 0x14200000u )
				^ ( ( x << 28 ) & 0x10000000u );
		} while( ++x != 0u );
	}
	t2 = Environment.TickCount;
	Console.WriteLine("その１ : {0}ミリ秒",t2-t1);

	t1 = Environment.TickCount;
	for( int i=0;i<max;i++ ) {
		do {
			// f_2(x) その２
			x ^= ( x <<  7 ) & 0x9d2c5680u;
			x ^= ( x << 14 ) & 0x94284000u;
			x ^= ( x << 28 ) & 0x10000000u;
		} while( ++x != 0u );
	}
	t2 = Environment.TickCount;
	Console.WriteLine("その２ : {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;
}

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

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

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

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

その１：ビットシフト４回，論理積４回，排他的論理和４回，代入１回。
その２：ビットシフト３回，論理積３回，排他的論理和３回，代入３回。

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

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

【追記】おまけ３

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

@plus315 飽きた疲れた。この図でだいたいの流れが伝わるといいな。　その２の図はその１とだいたい同じだし、もうゴールしてもいいよね…… pic.twitter.com/UKWiakepib
— ぷらす＠LINEARITY_CHK=ON (@plus315) 2014年10月26日

【追記】おまけ４

通常のメルセンヌ・ツイスタでは \(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;
}