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

線形合同法(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\) を次のように定義する。

(*1)
\[
f\left(x\right)\overset{\mathrm{def}}{=}ax+b
\]

ただし \(a,\,b\) は次式の関係を満たすものとする。

(*2)
\begin{eqnarray*}
a&\equiv& 1\pmod{4}\\
b&\equiv& 1\pmod{2}
\end{eqnarray*}


このとき、 定数 \(a,\,b\) および関数 \(f\) に関して、以下の関係や性質が成り立つ。

\(a,\,b\) の逆元 \(a^{-1},\,b^{-1}\)

式(*2)より \(a,\,b\) はどちらも \(M\left(=2^{N}\right)\) と互いに素なので、逆元が存在する。
よって以降、 \(a,\,b\) の逆元をそれぞれ \(a^{-1},\,b^{-1}\) とする。

(*3)
\[
\begin{array}{lclcl}
aa^{-1} &=& a^{-1}a &=& 1\\
bb^{-1} &=& b^{-1}b &=& 1
\end{array}
\]

\(f^{n}\left(x\right)\)

(*4)
\(n\ge 0\) なる \(n\) について、次式が成り立つ。
\[
f^{n}\left(x\right)=a^{n}x+b\sum_{k=0}^{n-1}a^{k}
\]

【証明】
帰納法で。


特に、 \(n=2\) のときは、

(*5)
\[
f^{2}\left(x\right)=a^{2}x+b\left(1+a\right)
\]

\(f^{n}\left(0\right)\)

(*6)
\(n\ge 0\) なる \(n\) について、次式が成り立つ。
\[
f^{n}\left(0\right)=b\sum_{k=0}^{n-1}a^{k}
\]

【証明】
式(*4)に \(x=0\) を代入。


式(*4)(*6)より、 \(n\ge 0\) のとき次式が成り立つ。

(*7)
\[
f^{n}\left(x\right)=a^{n}x+f^{n}\left(0\right)
\]

\(a^{n}\)

(*8)
\(n\ge 0\) なる \(n\) について、次式が成り立つ。
\[
a^{n}=\left(a-1\right)b^{-1}f^{n}\left(0\right)+1
\]

【証明】
式(*6)の両辺に \(\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*}

なお式(*8)において、 \(a-1\) は偶数ゆえに逆元が存在しないので、 \(f^{n}\left(0\right)=\cdots\) の形に変形しようなどと考えてはならない。
法 \(M\) と互いに素でない \(a-1\) をかけたときに同値性が崩れているのだ。

\(f^{-1}\left(x\right)\)

(*9)
\[
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\) となり式(*9)を得る。

\(0\) に \(f\) を何回施すと \(\alpha\) になるか?

ここまでは準備。
ここからが本題。

(★1)【問1】
ある整数 \(\alpha\) が与えられたとき、次式を満たす非負整数 \(n\left(< M\right)\) を求めよ。
\[
\alpha=f^{n}\left(0\right)
\]

【解1】

式(★1)(*6)より、式(*10)が成り立つ。

(*6)再掲
\[
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*}

(*2)再掲
\begin{eqnarray*}
a&\equiv& 1\pmod{4}\\
b&\equiv& 1\pmod{2}
\end{eqnarray*}

式(*2)より \(a\) は奇数なので \(1+a\) は偶数であり、ゆえに式(*14)の右辺も偶数である。
よって、 \(\alpha\) の偶奇は \(bn_{0}\) の偶奇に等しい。
式(*2)より \(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\) の最下位ビットであることと式(*2)より、次の関係が成り立つ。

(*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'\) の第1ビット \(\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'\) の第1ビットであることと式(*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''\) の第2ビット \(\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''\) の第2ビットであることと式(*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}
\]

以下同様に繰り返し計算していくと、下位ビットから順に1桁ずつ \(0\) になっていくので、いずれ左辺は \(0\) になり、このアルゴリズムを終了することができる。
そしてそれまでに得られた \(n_{0},\,n_{1},\,n_{2},\,\dots\) から \(n\) が得られる。

この方法で \(n\) を求めるプログラムは後述する。

【解2】

【解1】の途中、式(*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}\) を含む項を1箇所にまとめることができた。
以下【解1】のときと同じように、 \(\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*}

このとき式(*5)より、 \(f'\) は \(f^{2}\) に等しい。

(*5)再掲
\[
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}
\]

以下、【解1】と同様の繰り返し計算を行うことで、 \(n\) が求められる。
このとき、 \(f^{-1},\,f^{-2},\,f^{-4},\,f^{-8},\,\dots\) が必要になるが、その係数はいずれも \(a,\,b\) のみに依存する定数なので、予め求めておくことができる。
つまり、【解1】より速く \(n\) を求められる。

この方法で \(n\) を求めるプログラムは後述する。

【解3】

(*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*}

【解1】と同様に偶奇を調べると、 \(\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)は【解2】の式(*33)とほぼ同じ式だが、唯一の違いは \(f^{-1}\) の部分が \(f\) になった点である。

(*33)再掲
\[
f^{-\alpha_{0}}\left(\alpha\right)=b\left(1+a\right)\sum_{k=0}^{n'-1}a^{2k}
\]

式の形が【解2】とほぼ同じなので、以降の解き方も【解2】とほぼ同じである。
以下【解2】のときと同じように、 \(\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}
\]

以下、【解2】と同様の繰り返し計算を行うことで、 \(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\) が必要だった【解2】に比べて自然な解法と言えるだろう。タブンネ

この方法で \(n\) を求めるプログラムは後述する。

\(\alpha\) に \(f\) を何回施すと \(\beta\) になるか?

(★2)【問2】
ある整数 \(\alpha,\,\beta\) が与えられたとき、次式を満たす非負整数 \(n\left(< M\right)\) を求めよ。
\[
\beta=f^{n}\left(\alpha\right)
\]

【解1】

次の関係を満たすように \(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)の第1式を式(★2)に代入して、

(*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\) は【問1】の方法で求められるので、求めたい \(n\) も求められる。

【解2】

【解1】の方法は \(l,\,m\) を別個に求めているので \(O\left(N\right)\) の計算が2回必要である。
これはなんとなく無駄な気がするので、なんとか \(O\left(N\right)\) の計算1回で \(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{∵式(*7)}\right)\\
&=& a^{l}f^{n}\left(0\right)+\alpha\\
\Longleftrightarrow\;
\beta-\alpha
&=& a^{l}f^{n}\left(0\right)
\end{eqnarray*}

(*7)再掲
\[
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)\) なので【問1】の方法で求められる。
よって式(*44)の左辺を先に求めれば、再び【問1】の方法を使って \(n\) が求められる。
……しかし結局のところ【問1】の方法を2回使っていて、その上に余計な計算が増えている。
残念ながら【解2】は【解1】より劣る方法であると言わざるをえない。

【解3】

【解2】で、 \(l\) を求めるのに【問1】の方法を使うのは悪手っぽいことがわかった。
そこでいっそ \(l\) を消去してしまう方法を考えてみる。
そのために、式(*43)に式(*8)を適用する。

(*8)再掲
\[
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)
\]

これにより、【問1】の方法を使う回数を1回に減らすことができた。
しかし逆元を求める計算も \(O\left(N\right)\) であり、結局 \(O\left(N\right)\) の計算が2回必要なことは変わらない。
【解2】と同様、余計な計算が必要になる【解3】も【解1】より劣る方法だと言えるだろう。


結局のところ、単純に【問1】の方法で \(l,\,m\) を個別に求めて差分を取るのが1番いいようだ。
しんぷるいずべすとってコトネ。

【問1】実装

【問1】の3つの解法を実装してみた。

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;
		}
		// 【解1】の方法
		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;
		}
		// 【解2】の方法
		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;
		}
		// 【解3】の方法
		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

	}
}