author: fujidig, created at: 2017/12/9
数学とコンピュータⅡ Advent Calendar 2017 9日目の記事です。
メルセンヌツイスタは疑似乱数の一つです。 特定の形をした線形写像$f: \F_2^p \to \F_2^p$が一つ固定され、与えられた初期状態$\bm{x}_0 \in \F_2^p$に対し次の漸化式によって内部状態が更新されます。 \[ \bm{x}_{n+1} = f(\bm{x}_n) \]
内部状態からある線形写像$g: \F_2^p \to \F_2^w$により \[ \bm{y}_n = g(\bm{x}_n)\] により疑似乱数列が生成されます。
mt19937というメルセンヌツイスタの一例では$p = 19937, w = 32$が使われています。
線形写像$f, g$はビットシフトなどの方法でコンピュータ上で高速に実装が可能なものが採用されており、その計算は$O(1)$時間でできます。
$p$として$2^p-1$が素数あるもの、すなわち$p$がメルセンヌ指数であるものを採用していることがポイントです。 これにより最長の周期 (すなわち周期が$2^p-1$)を達することが比較的簡単にチェックすることができます。 実際、周期が最長となるパラメータはその方法によって見つけられています。 今回はその方法について見ていきます。
周期の最長性を確認する最もナイーブな方法は$\bm{0}$以外のある$\bm{x}_0$をとり, $\bm{x}_1, \bm{x}_2, \dots$を順番に計算していき$\bm{x}_{2^p-1}$で初めて$\bm{x}_0$に戻ってくることを確認する方法です。
これは$O(2^p)$時間かかりとても現実的な時間ではできません。
次の定理があります。
これを使えば$f$の行列の乗法位数が$2^p-1$になることをチェックすればいいわけです。 今、$2^p-1$は素数なので、$2^p-1$乗が$E$になることさえチェックできれば十分です。 $p$次行列の掛け算が$O(p^3)$時間でできることを考えれば、このチェックは$O(p^4)$時間でできます。
方法1より改善しましたが、まだ時間がかかります。
ここに原始多項式とは次で定義されるものです。
よって$f$の特性多項式$\varphi(t)$を計算し、$K[t]/(\varphi(t))$において$t^{2^p - 1} = 1$であることをチェックすればよいです。 次数が$p$以下の多項式同士の掛け算は$O(p^2)$時間でできることからその計算量は$O(p^3)$時間です。
また、計算量は変わりませんが次も使えます。
実は$O(p^2)$のアルゴリズムが存在します。 アルゴリズムを述べる前に準備をします。
$S^{\infty}$を$0, 1$からなる無限列全体のなす$\F_2$線形空間とします。すなわち \[S^{\infty} = \{ \chi = (\dots, x_5, x_4, x_3, x_2, x_1, x_0) \mid x_i \in \F_2 \}. \]
$D$ (delay operator)と$H$ (decimation operator)を次で定める$S^\infty$から$S^\infty$への線形作用素とします: \begin{align*} D(\dots, x_4, x_3, x_2, x_1, x_0) &= (\dots, x_5, x_4, x_3, x_2, x_1) \\ H(\dots, x_4, x_3, x_2, x_1, x_0) &= (\dots, x_8, x_6, x_4, x_2, x_0) \end{align*}
定理4
$\varphi(t)$が$\F_2$上の多項式で次数$p$がメルセンヌ指数だとする。
$\chi \in S^\infty$であって、$\varphi(D)\chi = 0, H\chi \ne \chi$であるものをとる。
このとき、$\varphi(t)$が原始的であるためには$H^p \chi = \chi$であることが必要十分。
この定理を用いると、$O(p^2)$時間のアルゴリズムの存在が示せます。ただし、ここまでには述べていなかった仮定が必要です。
定理5
疑似乱数生成器の状態空間$V$を$\F_2$上の$p$次元ベクトル空間とする。ただし$p$はメルセンヌ指数。
$f: V \to V$を線形な状態遷移写像とする。
$b: V \to \F_2$を線形写像とする ($b$はたとえば状態から$1$ビットlook upするもの)。
$f$と$b$は$O(1)$時間で計算できるものとする。
写像$\Phi: V \to \F_2^p$: \[\Phi: S \mapsto (bf^{p-1}(S), bf^{p-2}(S), \dots, bf(S), b(S))\] を全単射だと仮定する。 そしてこの逆写像は$O(p)$時間で計算可能とする。
このとき、$f$の特性多項式の原始性は$O(p^2)$時間でテストできる。
証明。 $\chi$を無限列$(\dots, bf^2(S), bf(S), b(S))$とする。 $H(\chi) \ne \chi$となる$S$を選ぶ。
$H^m(\chi)$の最初の$p$ビットから$H^{m+1}(\chi)$の最初の$p$ビットを$O(p)$時間で求められることを示す。
$H^m(\chi)$の最初の$p$ビットから、$H^m(\chi)$を生成する状態$S_m$を得ることができる ($\Phi$の逆関数を使って$O(p)$時間で)。 この状態から$H^m(\chi)$の最初の$2p$ビットを$O(2p)$時間で計算できる ($f$と$b$が$O(1)$なので)。 ここでこの$2p$ビットの列から奇数項を間引き、$p$ビットの列を作る。 すると$H^{m+1}(\chi)$の最初の$p$ビットを$O(p)$時間で得る。
以上の方法を繰り返せば$O(p^2)$時間で$H^p(\chi)$の最初の$p$ビットを得られる。 あとは、定理4より、$H^p(\chi) = \chi$が成立することをチェックすればよい。 □
この定理のポイントは$\Phi$の逆関数の計算を$O(p)$時間でできるような$b: V \to \F_2$があると仮定した点です。 メルセンヌツイスタの状態遷移写像$f$の詳細な形を見てその仮定が満たされることはチェック可能ですが、この記事では省くことにします。
[1] 松本眞 (2004), 「有限体の擬似乱数への応用」 (http://www.math.sci.hiroshima-u.ac.jp/~m-mat/TEACH/0407-2.pdf)
[2] M. Matsumoto and T. Nishimura (1998), "Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator", ACM Trans. on Modeling and Computer Simulation Vol. 8, No. 1, January pp.3-30