線形合同法（LCG）に関する考察
LCGの復習的な。
新情報はたぶんないです。
はじめに
この記事で使われる文字（定数・変数）はすべて整数であり、特に、明記した場合や添字などの場合を除き、 \(M\) 未満の非負整数とする。
また この記事中の数式は、特に明記しない限り、法を \(M\overset{\mathrm{def}}{=}2^{N}\;\left(2\le N\in\mathbb{N}\right)\) とする合同式として扱う。
つまり、単に \(x=y\) と書かれた式は \(x\equiv y\;\pmod{M}\) の意味として扱う。
なお、省略された積記号は通常の算術積 \(\cdot\) である（論理積 \(\cap\) として扱う記事が続いたので念のため）。
\(f\) の定義
定数 \(a,\,b\) について、関数 \(f\) を次のように定義する。
（＊１）
\[
f\left(x\right)\overset{\mathrm{def}}{=}ax+b
\]
ただし \(a,\,b\) は次式の関係を満たすものとする。
（＊２）
\begin{eqnarray*}
a&\equiv& 1\pmod{4}\\
b&\equiv& 1\pmod{2}
\end{eqnarray*}
このとき、 定数 \(a,\,b\) および関数 \(f\) に関して、以下の関係や性質が成り立つ。
\(a,\,b\) の逆元 \(a^{-1},\,b^{-1}\)
式（＊２）より \(a,\,b\) はどちらも \(M\left(=2^{N}\right)\) と互いに素なので、逆元が存在する。
よって以降、 \(a,\,b\) の逆元をそれぞれ \(a^{-1},\,b^{-1}\) とする。
（＊３）
\[
\begin{array}{lclcl}
aa^{-1} &=& a^{-1}a &=& 1\\
bb^{-1} &=& b^{-1}b &=& 1
\end{array}
\]
\(f^{n}\left(x\right)\)
（＊４）
\(n\ge 0\) なる \(n\) について、次式が成り立つ。
\[
f^{n}\left(x\right)=a^{n}x+b\sum_{k=0}^{n-1}a^{k}
\]
【証明】
帰納法で。
特に、 \(n=2\) のときは、
（＊５）
\[
f^{2}\left(x\right)=a^{2}x+b\left(1+a\right)
\]
\(f^{n}\left(0\right)\)
（＊６）
\(n\ge 0\) なる \(n\) について、次式が成り立つ。
\[
f^{n}\left(0\right)=b\sum_{k=0}^{n-1}a^{k}
\]
【証明】
式（＊４）に \(x=0\) を代入。
式（＊４）（＊６）より、 \(n\ge 0\) のとき次式が成り立つ。
（＊７）
\[
f^{n}\left(x\right)=a^{n}x+f^{n}\left(0\right)
\]
\(a^{n}\)
（＊８）
\(n\ge 0\) なる \(n\) について、次式が成り立つ。
\[
a^{n}=\left(a-1\right)b^{-1}f^{n}\left(0\right)+1
\]
【証明】
式（＊６）の両辺に \(\left(a-1\right)b^{-1}\) をかけて整理する。
\begin{eqnarray*}
\left(a-1\right)b^{-1}\cdot f^{n}\left(0\right)
&=& \left(a-1\right)b^{-1}\cdot b\sum_{k=0}^{n-1}a^{k}\\
&=& \left(a-1\right)\sum_{k=0}^{n-1}a^{k}\\
&=& a^{n}-1\\
\Longleftrightarrow\quad
a^{n}
&=& \left(a-1\right)b^{-1}f^{n}\left(0\right)+1
\end{eqnarray*}
なお式（＊８）において、 \(a-1\) は偶数ゆえに逆元が存在しないので、 \(f^{n}\left(0\right)=\cdots\) の形に変形しようなどと考えてはならない。
法 \(M\) と互いに素でない \(a-1\) をかけたときに同値性が崩れているのだ。
\(f^{-1}\left(x\right)\)
（＊９）
\[
f^{-1}\left(x\right)=a^{-1}x-a^{-1}b
\]
【証明】
\(f^{-1}\left(x\right)\overset{\mathrm{def}}{=}cx+d\) とおくと、 \(f^{-1}\left(f\left(x\right)\right)=c\left(ax+b\right)+d=acx+\left(bc+d\right)\) となる。
左辺は \(x\) に等しいので、恒等関係から \(ac=1,\,bc+d=0\) が成り立つ。
変形すると \(c=a^{-1},\,d=-bc=-a^{-1}b\) となり式（＊９）を得る。
\(0\) に \(f\) を何回施すと \(\alpha\) になるか？
ここまでは準備。
ここからが本題。
（★１）【問１】
ある整数 \(\alpha\) が与えられたとき、次式を満たす非負整数 \(n\left(< M\right)\) を求めよ。
\[
\alpha=f^{n}\left(0\right)
\]
【解１】
式（★１）（＊６）より、式（＊10）が成り立つ。
（＊６）再掲
\[
f^{n}\left(0\right)=b\sum_{k=0}^{n-1}a^{k}
\]
（＊10）
\[
\alpha=f^{n}\left(0\right)=b\sum_{k=0}^{n-1}a^{k}
\]
式（＊10）について、 \(n\) の偶奇で場合分けをする。
- \(n=2n'\;\left(n'\ge 0\right)\) のとき
（＊11）
\begin{eqnarray*}
\alpha
&=& b\sum_{k=0}^{n-1}a^{k}\\
&=& b\sum_{k=0}^{2n'-1}a^{k}\\
&=& b\left(0+a^{0}\sum_{k=0}^{2n'-1}a^{k}\right)
\end{eqnarray*}
- \(n=2n'+1\;\left(n'\ge 0\right)\) のとき
（＊12）
\begin{eqnarray*}
\alpha
&=& b\sum_{k=0}^{n-1}a^{k}\\
&=& b\sum_{k=0}^{2n'}a^{k}\\
&=& b\left(a^{0}+\sum_{k=1}^{2n'}a^{k}\right)\\
&=& b\left(1+a^{1}\sum_{k=0}^{2n'-1}a^{k}\right)
\end{eqnarray*}
式（＊11）（＊12）をまとめると、
（＊13）
\[
\alpha=
\begin{cases}
\displaystyle
b\left(0+a^{0}\sum_{k=0}^{2n'-1}a^{k}\right)&\left(n=2n' \text{のとき}\right)\\
\displaystyle
b\left(1+a^{1}\sum_{k=0}^{2n'-1}a^{k}\right)&\left(n=2n'+1 \text{のとき}\right)\\
\end{cases}
\]
よって \(n\) の最下位ビットを \(n_{0}\,\left(\in\left\{0,1\right\}\right)\) として \(n\overset{\mathrm{def}}{=}2n'+n_{0}\;\left(n'\ge 0\right)\) とおくと、次式が成り立つ。
（＊14）
\begin{eqnarray*}
\alpha
&=& b\left(n_{0}+a^{n_{0}}\sum_{k=0}^{2n'-1}a^{k}\right)\\
&=& b\left(n_{0}+a^{n_{0}}\sum_{k=0}^{n'-1}\left(a^{2k}+a^{2k+1}\right)\right)\\
&=& b\left(n_{0}+a^{n_{0}}\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}\right)\\
\Longleftrightarrow\;
\alpha-bn_{0}&=& ba^{n_{0}}\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}
\end{eqnarray*}
（＊２）再掲
\begin{eqnarray*}
a&\equiv& 1\pmod{4}\\
b&\equiv& 1\pmod{2}
\end{eqnarray*}
式（＊２）より \(a\) は奇数なので \(1+a\) は偶数であり、ゆえに式（＊14）の右辺も偶数である。
よって、 \(\alpha\) の偶奇は \(bn_{0}\) の偶奇に等しい。
式（＊２）より \(b\) は奇数なので、 \(\alpha\) の最下位ビット \(\alpha_{0}\) は \(n_{0}\) に等しい。
（＊15）
\[
\alpha_{0}=n_{0}
\]
よって式（＊14）（＊15）より、
（＊16）
\[
\alpha-b\alpha_{0}=ba^{\alpha_{0}}\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}
\]
ここで、 \(\alpha',\,a',\,b'\) を次のように定義する。
（＊17）
\begin{eqnarray*}
\alpha' &\overset{\mathrm{def}}{=}& \alpha-b\alpha_{0}\\
a' &\overset{\mathrm{def}}{=}& a^{2}\\
b' &\overset{\mathrm{def}}{=}& ba^{\alpha_{0}}\left(1+a\right)
\end{eqnarray*}
このとき、 \(\alpha_{0}\) が \(\alpha\) の最下位ビットであることと式（＊２）より、次の関係が成り立つ。
（＊18）
\begin{eqnarray*}
\alpha' &\equiv& 0\pmod{2}\\
a' &\equiv& 1\pmod{8}\\
b' &\equiv& 2\pmod{4}
\end{eqnarray*}
式（＊16）（＊17）より、次式が成り立つ。
（＊19）
\[
\alpha'=b'\sum_{k=0}^{n'-1}\left(a'\right)^{k}
\]
\(n'\) の最下位ビットを \(n'_{0}\,\left(\in\left\{0,1\right\}\right)\) として \(n'\overset{\mathrm{def}}{=}2n''+n'_{0}\;\left(n''\ge 0\right)\) とおくと、次式が成り立つ。
（＊20）
\begin{eqnarray*}
\alpha'
&=& b'\left(n'_{0}+\left(a'\right)^{n'_{0}}\sum_{k=0}^{2n''-1}\left(a'\right)^{k}\right)\\
&=& b'\left(n'_{0}+\left(a'\right)^{n'_{0}}\sum_{k=0}^{n''-1}\left(\left(a'\right)^{2k}+\left(a'\right)^{2k+1}\right)\right)\\
&=& b'\left(n'_{0}+\left(a'\right)^{n'_{0}}\left(1+a'\right)\sum_{k=0}^{n''-1}\left(a'\right)^{2k}\right)\\
\Longleftrightarrow\;
\alpha'-b'n'_{0}&=& b'\left(a'\right)^{n'_{0}}\left(1+a'\right)\sum_{k=0}^{n''-1}\left(a'\right)^{2k}
\end{eqnarray*}
式（＊18）より、式（＊20）の右辺は \(\left(\text{奇数の\(4\)倍}\right)\cdot\sum\left(\text{略}\right)\) の形であり、また左辺は \(\alpha'\) が偶数， \(b'\) が奇数の \(2\) 倍である。
よって、 \(\alpha'\) の第１ビット \(\alpha'_{1}\) は \(n'_{0}\;\left(=n_{1}\right)\) に等しい。
（＊21）
\begin{eqnarray*}
\alpha'_{1}
&=& n'_{0}\\
&=& n_{1}
\end{eqnarray*}
よって式（＊20）（＊21）より、
（＊22）
\[
\alpha'-b'\alpha'_{1}=b'\left(a'\right)^{\alpha'_{1}}\left(1+a'\right)\sum_{k=0}^{n''-1}\left(a'\right)^{2k}
\]
ここで、 \(\alpha'',\,a'',\,b''\) を次のように定義する。
（＊23）
\begin{eqnarray*}
\alpha'' &\overset{\mathrm{def}}{=}& \alpha'-b'\alpha'_{1}\\
a'' &\overset{\mathrm{def}}{=}& \left(a'\right)^{2}\\
b'' &\overset{\mathrm{def}}{=}& b'\left(a'\right)^{\alpha'_{1}}\left(1+a'\right)
\end{eqnarray*}
このとき、 \(\alpha'_{1}\) が \(\alpha'\) の第１ビットであることと式（＊18）より、式（＊24）が成り立つ。
（＊18）再掲
\begin{eqnarray*}
\alpha' &\equiv& 0\pmod{2}\\
a' &\equiv& 1\pmod{8}\\
b' &\equiv& 2\pmod{4}
\end{eqnarray*}
（＊24）
\begin{eqnarray*}
\alpha'' &\equiv& 0\pmod{4}\\
a'' &\equiv& 1\pmod{16}\\
b'' &\equiv& 4\pmod{8}
\end{eqnarray*}
式（＊22）（＊23）より、次式が成り立つ。
（＊25）
\[
\alpha''=b''\sum_{k=0}^{n''-1}\left(a''\right)^{k}
\]
\(n''\) の最下位ビットを \(n''_{0}\,\left(\in\left\{0,1\right\}\right)\) として \(n''\overset{\mathrm{def}}{=}2n'''+n''_{0}\;\left(n'''\ge 0\right)\) とおくと、次式が成り立つ。
（＊26）
\begin{eqnarray*}
\alpha''
&=& b''\left(n''_{0}+\left(a''\right)^{n''_{0}}\sum_{k=0}^{2n'''-1}\left(a''\right)^{k}\right)\\
&=& b''\left(n''_{0}+\left(a''\right)^{n''_{0}}\sum_{k=0}^{n'''-1}\left(\left(a''\right)^{2k}+\left(a''\right)^{2k+1}\right)\right)\\
&=& b''\left(n''_{0}+\left(a''\right)^{n''_{0}}\left(1+a''\right)\sum_{k=0}^{n'''-1}\left(a''\right)^{2k}\right)\\
\Longleftrightarrow\;
\alpha''-b''n''_{0}&=& b''\left(a''\right)^{n''_{0}}\left(1+a''\right)\sum_{k=0}^{n'''-1}\left(a''\right)^{2k}
\end{eqnarray*}
式（＊24）より、式（＊26）の右辺は \(\left(\text{奇数の\(8\)倍}\right)\cdot\sum\left(\text{略}\right)\) の形であり、また左辺は \(\alpha''\) が \(4\) の倍数， \(b''\) が奇数の \(4\) 倍である。
よって、 \(\alpha''\) の第２ビット \(\alpha''_{2}\) は \(n''_{0}\;\left(=n'_{1}=n_{2}\right)\) に等しい。
（＊27）
\begin{eqnarray*}
\alpha''_{2}
&=& n''_{0}\\
&=& n_{2}
\end{eqnarray*}
よって式（＊26）（＊27）より、
（＊28）
\[
\alpha''-b''\alpha''_{2}=b''\left(a''\right)^{\alpha''_{2}}\left(1+a''\right)\sum_{k=0}^{n'''-1}\left(a''\right)^{2k}
\]
ここで、 \(\alpha''',\,a''',\,b'''\) を次のように定義する。
（＊29）
\begin{eqnarray*}
\alpha''' &\overset{\mathrm{def}}{=}& \alpha''-b''\alpha''_{2}\\
a''' &\overset{\mathrm{def}}{=}& \left(a''\right)^{2}\\
b''' &\overset{\mathrm{def}}{=}& b''\left(a''\right)^{\alpha''_{2}}\left(1+a''\right)
\end{eqnarray*}
このとき、 \(\alpha''_{2}\) が \(\alpha''\) の第２ビットであることと式（＊24）より、式（＊30）が成り立つ。
（＊24）再掲
\begin{eqnarray*}
\alpha'' &\equiv& 0\pmod{4}\\
a'' &\equiv& 1\pmod{16}\\
b'' &\equiv& 4\pmod{8}
\end{eqnarray*}
（＊30）
\begin{eqnarray*}
\alpha''' &\equiv& 0\pmod{8}\\
a''' &\equiv& 1\pmod{32}\\
b''' &\equiv& 8\pmod{16}
\end{eqnarray*}
式（＊28）（＊29）より、次式が成り立つ。
（＊31）
\[
\alpha'''=b'''\sum_{k=0}^{n'''-1}\left(a'''\right)^{k}
\]
以下同様に繰り返し計算していくと、下位ビットから順に１桁ずつ \(0\) になっていくので、いずれ左辺は \(0\) になり、このアルゴリズムを終了することができる。
そしてそれまでに得られた \(n_{0},\,n_{1},\,n_{2},\,\dots\) から \(n\) が得られる。
この方法で \(n\) を求めるプログラムは後述する。
【解２】
【解１】の途中、式（＊16）において、両辺に \(\left(a^{-1}\right)^{\alpha_{0}}\) をかけると、式（＊32）のようになる。
（＊16）再掲
\[
\alpha-b\alpha_{0}=ba^{\alpha_{0}}\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}
\]
（＊32）
\begin{eqnarray*}
\left(a^{-1}\right)^{\alpha_{0}}\cdot\left(\alpha-b\alpha_{0}\right)
&=&\left(a^{-1}\right)^{\alpha_{0}}\cdot ba^{\alpha_{0}}\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}\\
&=&\left(a^{-1}a\right)^{\alpha_{0}}b\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}\\
&=& b\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}\\
\Longleftrightarrow\;
\left(a^{-1}\right)^{\alpha_{0}}\alpha-\left(a^{-1}\right)^{\alpha_{0}}b\alpha_{0}
&=& b\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}\\
\Longleftrightarrow\;
\left(f^{-1}\right)^{\alpha_{0}}\left(\alpha\right)
&=& b\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}\\
\end{eqnarray*}
記述の短縮のために \(f^{-m}\left(x\right)\overset{\mathrm{def}}{=}\left(f^{-1}\right)^{m}\left(x\right)\) とすると、
（＊33）
\[
f^{-\alpha_{0}}\left(\alpha\right)=b\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}
\]
この変形によって、実質的な場合分けに当たる \(\alpha_{0}\) を含む項を１箇所にまとめることができた。
以下【解１】のときと同じように、 \(\alpha',\,a',\,b'\) および関数 \(f'\) を次のように定義する。
（＊34）
\begin{eqnarray*}
\alpha' &\overset{\mathrm{def}}{=}& f^{-\alpha_{0}}\left(\alpha\right)\\
a' &\overset{\mathrm{def}}{=}& a^{2}\\
b' &\overset{\mathrm{def}}{=}& b\left(1+a\right)\\
f'\left(x\right) &\overset{\mathrm{def}}{=}& a'x+b'\\
&=& a^{2}x+b\left(1+a\right)
\end{eqnarray*}
このとき式（＊５）より、 \(f'\) は \(f^{2}\) に等しい。
（＊５）再掲
\[
f^{2}\left(x\right)=a^{2}x+b\left(1+a\right)
\]
ところで式（＊34）を式（＊33）に代入すると、次式が成り立つ。
（＊35）
\[
\alpha'=b'\sum_{k=0}^{n'-1}\left(a'\right)^{k}
\]
以下、【解１】と同様の繰り返し計算を行うことで、 \(n\) が求められる。
このとき、 \(f^{-1},\,f^{-2},\,f^{-4},\,f^{-8},\,\dots\) が必要になるが、その係数はいずれも \(a,\,b\) のみに依存する定数なので、予め求めておくことができる。
つまり、【解１】より速く \(n\) を求められる。
この方法で \(n\) を求めるプログラムは後述する。
【解３】
（＊10）再掲
\[
\alpha=b\sum_{k=0}^{n-1}a^{k}
\]
式（＊10）において、 \(n\) の最下位ビットを \(n_{0}\,\left(\in\left\{0,1\right\}\right)\) として \(n\overset{\mathrm{def}}{=}2n'-n_{0}\;\left(n'\ge 1\right)\) とおくと、次式が成り立つ。
（＊36）
\begin{eqnarray*}
\alpha
&=& b\sum_{k=0}^{n-1}a^{k}\\
% &=& b\sum_{k=0}^{2n'-n_{0}-1}a^{k}\\
% &=& \left(a^{-1}\right)^{n_{0}}b\sum_{k=n_{0}}^{2n'-1}a^{k}\\
&=& \left(a^{-1}\right)^{n_{0}}b\left(-n_{0}+\sum_{k=0}^{2n'-1}a^{k}\right)\\
\Longleftrightarrow\;
a^{n_{0}}\alpha+bn_{0}
&=& b\sum_{k=0}^{2n'-1}a^{k}\\
&=& b\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}\\
\Longleftrightarrow\;
f^{n_{0}}\left(\alpha\right)
&=& b\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}
\end{eqnarray*}
【解１】と同様に偶奇を調べると、 \(\alpha\) の最下位ビット \(\alpha_{0}\) が \(n_{0}\) に等しいことが分かる。
よって、
（＊37）
\[
f^{\alpha_{0}}\left(\alpha\right)=b\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}
\]
式（＊37）は【解２】の式（＊33）とほぼ同じ式だが、唯一の違いは \(f^{-1}\) の部分が \(f\) になった点である。
（＊33）再掲
\[
f^{-\alpha_{0}}\left(\alpha\right)=b\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}
\]
式の形が【解２】とほぼ同じなので、以降の解き方も【解２】とほぼ同じである。
以下【解２】のときと同じように、 \(\alpha',\,a',\,b'\) および関数 \(f'\) を次のように定義する。
（＊38）
\begin{eqnarray*}
\alpha' &\overset{\mathrm{def}}{=}& f^{\alpha_{0}}\left(\alpha\right)\\
a' &\overset{\mathrm{def}}{=}& a^{2}\\
b' &\overset{\mathrm{def}}{=}& b\left(1+a\right)\\
f'\left(x\right) &\overset{\mathrm{def}}{=}& a'x+b'\\
&=& a^{2}x+b\left(1+a\right)
\end{eqnarray*}
式（＊37）（＊38）より、次式が成り立つ。
（＊39）
\[
\alpha'=b'\sum_{k=0}^{n'-1}\left(a'\right)^{k}
\]
以下、【解２】と同様の繰り返し計算を行うことで、 \(n_{0},\,n_{1},\,n_{2},\,\dots\) が順に求められる。
これらを順に並べてできる数を \(m\;\left(\overset{\mathrm{def}}{=}\sum\limits_{i=0}^{N}2^{i}n_{i}\right)\) とする。
このとき \(m\) は、題意の \(\alpha=f^{m}\left(0\right)\) を満たす \(m\) ではなく、 \(f^{m}\left(\alpha\right)=0\) を満たす \(m\) である。
言い換えると、 \(\alpha=f^{M-m}\left(0\right)\) を満たす \(m\) である。
求めるべき値は \(\alpha=f^{n}\left(0\right)\) を満たす \(n\) なので、 \(n=M-m\) によって \(n\) が得られる。
なお 計算過程で \(f,\,f^{2},\,f^{4},\,f^{8},\,\dots\) が必要になるが、その係数はいずれも \(a,\,b\) のみに依存する定数なので、予め求めておくことができる。
\(f^{-1},\,f^{-2},\,f^{-4},\,f^{-8},\,\dots\) が必要だった【解２】に比べて自然な解法と言えるだろう。タブンネ。
この方法で \(n\) を求めるプログラムは後述する。
\(\alpha\) に \(f\) を何回施すと \(\beta\) になるか？
（★２）【問２】
ある整数 \(\alpha,\,\beta\) が与えられたとき、次式を満たす非負整数 \(n\left(< M\right)\) を求めよ。
\[
\beta=f^{n}\left(\alpha\right)
\]
【解１】
次の関係を満たすように \(l,\,m\) を定める。
（＊40）
\begin{eqnarray*}
\alpha &\overset{\mathrm{def}}{=}& f^{l}\left(0\right)\\
\beta &\overset{\mathrm{def}}{=}& f^{m}\left(0\right)
\end{eqnarray*}
式（＊40）の第１式を式（★２）に代入して、
（＊41）
\begin{eqnarray*}
\beta
&=& f^{n}\left(\alpha\right)\\
&=& f^{n}\left(f^{l}\left(0\right)\right)\\
&=& f^{n+l}\left(0\right)
\end{eqnarray*}
よって式（＊40）（＊41）より \(m=n+l\) なので、
（＊42）
\[
n=m-l
\]
\(m,\,l\) は【問１】の方法で求められるので、求めたい \(n\) も求められる。
【解２】
【解１】の方法は \(l,\,m\) を別個に求めているので \(O\left(N\right)\) の計算が２回必要である。
これはなんとなく無駄な気がするので、なんとか \(O\left(N\right)\) の計算１回で \(n\) を求められないだろうか？
模索してみる。
（※ちなみに結論を先に言うと無理でした。）
式（＊41）をさらに変形すると、
（＊43）
\begin{eqnarray*}
\beta
&=& f^{n+l}\left(0\right)\\
&=& f^{l}\left(f^{n}\left(0\right)\right)\\
&=& a^{l}f^{n}\left(0\right)+f^{l}\left(0\right)\qquad\left(\text{∵式（＊７）}\right)\\
&=& a^{l}f^{n}\left(0\right)+\alpha\\
\Longleftrightarrow\;
\beta-\alpha
&=& a^{l}f^{n}\left(0\right)
\end{eqnarray*}
（＊７）再掲
\[
f^{n}\left(x\right)=a^{n}x+f^{n}\left(0\right)
\]
式（＊43）の両辺に \(\left(a^{-1}\right)^{l}\) をかけて、
（＊44）
\[
\left(a^{-1}\right)^{l}\left(\beta-\alpha\right)=f^{n}\left(0\right)
\]
式（＊44）において、 \(\alpha,\,\beta\) は既知であり、 \(l\) も式（＊40）より \(\alpha=f^{l}\left(0\right)\) なので【問１】の方法で求められる。
よって式（＊44）の左辺を先に求めれば、再び【問１】の方法を使って \(n\) が求められる。
……しかし結局のところ【問１】の方法を２回使っていて、その上に余計な計算が増えている。
残念ながら【解２】は【解１】より劣る方法であると言わざるをえない。
【解３】
【解２】で、 \(l\) を求めるのに【問１】の方法を使うのは悪手っぽいことがわかった。
そこでいっそ \(l\) を消去してしまう方法を考えてみる。
そのために、式（＊43）に式（＊８）を適用する。
（＊８）再掲
\[
a^{n}=\left(a-1\right)b^{-1}f^{n}\left(0\right)+1
\]
（＊45）
\begin{eqnarray*}
\beta-\alpha
&=& a^{l}f^{n}\left(0\right)\\
&=& \left(\left(a-1\right)b^{-1}f^{l}\left(0\right)+1\right)f^{n}\left(0\right)\\
&=& \left(\left(a-1\right)b^{-1}\alpha+1\right)f^{n}\left(0\right)
\end{eqnarray*}
\(a-1\) が偶数ゆえに \(\left(a-1\right)b^{-1}\alpha+1\) は奇数なので、逆元 \(\left(\left(a-1\right)b^{-1}\alpha+1\right)^{-1}\) が存在して、
（＊46）
\[
\left(\left(a-1\right)b^{-1}\alpha+1\right)^{-1}\left(\beta-\alpha\right)=f^{n}\left(0\right)
\]
これにより、【問１】の方法を使う回数を１回に減らすことができた。
しかし逆元を求める計算も \(O\left(N\right)\) であり、結局 \(O\left(N\right)\) の計算が２回必要なことは変わらない。
【解２】と同様、余計な計算が必要になる【解３】も【解１】より劣る方法だと言えるだろう。
結局のところ、単純に【問１】の方法で \(l,\,m\) を個別に求めて差分を取るのが１番いいようだ。
しんぷるいずべすとってコトネ。
【問１】実装
【問１】の３つの解法を実装してみた。
using System; using System.Collections.Generic; using System.Text; // C# namespace Tmp { public class LCGTest { #region static constructor static LCGTest() { // a の逆元 c を求める c = 0ul; ulong x = 1ul; ulong multiplier = a; for( ulong mask=1ul;mask!=0ul;mask<<=1,multiplier<<=1 ) { if( ( x & mask ) != 0ul ) { c |= mask; x -= multiplier; } } // f^{-1}(x) = c * x + d を満たす d を求める d = 0ul - c * b; // f^{n}(x) = aArr[n] * x + bArr[n] を満たす aArr, bArr を求める // f^{M-n}(x) = cArr[n] * x + dArr[n] を満たす cArr, dArr を求める aArr[0] = a; bArr[0] = b; cArr[0] = c; dArr[0] = d; for( int i=1;i<N;i++ ) { aArr[i] = aArr[i-1] * aArr[i-1]; bArr[i] = bArr[i-1] * ( aArr[i-1] + 1 ); cArr[i] = cArr[i-1] * cArr[i-1]; dArr[i] = dArr[i-1] * ( cArr[i-1] + 1 ); } } #endregion #region static field private const int N = 64; private const ulong a = 0x5d588b656c078965ul; private const ulong b = 0x269ec3ul; private static readonly ulong c; // = a^{-1} private static readonly ulong d; // f^{-1}(x) = c * x + d private static readonly ulong[] aArr = new ulong[N]; private static readonly ulong[] bArr = new ulong[N]; private static readonly ulong[] cArr = new ulong[N]; private static readonly ulong[] dArr = new ulong[N]; #endregion #region static method public static ulong Next(ulong x) { return x * a + b; } public static ulong Prev(ulong x) { return x * c + d; } // 【解１】の方法 public static ulong CalcN1(ulong x) { ulong n = 0ul; ulong multiplier = a; ulong adder = b; for( ulong mask=1ul;mask!=0ul;mask<<=1 ) { if( ( x & mask ) != 0ul ) { n |= mask; x -= adder; adder *= multiplier; } adder *= multiplier + 1; multiplier *= multiplier; } return n; } // 【解２】の方法 public static ulong CalcN2(ulong x) { ulong n = 0ul; for( int i=0;i<N;i++ ) { ulong mask = 1ul << i; if( ( x & mask ) != 0ul ) { n |= mask; x = x * cArr[i] + dArr[i]; // x = f^{-2^{i}}(x) } } return n; } // 【解３】の方法 public static ulong CalcN3(ulong x) { ulong m = 0ul; for( int i=0;i<N;i++ ) { ulong mask = 1ul << i; if( ( x & mask ) != 0ul ) { m |= mask; x = x * aArr[i] + bArr[i]; // x = f^{2^{i}}(x) } } ulong n = 0ul - m; // n = M - m = 2^{N} - m ≡ 0 - m (mod M) return n; } #endregion } }
