## イントロ
[pokefinderの実装](https://github.com/Admiral-Fish/PokeFinder/blob/dadea37e8259603dd44e9280bf28bd7c8181fe38/Source/Core/RNG/Xoroshiro.cpp#L61)を見に行ったらなんか短かった。
試してみたら、SIMDを抜きにしても前に書いた記事の実装の倍くらい速かった。

調べてみたら[Xoroshiroのオリジナル実装](https://prng.di.unimi.it/xoroshiro128plus.c)でも同じようにしてジャンプ関数が提供されていた(ただし固定長のみ)。

## 考察しようと思ったモチベーション
2倍程度の速度向上だけであれば後回しにしていたと思う。
しかし、既存の実装のネックであった空間計算量（128bit * 128次配列 * 64）が大きく改善されるので、取り組まないわけにはいかなかった。

## どういう原理なのか

任意のベクトル$x$に対し、$f^0(x), f^1(x), \cdots, f^{127}(x)$は一次独立になる[^1]、すなわち基底になる。つまり、Xoroshiro128pやXorShiftの任意の状態ベクトルは、任意に固定したベクトル$x$について$c_0f^0(x)+c_1 f^1(x)+ \cdots + c_{127}f^{127}(x) \ (c_i \in \{0,1\} )$の形式で表現される。

基準として取ったベクトルを$x_0$とし、$x_n = f^n(x_0)$を表現したときの係数を$\{ c_i \}$とする。この$\{ c_i \}$は$x_0$の取り方には依らず、ただ$n$の値のみによって定まる。これは基準とするベクトルを$x_0^{\prime} = f^k(x_0)$に取り換えたとしても、$f$が線形写像であることにより

$\begin{align}
x_n^{\prime}
&= f^n(x_0^{\prime}) \\
&= f^n(f^k(x_0)) \\
&= f^k(f^n(x_0)) \\
&= f^k(c_0f^0(x_0)+c_1 f^1(x_0)+ \cdots + c_{127}f^{127}(x_0)) \\
&= c_0f^k(x_0)+c_1 f^{k+1}(x_0)+ \cdots + c_{127}f^{k+127}(x_0) \\
&= c_0f^0(x_0^{\prime})+c_1 f^1(x_0^{\prime})+ \cdots + c_{127}f^{127}(x_0^{\prime}) \\
\end{align}$

が成り立つからである。

この$\{c_i\}$を$128=2^7, 2^8, 2^9, \cdots$についてあらかじめ計算しておけば、任意の$n$に対して${\rm log}_2(n)$回の小ジャンプの組み合わせてジャンプが可能になる。

## 係数の求め方
各$2^k$に対応する係数$\{c_i\}$は以下のようにして求められる。

- 適当なベクトル$x$から基底$x_0 = f^0(x),\ x_1 = f(x),\ \cdots,\ x_{127} = f^{127}(x)$を計算する。
- $x_i$を縦ベクトルと見なし、横に並べて$128$次正方行列$T$を作る。
- $x_{2^k}$を横ベクトルと見なし、${c_i}$も横ベクトル$c$と見なせば、$x_{2^k} = cT$の関係が成り立っており、$T$は正則。
- $T^{-1}$を求め、$x_{2^k}T^{-1}$を計算すれば$c = \{c_i\}$が得られる。

### Tips
- $x_{2^{k+1}}$は$x_{2^k}$に対して$2^k$ジャンプを適用させれば得られる。
- $f^{128}$の小ジャンプは普通に128回回したほうが速い。
- $f^{128}$から$f^{256}$までの計算回数も128回なので、$n < 256$であれば線形に進めたほうが速い。
- pokefinderでの実装では$f^{128}$の小ジャンプは線形に進めているが、$f^{128}$から$f^{256}$へのジャンプは線形結合による計算を行ってしまっている。

## 計算するプログラムと計算結果
https://gist.github.com/yatsuna827/fe396453d474ba7727497fae4418a792

## その他
- こういうのこそPythonで書いた方が楽なのかもしれん。
- MTやSFMTでも理論上は同じやり方でジャンプ関数が実装できるけど、次元が大きいので逆行列を求める部分に工夫が必要になるし、何より19937消費以内ならfor文で回すほうが速い。

[^1]:松本眞 (2004)「 [有限体の擬似乱数への応用](http://www.math.sci.hiroshima-u.ac.jp/m-mat/TEACH/0407-2.pdf) 」の定理4.20による。$c_0f^0(x)+c_1 f^1(x)+ \cdots + c_{127}f^{127}(x)$が原始多項式になることは命題4.10を見よ。