メルセンヌ・ツイスタ（Mersenne Twister, MT）をご存知でしょうか？
ポケモンやその他のゲームで乱数調整などをする人ならば名前を聞いたり、あるいは使いこなしてる方もいるかもしれない、擬似乱数生成アルゴリズムの一つですね。
私はPythonをよく使うのですが、numpy.randomモジュールのrand関数の中に入っているのがこのMTだったりもします。
擬似乱数生成法のスタンダードの一つといってもいいのかなと思うくらいには浸透していると思います。
研究室での活動の中でMTを使いたい要請があり、少し勉強していたのですが、MTの周期などについて解説記事や文献が少し見つけにくく、特に周期に関する証明を見つけるのにかなり苦労しました。
せっかく勉強したんだから、ここは一つ記事にまとめようと思った次第です。
この記事では、MTの前身である線形帰還シフトレジスタ（Linear Feedback Shift Register, LFSR）の周期とその証明を示し、MTの周期・高次元分布について概観します。
ちょっとだけ数学的に真面目にやるので、線形代数についての基礎的な知識はあったほうが読みやすいかも知れません。
MTとは
本編に入る前にMTってそもそもなんぞやということについてもう少しだけ述べておきます。MTは
- $2^{19937} - 1$という超長周期
- 32ビット精度で623次元に均等分布（後述）
- 高速な実装が可能
という点で非常に優れた乱数です。
一応デメリットも挙げておきますと、ワーキングメモリが比較的大きく、そこまで高いクオリティの乱数が必要とされない場合にはオーバースペックとなりがちなことです。
大きいと言っても、一般的なCPUでは気になるほどのサイズではないのですが、例えば組み込みで小規模に実装しようと考えた際などには、本当にMTが必要か検討する必要があるかもしれません。
メルセンヌ・ツイスターは必要か？も一理あり、大規模シミュレーションなどでもない限りは別の乱数でもいいのかも知れませんが、今のPCのスペックを考えると誤差という気もします。
参考サイト
あなたの使っている乱数、大丈夫？
こちらも製作者の松本さんによる、広島大学での講演スライドを公開してくださっているものです。
一般向けに非常に分かりやすく書かれています。
下の論文の裏話などもあって読んでて面白いです。
Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator
元論文です。歴史的経緯とかが簡潔にまとめられています。MTのメリットなどは前身のTwisted GFSR Generatorsの方が分かりやすかったりします。
証明は私は追えませんでした。イメージが豊富なので専門家には優しいのだと思います。
有限体の擬似乱数への応用
詳細は分からないですが、おそらくどこかの大学で松本先生が講義された際のレジュメっぽい雰囲気があります。
数学的にきっちりまとめてくださっているので、今回の記事作成にあたって非常に参考になりました。
線形合同法（LCG）
さて、本題のMTに入る前にMT以前に用いられていた乱数である線形合同法（Linear congruential generator, LCG）について簡単に見ておきます。
LCGは漸化式
$$
x_{i+1} = a x_i + b \mathrm{~mod~} m \quad (a \ge 1, b \ge 0 ,m \ge 2 \mathrm{:~定数})
$$
により乱数の系列を作るものです。
$a=3, b=1, m=9, x_0=1$として実際に数列を生成してみると、
\begin{align}
x_0 =& 1\\
x_1 =& (4 \times 1 + 7) \mathrm{~mod~} 9 = 2\\
x_2 =& (4 \times 2 + 7) \mathrm{~mod~} 9 = 6\\
x_3 =& (4 \times 6 + 7) \mathrm{~mod~} 9 = 4\\
x_4 =& (4 \times 4 + 7) \mathrm{~mod~} 9 = 5\\
x_5 =& (4 \times 5 + 7) \mathrm{~mod~} 9 = 0\\
x_6 =& (4 \times 0 + 7) \mathrm{~mod~} 9 = 7\\
x_7 =& (4 \times 7 + 7) \mathrm{~mod~} 9 = 8\\
x_8 =& (4 \times 8 + 7) \mathrm{~mod~} 9 = 3\\
x_9 =& (4 \times 3 + 7) \mathrm{~mod~} 9 = 1\\
x_{10} =& (4 \times 1 + 7) \mathrm{~mod~} 9 = 2\\
&\hspace{2cm} \vdots
\end{align}
のようになり、確かに0から9までの数字がランダムに見えるような順番で出現することが分かります。
このような、漸化式によって得られる数列を擬似乱数列と言います。
擬似乱数の周期を、次のように定義します。
定義（擬似乱数の周期）
擬似乱数列$\{ x_i \}_{i=0,1,\cdots}$について、ある$d \ge 0$が存在して、任意の$j \ge d$について$x_j = x_{j+r}$が成り立つような$r ~ (0 < r \in \mathbb{Z})$のうち、最小のものを${ x_i }$の周期と言う。
このとき、LCGを含むような一階漸化式が生成する擬似乱数（一つ前の乱数のみ参照するもの）の周期について次のことが分かります。
命題（一階漸化式が生成する擬似乱数列の周期）
擬似乱数列$\{x_i \}$を
$$
x_{i+1} = f(x_i) \mathrm{~mod~} m
$$
により生成する。($f: \mathbb{Z}\rightarrow\mathbb{Z}$)このとき、${x_i}$の周期は$m$を超えない。
（証明）
$x_i$の値は0から$m-1$までの$m$個の整数のいずれかであるので、$x_0$から$x_m$までの$m+1$個の数字のなかに、$x_{i} = x_j$となるような$i,j~(i<j)$が必ず存在します。
$k > i$について、$x_k$は$x_i$に$f$を作用させて$\mathrm{mod~}m$を取る操作を$(k-i)$回くりかえしたものです。
同様に、$x_{k+(j-i)} = x_{j+(k-i)}$は$x_j$に$f$を作用させて$\mathrm{mod~}m$を取る操作を$(k-i)$回くりかえしたものです。
$x_i = x_j$よりこの2数は同じ値で、$x_{k} = x_{k+(j-i)}$が成り立ちます。
上の議論は$i$より大きい任意の$k$について可能ですので、定義より$(j-i)$は$\{ x_i \}$の周期よりも大きいです。
$0\le i < j \le m$より、周期が$m$を超えないことが従います。
（証明終）
面倒くさく書いていますが、要は「一階の漸化式を用いて擬似乱数列を作る場合、同じ数字が出現したらそのあとは前に出現したパターンが繰り返し出現する」ということを言っているだけです。
この節の最初に示した例で言えば、$x_0$、$x_9$がともに1ですので、$x_9$以降同じパターンの繰り返し（1, 2, 6, 4, $\cdots$）になります。
さて、LCGは一階漸化式を用いているため、最長の周期は$m$となります。
乱数の周期は長いことが望ましいのですが、LCGの場合、周期（$< m$）と出現する数字の範囲（$0 \le x_i \le m-1$）が相関しているのが非常にイケてないところです。
つまり、0から255の範囲の乱数が欲しいと言われると、周期は最大でも256となってしまう、ということになります。
もちろん、「より大きな乱数を用意してそれをさらにmodをとって小さな値の乱数を得る」といった、ちょっとした工夫によって周期を伸ばすことはできるものの、LCGの下位ビットの乱数性は非常に悪いこと、また$m$の値を大きくすれば計算量・消費メモリが増大することなどからあまりうまい方法ではありません。
さらにとどめを刺すのが上でも少し言葉が出てきた高次元での均等分布です。
高次元での均等分布とは、連続する$k$個の数を組にしたベクトルを考えたときに、それが$k$次元空間を一様に満たすということを表します。
定義は次の通りです。
定義（高次元均等分布）
（有限）集合$X$上に値を取る擬似乱数列$\{ x_i \}_{i=0,1,\cdots}$を考える。
いま、$ \{ x_i \} $の連続する$k$個の組の列$ \{ [ x_i, x_{i+1}, \cdots ,x_{i+k-1} ] \in X^k \}_{i=0,1,\cdots}$を一周期にわたってとったとき、$X^k$上の全ての点が同じ回数ずつ出現するならば、$\{ x_i \}$は$k$次元に均等分布しているという。
（ただし、全て0という組のみは特別に他より1回少なくても良いことにする。）
この定義に従えば、例えば、2次元で均等分布が保証されている場合には、1の後に何が出るのか全くわからない（特定の数字の後の数字はどれも等しい数で出現する）ですが、2次元で均等分布していない場合には、連続する2数の組が均等に出現しないため、例えば1の次には100が出やすい、といったことが起きます。
高次元で均等分布していない例につきましては、SFC版風来のシレンの乱数生成アルゴリズムの話 考察編の最後のグラフが非常に参考になります。
LCGについては次のことがわかります。
命題（LCGの高次元分布）
次が成り立つ。
- 周期$m$のLCGは1次元に均等分布する。
- LCGは2次元以上に均等分布しない。
（証明）
1.
上の命題（一階漸化式が生成する擬似乱数列の周期）の証明から、周期$m$のLCGでは$x_0$から$x_{m-1}$の数字は、0からm-1までの相異なる数字であることが分かります。このとき、明らかに全ての数字は一度ずつ出現します。したがって一周期にわたってとったときに$X = \{ 0, 1, \cdots ,m-1\}$の全ての点が一度ずつ出現するため、周期$m$のLCGは1次元に均等分布しています。
$X^2$の要素の数は$m^2$です。そのため、$(0, 0)$という組を除く$X^2$の要素が同じ回数ずつ出現するには周期は少なくとも$m^2-1$でなければなりません。
しかし、LCGの周期は$m$を超えないので、LCGは2次元に均等分布しません。
（LCGでは特定の数字の後には同じ数字しか出現しないことからも明らかです。）
（証明終）
LCGでは2次元の均等分布をしておらず、複数個の乱数を組にして使用するような用途には全く向いていません。
複数個の乱数を用いる場合には互いに異なるseed値を用いる必要がありますが、LCGで最大周期で良い乱数性を持つようなseedがそこまで用意できないという事情も相まって、LCGをそのまま用いることは特に近年ではあまりないと思われます。
ちなみにLCGが最大周期を持つ条件はOn the Period Length of Pseudorandom Number Sequencesなどを参照していただければ書いてありますが、本旨とは外れる上に証明もかなり面倒なので割愛します。
線形帰還シフトレジスタ（LFSR）
ここまででLCGはあまりよくない乱数だということが分かりましたが、では、より良い乱数を作るためにはどうすればいいでしょうか？
一階の漸化式が生成する乱数列の場合にはどうあがいても周期はそこまで伸びないことを見たので、少なくともより高階の漸化式を用いる必要があるでしょう。
ついでに$\mathrm{~mod~}$をとっていた$m$も計算を簡単にするために$m=2$としておくと、次のような生成法が考えられます。
\begin{align}
x_{i+2} = (x_{i+1} + x_i) \mathrm{~mod~} 2
\end{align}
このように、一階以上の線形な漸化式を用いて$2$でmodをとって生成する乱数列を線形帰還シフトレジスタ（Linear feedback shift register, LFSR）といいます。
次の乱数の値を得るために、一つ前の数字ともう一つ前の数字を足していて、Fibonacci数のような雰囲気をもった擬似乱数列です。
具体的に$x_0 = 0, x_1 = 1$として実際に数列を見てみると、
\begin{align}
x_0 =& 0\\
x_1 =& 1\\
x_2 =& (0 + 1) \mathrm{~mod~} 2 = 1\\
x_3 =& (1 + 1) \mathrm{~mod~} 2 = 0\\
x_4 =& (1 + 0) \mathrm{~mod~} 2 = 1\\
x_5 =& (0 + 1) \mathrm{~mod~} 2 = 1\\
\vdots &
\end{align}
となり、どうやら(0, 1, 1)という系列を繰り返し、周期は3となっているようです。
$m=2$から一階漸化式を用いた場合の最長周期は2であるので、一階漸化式を用いた場合よりも周期は伸びています。
また、高次元分布についても確認しておきましょう。
まず、1次元の均等分布は一周期にわたって0と1の個数を数えることで確認できます。
0が1つ、1が2つとなっていて、0の個数は他の数より一つ少なくても良いという定義でしたので、1次元には均等分布しているようです。
2次元の均等分布は、一周期にわたって隣接する2つの数を組みにしたもの、すなわち（$[0,1], [1,1], [1,0]$）から確認します。$[0, 0]$という組のみ出現していませんが、それ以外の0と1の組み合わせは全て一度ずつ出現していて、2次元で均等分布していることが確認できます。（[$0, 0]$のみは他よりも一度出現回数が少なくても良いのでした。）
さて、今は乱数を生成するために2つ前の値を参照する、つまり二階の漸化式を用いて乱数を生成していたのでした。
このとき、周期は$3 = 2^2 - 1$まで増加し、また1ビットが2次元に均等分布するようになりました。
実は、$n$階の漸化式を用いたLFSRについて次が言えます。
命題（$n$階の漸化式を用いたLFSRの周期と高次元分布）
$n$階の漸化式を用いたLFSRを
\begin{align}
&x_{i+n} = (a_{n-1} x_{i+(n-1)} + a_{n-2} x_{i+(n-2)} + \cdots + a_0 x_i) \mathrm{~mod~} 2\\
&(a_{n-1}, a_{n-2}, \cdots , a_1 \in \{ 0, 1\}, a_0 = 1)
\end{align}
と定義する。（ただし$n\ge 1$）このとき、次が成り立つ。
- 周期は$2^n-1$を超えない。
- 1ビット精度で$n+1$次元以上に均等分布しない。
（証明）
1.
まず、$x_{i+(n-1)} = x_{i+(n-2)} = \cdots = x_i = 0$となったとき、その後の乱数の値は全て0となることに注意すると、もしもそのようなパターンが一度でも出現した場合には周期が1となることが分かります。
対偶を取れば、周期が1よりも大きい場合には$x_i$から$x_{i+(n-1)}$まで全て0というパターンは決して出現しません。
$n$個の乱数を組にしたもの$[x_{i+(n-1)}, x_{i+(n-2)}, \cdots , x_{i}]$を考えると、全て0というパターンは除いて、取り得るパターンは$2^n - 1$個です。したがって、ある$0 \le i < j \le 2^n-1$で乱数の組が全く同じになるもの、すなわち
$$
[x_{i+(n-1)}, x_{i+(n-2)}, \cdots , x_{i}] = [x_{j+(n-1)}, x_{j+(n-2)}, \cdots , x_{j}]
$$
となるような$i, j$が存在します。以降は上の命題（一回漸化式が生成する乱数列の周期）の証明と同様に考えて、周期が$2^n - 1$を超えないことが証明できます。
1ビット精度で$n+1$次元に均等分布する場合、少なくとも$2^{n+1}-1$のパターンが1度ずつ出現する必要があるので、周期は$2^{n+1}-1$よりも長い必要があります。1から周期は$2^n-1$を超えないことが示されたので、$n+1$次元に均等分布しません。より大きな次元に関しても同様に示すことができます。
（証明終）
この命題から、高階の漸化式を用いればより長周期かつ高次元で均等分布するような乱数列を生成できそうです。
実際により高階の漸化式を用いた乱数列を生成してみましょう。次は、三階の漸化式
$$
x_{i+3} = x_{i+2} + x_{i} \mathrm{~mod~} 2
$$
を用いて乱数を生成します。（$x_{i+3}$を生成するために3つ前の乱数の値$x_i$を参照しているので三階の漸化式）
$x_0 = 1, x_1 = 1, x_2 = 1$として乱数列を生成すると、
\begin{align}
x_0 =& 1\\
x_1 =& 1\\
x_2 =& 1\\
x_3 =& (1+1) \mathrm{~mod~} 2 = 0\\
x_4 =& (0+1) \mathrm{~mod~} 2 = 1\\
x_5 =& (1+1) \mathrm{~mod~} 2 = 0\\
x_6 =& (0+0) \mathrm{~mod~} 2 = 0\\
x_7 =& (0+1) \mathrm{~mod~} 2 = 1\\
x_8 =& (1+0) \mathrm{~mod~} 2 = 1\\
x_9 =& (1+0) \mathrm{~mod~} 2 = 1\\
x_{10} =& (1+1) \mathrm{~mod~} 2 = 0\\
&\vdots
\end{align}
となり、どうやら周期は7で(1, 1, 1, 0, 1, 0, 0)というパターンの繰り返しのようです。周期は先ほどの命題での上限だった$7 = 2^3 - 1$となっており、また、3次元までの全ての次元で均等分布しているため、こちらも上限に達していることが分かります。
この調子で四階の漸化式
$$
x_{i+4} = x_{i+2} + x_{i} \mathrm{~mod~} 2
$$
でも作ってみます。実際に計算すると、
\begin{align}
x_0 =& 1\\
x_1 =& 1\\
x_2 =& 1\\
x_3 =& 1\\
x_4 =& (1+1) \mathrm{~mod~} 2 = 0\\
x_5 =& (1+1) \mathrm{~mod~} 2 = 0\\
x_6 =& (0+1) \mathrm{~mod~} 2 = 1\\
x_7 =& (0+1) \mathrm{~mod~} 2 = 1\\
x_8 =& (1+0) \mathrm{~mod~} 2 = 1\\
x_9 =& (1+0) \mathrm{~mod~} 2 = 1\\
x_{10} =& (1+1) \mathrm{~mod~} 2 = 0\\
&\vdots
\end{align}
となり周期は6となっています。残念ながら先ほどの三階の漸化式を用いた場合よりも周期が短くなってしまいました。
先ほどの命題では周期や高次元分布の上限を与えましたが、どんな$n$階の漸化式でも上限を達成するわけではなく、上限を達成するには何かの条件を満たす必要ががありそうです。
これからそのための条件を探りにいくのですが、一つの記事に書くには少し長くなるので、次回に回します。
まとめ
2回に分けてMTについてちょっとだけ真面目に扱おうという記事ですが、今回は準備としてLCGとLFSRについて扱いました。
LCGは一階の漸化式を用いて乱数列を生成する方法で、その周期はmodをとっている数$m$を超えず、また、2次元以上には均等分布しないことを確認しました。
このような悪い性質を示すのはLCGが一階の漸化式を用いることだったので、より高階の漸化式を用いるLFSRを用いたところ、周期は最長で$2^n-1$、1ビット精度で最大$n$次元までの次元に均等分布できることが分かりました。
しかし、$n$階の漸化式を用いた場合には必ずこの上限を達成できるわけではなく、最大値を達成するにはさらに条件が必要です。