Xorshift の逆関数に関する考察／正逆双方向に計算できる Xorshift
はじめに
取説（というかライセンス表記）のないSwitchでポケモン新作が出たら乱数生成器や初期化方法を特定するのがめんどくさくなりそうだから、あらかじめいろんな疑似乱数の逆関数を求めておけばいいんじゃね？ な企画第１弾。
ということでとりあえず目についた“Xorshift”の逆関数を求めてみる。なお、第２弾以降があるかは時間と気分次第。たぶんないだろうけど。
準備
\(N\) bitのブール代数を考える。以下で扱う文字（定数・変数）はすべて整数とする。
定数（パラメータ変数） \(a\left(\ge 0\right),\,m\left(\ge 0\right),\,n\left(>0\right)\) および変数 \(x\left(\ge 0\right)\) について、関数 \(S_{R},\,S_{L},\,h_{R},\,h_{L},\,g_{R},\,g_{L},\,f_{R},\,f_{L}\) を次のように定義する。
（＊１）
\begin{eqnarray*}
S_{R}\left(x\right)&\stackrel{\mathrm{def}}{=}&x >> n\\
S_{L}\left(x\right)&\stackrel{\mathrm{def}}{=}&x << n\\
h_{R}\left(x;n,a,m\right)&=&S_{R}^{m}\left(x\right)\bigcap_{i=0}^{m-1}S_{R}^{i}\left(a\right)\\
h_{L}\left(x;n,a,m\right)&=&S_{L}^{m}\left(x\right)\bigcap_{i=0}^{m-1}S_{L}^{i}\left(a\right)\\
g_{R}\left(x;n,a\right)
&\stackrel{\mathrm{def}}{=}&h_{R}\left(x;n,a,1\right)\\
&=&a\,S_{R}\left(x\right)\\
g_{L}\left(x;n,a\right)
&\stackrel{\mathrm{def}}{=}&h_{L}\left(x;n,a,1\right)\\
&=&a\,S_{L}\left(x\right)\\
f_{R}\left(x;n,a\right)
&\stackrel{\mathrm{def}}{=}&x\oplus g_{R}\left(x;n,a\right)\\
&=&x\oplus a\,S_{R}\left(x\right)\\
f_{L}\left(x;n,a\right)
&\stackrel{\mathrm{def}}{=}&x\oplus g_{L}\left(x;n,a\right)\\
&=&x\oplus a\,S_{L}\left(x\right)
\end{eqnarray*}
（※省略された積記号は 算術積\(\cdot\) ではなくすべて論理積 \(\cap\) である。算術積を使う場合は、省略せずに常に \(x\cdot y\) のように積記号を明示するものとする。以降も同様。）
なお、 \(a,\,m,\,n\) はパラメータ変数であって変数でないので、 \(h_{R},\,h_{L},\,g_{R},\,g_{L},\,f_{R},\,f_{L}\) はいずれも四変数関数や三変数関数でなく一変数関数である。
そのため \(a,\,m,\,n\) が文脈などから明らかな場合は、これらのうち一部または全部を関数の引数から省略して書くことにする。
また 以下しばらくは、式中の関数の添字の \(R,\,L\) を入れ替えても成り立つので、略記することにする（ただしその場合、 \(R\) または \(L\) は、一連の数式中ではすべて同じものを使わなければならないことに注意）。
したがって、以降では（＊１）の定義代わりに次式（＊２）を使う。
（＊２）
\begin{eqnarray*}
S\left(x\right)&\stackrel{\mathrm{def}}{=}&
\begin{cases}
x >> n & \left(R \text{のとき}\right)\\
x << n & \left(L \text{のとき}\right)
\end{cases}\\
h\left(x\right)&=&S^{m}\left(x\right)\bigcap_{i=0}^{m-1}S^{i}\left(a\right)\\
g\left(x\right)&\stackrel{\mathrm{def}}{=}&a\,S\left(x\right)\\
f\left(x\right)&\stackrel{\mathrm{def}}{=}&x\oplus a\,S\left(x\right)
\end{eqnarray*}
上記の関数 \(S,\,g,\,f\) の定義は、メルセンヌ・ツイスタのtemperingの逆関数に関する考察 - Plus Le Blogと同じものである（ただし \(n>0\) の制限がついたため関数 \(S\) の見た目が若干変わっている）。
（上記リンクは以降でもたびたび使うので、以降は「★」と記す。）
\(h\) は今回初出の関数だが、★の（＊５）より次の関係が成り立つ（これが \(h\) の本来の定義と考えてもよい）。
★（＊５）再掲
\[
g^{k}\left(x\right)=S^{k}\left(x\right)\bigcap_{i=0}^{k-1}S^{i}\left(a\right)
\]
（＊３）
\begin{eqnarray*}
h\left(x;n,a,m\right)
&=&S^{m}\left(x\right)\bigcap_{i=0}^{m-1}S^{i}\left(a\right)\\
&=&g^{m}\left(x;n,a\right)
\end{eqnarray*}
このとき★の（＊８）より、関数 \(f\) には次の逆関数が存在する（もちろん、変数 \(x\) についての逆関数である）。
★（＊８）再掲
\[
f^{-1}\left(x\right)=\bigoplus_{k=0}^{t}g^{k}\left(x\right)
\]
ただし、 \(t\stackrel{\mathrm{def}}{=}\max\left(\left\lceil \frac{N}{\left|n\right|}\right\rceil-1,0\right)\) である。
（＊４）
\begin{eqnarray*}
f^{-1}\left(x;n,a\right)
&=&\bigoplus_{k=0}^{t}g^{k}\left(x;n,a\right)\\
&=&\bigoplus_{k=0}^{t}h\left(x;n,a,k\right)\\
&=&\bigoplus_{k=0}^{t}S^{k}\left(x\right)\bigcap_{i=0}^{k-1}S^{i}\left(a\right)
\end{eqnarray*}
ただし、 \(t\stackrel{\mathrm{def}}{=}\max\left(\left\lceil \frac{N}{\left|n\right|}\right\rceil-1,0\right)\) である。
また ★の（＊18）（＊19）を使えば、上記（＊４）より少ない演算回数で \(f\) の逆関数を求められる。
★（＊18）再掲
\[
\left\{
\begin{array}{lcll}
x_{0} & \stackrel{\mathrm{def}}{=} & x & \\
x_{i+1} & \stackrel{\mathrm{def}}{=} & \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.
\]★（＊19）再掲
\[
f^{-1}\left(x\right)=x_{t}
\]
ただし、 \(t\stackrel{\mathrm{def}}{=}\max\left(\left\lceil\log_{2}\frac{N}{\left|n\right|}\right\rceil,0\right)\) である。
（＊５）
数列 \(\left\{x_{i}\right\}\) を次のように定義する。\[
\left\{
\begin{array}{lcll}
x_{0} & \stackrel{\mathrm{def}}{=} & x & \\
x_{i+1} & \stackrel{\mathrm{def}}{=} & x_{i}\oplus h\left(x_{i};n,a,2^{i}\right) & \left(i\ge 0 \text{のとき}\right)
\end{array}
\right.
\]
このとき、 \(f\) の逆関数は次式で表せる。
（＊６）
\[
f^{-1}\left(x\right)=x_{T}
\]
ただし、 \(T\stackrel{\mathrm{def}}{=}\max\left(\left\lceil\log_{2}\frac{N}{\left|n\right|}\right\rceil,0\right)\) である。
さらに、 \(a=\bar{0}\) ならば逆関数はより簡単に書ける。
（＊７）
\[
f\left(x\right)\stackrel{\mathrm{def}}{=}x\oplus S\left(x\right)
\]（＊８）
\begin{eqnarray*}
f^{-1}\left(x\right)
&=&\bigoplus_{k=0}^{t}S^{k}\left(x\right)\\
&=&\bigoplus_{k=0}^{T}S^{2^{k}}\left(x\right)
\end{eqnarray*}
ただし \(t\) および \(T\) は以下の通り。
\begin{eqnarray*}
t&\stackrel{\mathrm{def}}{=}&\max\left(\left\lceil \frac{N}{\left|n\right|}\right\rceil-1,0\right)\\
T&\stackrel{\mathrm{def}}{=}&\max\left(\left\lceil\log_{2}\frac{N}{\left|n\right|}\right\rceil,0\right)
\end{eqnarray*}
以上をプログラムで書くと、以下のようになる。
// C# private const int N = 32; #region （＊１） f_R, f_L public static uint f_R(uint x,int n,uint a) { x ^= ( x >> n ) & a; return x; } public static uint f_L(uint x,int n,uint a) { x ^= ( x << n ) & a; return x; } #endregion #region （＊４） f_R, f_L の逆関数（※この実装は１つ次の実装より演算量が多いのでコメントアウトしている） //public static uint f_R_inv(uint x,int n,uint a) { // uint x0 = x; // uint a0 = a; // int n0 = n; // while( n < N && a != 0u ) { // x ^= ( x0 >> n ) & a; // a &= ( a0 >> n ); // n += n0; // } // return x; //} //public static uint f_L_inv(uint x,int n,uint a) { // uint x0 = x; // uint a0 = a; // int n0 = n; // while( n < N && a != 0u ) { // x ^= ( x0 << n ) & a; // a &= ( a0 << n ); // n += n0; // } // return x; //} #endregion #region （＊５）（＊６） f_R, f_L の逆関数 public static uint f_R_inv(uint x,int n,uint a) { while( n < N && a != 0u ) { x ^= ( x >> n ) & a; a &= ( a >> n ); n <<= 1; } return x; } public 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; } #endregion // ---------------- #region （＊７） a == ~0 の場合の f_R, f_L public static uint f_R(uint x,int n) { x ^= ( x >> n ); return x; } public static uint f_L(uint x,int n) { x ^= ( x << n ); return x; } #endregion #region （＊８） a == ~0 の場合の f_R, f_L の逆関数（上側の式による実装：※この実装は１つ次の実装より演算量が多いのでコメントアウトしている） //public static uint f_R_inv(uint x,int n) { // uint x0 = x; // int n0 = n; // while( n < N ) { // x ^= ( x0 >> n ); // n += n0; // } // return x; //} //public static uint f_L_inv(uint x,int n) { // uint x0 = x; // int n0 = n; // while( n < N ) { // x ^= ( x0 << n ); // n += n0; // } // return x; //} #endregion #region （＊８） a == ~0 の場合の f_R, f_L の逆関数（下側の式による実装） public static uint f_R_inv(uint x,int n) { while( n < N ) { x ^= ( x >> n ); n <<= 1; } return x; } public static uint f_L_inv(uint x,int n) { while( n < N ) { x ^= ( x << n ); n <<= 1; } return x; } #endregion
（※C#以外で実装する際は、言語によっては右ビットシフトが論理シフトでなく算術シフトの場合があるので注意。）
Xorshift
上で定義・実装した \(f_{R},\,f_{L},\,f_{R}^{-1},\,f_{L}^{-1}\) を使って、Wikipediaにある32bit, 64bit, 96bit, 128bit版XorshiftのC言語での実装とその逆関数をC#で書くと次のようになる。
// C# public abstract class XorshiftBase { #region constructor protected XorshiftBase(int a,int b,int c) { this._a = a; this._b = b; this._c = c; } #endregion #region static field private const int N = 32; #endregion #region field protected readonly int _a; protected readonly int _b; protected readonly int _c; #endregion #region static method protected static uint f_R(uint seed,int shift) { seed ^= ( seed >> shift ); return seed; } protected static uint f_L(uint seed,int shift) { seed ^= ( seed << shift ); return seed; } protected static uint f_R_inv(uint seed,int shift) { while( shift < N ) { seed ^= ( seed >> shift ); shift <<= 1; } return seed; } protected static uint f_L_inv(uint seed,int shift) { while( shift < N ) { seed ^= ( seed << shift ); shift <<= 1; } return seed; } #endregion #region abstract method public abstract uint Next(); public abstract uint Prev(); #endregion } public class Xorshift32 : XorshiftBase { // constructor public Xorshift32() : base(13,17,5) { } public Xorshift32(int a,int b,int c) : base(a,b,c) { } // filed private uint _x = 2463534242u; // override method public override sealed uint Next() { _x = f_L(_x,_a); _x = f_R(_x,_b); _x = f_L(_x,_c); return _x; } public override sealed uint Prev() { _x = f_L_inv(_x,_c); _x = f_R_inv(_x,_b); _x = f_L_inv(_x,_a); return _x; } } public class Xorshift64 : XorshiftBase { // constructor public Xorshift64() : base(13,7,17) { } public Xorshift64(int a,int b,int c) : base(a,b,c) { } // static field private const int N = 64; // field private ulong _x = 88172645463325252ul; // static method protected static ulong f_R(ulong seed,int shift) { seed ^= ( seed >> shift ); return seed; } protected static ulong f_L(ulong seed,int shift) { seed ^= ( seed << shift ); return seed; } protected static ulong f_R_inv(ulong seed,int shift) { while( shift < N ) { seed ^= ( seed >> shift ); shift <<= 1; } return seed; } protected static ulong f_L_inv(ulong seed,int shift) { while( shift < N ) { seed ^= ( seed << shift ); shift <<= 1; } return seed; } // override method public override sealed uint Next() { _x = f_L(_x,_a); _x = f_R(_x,_b); _x = f_L(_x,_c); return (uint)_x; } public override sealed uint Prev() { _x = f_L_inv(_x,_c); _x = f_R_inv(_x,_b); _x = f_L_inv(_x,_a); return (uint)_x; } } public class Xorshift96 : XorshiftBase { // constructor public Xorshift96() : base(3,19,6) { } public Xorshift96(int a,int b,int c) : base(a,b,c) { } // field private uint _x = 123456789u; private uint _y = 362436069u; private uint _z = 521288629u; // override method public override sealed uint Next() { uint t = f_L(_x,_a) ^ f_R(_y,_b) ^ f_L(_z,_c); _x = _y; _y = _z; _z = t; return _z; } public override sealed uint Prev() { uint t = _z; _z = _y; _y = _x; t ^= f_L(_z,_c) ^ f_R(_y,_b); _x = f_L_inv(t,_a); return _z; } } public class Xorshift128 : XorshiftBase { // constructor public Xorshift128() : base(11,8,19) { } public Xorshift128(int a,int b,int c) : base(a,b,c) { } // field private uint _x = 123456789u; private uint _y = 362436069u; private uint _z = 521288629u; private uint _w = 88675123u; // override method public override sealed uint Next() { uint t = f_L(_x,_a); _x = _y; _y = _z; _z = _w; _w = f_R(_w,_c) ^ f_R(t,_b); return _w; } public override sealed uint Prev() { uint t = _w; _w = _z; _z = _y; _y = _x; t ^= f_R(_w,_c); t = f_R_inv(t,_b); _x = f_L_inv(t,_a); return _w; } }
