📝
【技術】擬似乱数への道のり(1)
はじめに
ある仕事の担当をランダムに決めたい場合やプログラムに対するテストケースを無作為に作りたい場合など、ランダム性が求められることはよくあります。しかし、コンピュータは決定的な動作をするのが得意であるため、真の乱数を生成することは困難です[1]。そこで、厳密な乱数を扱うのは諦めて、一見乱数のように見える擬似乱数というものを扱います。擬似乱数生成器は決まった数式と初期値（シード）をもちます。数式と初期値にしたがって擬似乱数が生成されていくので、擬似乱数生成器 , 初期値 として、 番目の擬似乱数 は と表せます。本稿では、擬似乱数に求められる性質を考えながら、いくつかの擬似乱数生成アルゴリズムを紹介します。
乱択アルゴリズム
あるアルゴリズムにおいて入力が決まれば常に同じ結果を出力するとき、そのアルゴリズムを決定的といいます。これとは対照的に、同じ入力であっても実行するたびに異なる値になる可能性があるとき、そのアルゴリズムを乱択アルゴリズムといいます[1]。
身近な乱択アルゴリズムの例として、円周率を求めるモンテカルロ法があります。1 辺の長さが 2 の正方形の中に半径の長さが 1 の円を書きます。正方形の中にでたらめに沢山点をうち、円の中に入った個数を数えます。打った点の総数と円の中に入った点の数との比から円周率を求める方法です[1]。
1 辺の長さが 2 の正方形を 1/4 にして、ランダムに 20 個の点を打つと下図のようになりました。
円内には 16 個の点が入っているので、
のように円周率を（大雑把に）求めることができました。
円周率を求めるモンテカルロ法の他にも、素数判定を効率的に行うミラーラビン素数判定法[2] など様々な乱択アルゴリズムが存在します。
乱数の有用性を確認したところで、実際に擬似乱数を生成するアルゴリズムを紹介していきます。
平方採中法
初期に提案された擬似乱数に、平方採中法（平方中抜き法, middle-square method）があります[3]。平方採中法では、 桁のシード に対し、 の上位 桁と下位 桁を削除したものを次の擬似乱数とします。すなわち、平方採中法によって得られる 進数 桁の擬似乱数列を のように書くことにすると、
と表せます。ただし、（ に対し、） は を で割った余りとします。 を で割った余りとは、 を満たす を指します。
例えば、10 進数 4 桁で とすると、
となり、なんだか良さげな列が生成されているように見えます。しかし、平方採中法には「0 が続く列になってしまいやすい」などの問題点があります[4]。実際、先ほど示した系列も となり、以降 0 が続く列となってしまいます。このような欠点があるため、現在では平方採中法が利用されることはほぼありません。
擬似乱数の周期
多くの擬似乱数生成器は内部状態を持ち、内部状態にしたがって次の擬似乱数を生成します。内部状態はコンピュータ上に置かれるためそのサイズは有限であり、その種類数は有限です。したがって、擬似乱数を生成し続ければいつか同じ内部状態となります。このことから、内部状態を持つ擬似乱数生成器が生成する擬似乱数列はある周期をもつことがわかります。ただし、擬似乱数列 において、 となる が存在し、またそのような のうち最小であるとき、この擬似乱数列は周期 を持つということにします。例えば、平方採中法で、10進数4桁で とすると、
となります。平方採中法では は のみによって定まることを思えば、 のループとなっており、周期が 4 であることがわかります。0 が続く列も周期 1 と捉えることができます。このように周期が短いのは擬似乱数として好ましくありません。
次に紹介する線形合同法はこの周期に関してかなり良い性質を持っています。
線形合同法
整数の3つ組 が を満たすとします。このとき、初期値
として， を次の漸化式で定めます。
線形合同法（linear congruential generators）はこの漸化式にしたがって、擬似乱数を生成するアルゴリズムです[5][6]。
例として、 の場合を考えてみると、
のようになります。
線形合同法で生成される擬似乱数 は を明らかに満たします。したがって、周期は 以下です。実は、線形合同法について周期 を達成するための必要十分条件が知られていて、次のようになります[6]。
定理（線形合同法の最大周期）
パラメータ の線形合同法の周期を とする。このとき、次の必要十分条件が成り立つ。ただし、 は が互いに素であることを意味し、 は が を割り切ることを意味するものとする。
これらを満たす例としては などが挙げられます。初期値 として、擬似乱数列を計算してみると、
となり、確かに周期 8 であることがわかります。
が与えられたとき、最大周期になるように を設定するのは（可能であるなら、計算量を無視すれば）簡単で、
のようにすれば条件を満たす最小の数 が得られるので、これが 未満であれば、 なる とともに用いればよいです。 ただし、 が偶数で最大周期を達成しようとすると偶数奇数が順に出るようになってしまうので、注意が必要です。
が偶数のとき、条件(i)(ii) から は奇数。偶数で割った余りを取るということは何らかの偶数を引くことであることに注意すると、
となるため。
大抵のプログラミング言語に用意されている rand() 関数は線形合同法を使って実装されることがあります[1]。例えば、C++ の rand() 関数は線形合同法で実装されており、その周期は 以下です[7]。
パラメータ の線形合同法によって生成される擬似乱数 は の 通りの値のみを取りえます。したがって、周期が のとき擬似乱数は 上に一様分布します。この「取りうる値上に一様分布する」という性質は擬似乱数にとって重要な性質です。パラメータによっては一様分布する線形合同法は優れた擬似乱数生成器のように思えます。しかし、連続するいくつかの擬似乱数を取ってきて多次元で考えると、線形合同法では一様分布しなくなってしまいます[7]。
次に紹介するメルセンヌ・ツイスタは という長大な周期と、623 次元という高次元に一様分布するという性質を持つ擬似乱数生成器です[8]。
メルセンヌ・ツイスタ
まず、 次元に一様分布するとはどういうことかを確認します。[8] では、これは -distribution として次のように定義されています。
定義 (-distribution)
周期 の擬似乱数生成器が生成する擬似乱数列を とする。また、 を ビット擬似乱数 の先頭 ビットを取ってきたものとする。 について， ビットベクトル
を作る． ビットベクトルは の値を取りうるが、そのどれもが同じ回数だけ上記の 個のベクトルに登場する（ただし、 は他のベクトルより登場回数が1少ない）とき、擬似乱数列 は -distribution を ビットで満足するという。また、ある について，満たす -distribution の最大値を とし、 を調べることを -distribution test と呼ぶ。
この定義は次のように解釈できます[8]。
ビットベクトルを数値とみなしたもの を で割ることで，実数 に対応させたものを と書くことにします。ここで、ベクトル は、 次元単位超立方体内部の点となります。ベクトル は， 次元単位超立方体の各軸を 個に分割してできる格子点 のいずれかにあたります。-distribution を ビットで満足するとき、周期 個のベクトルが全て同じ回数（ は 1 回少ない）だけ各格子点上に位置するということを意味しています。
-distribution を ビットで満足するときは、下の図に示すような格子点に一様分布するということになります。
メルセンヌ・ツイスタの構成
次元に一様分布することの意味が確認できたので、ここからはメルセンヌ・ツイスタ（mersenne twister）について説明していきます。メルセンヌ・ツイスタは以下のように構成されます[8][9][10]。
ビットの整数からなる乱数列を生成する場合を考えます。また、整数を各ビットで分けて 次元行ベクトル として考えます。ここで、メルセンヌ・ツイスタの式は、
となります（）。ただし、は の上位 ビットと の下位 ビットを連結したベクトルであり、 は次に示す 正方行列
です。また、 の計算は有限体 で行います。すなわち、足し算はXOR 演算であり、掛け算はAND 演算です。 の計算は、愚直に行うと 回ほどの演算が必要になりますが、 の構造を利用して次のように高速に計算することが可能です。
ただし、
です。
ここで、 の分布を改善するために， の正則行列 を右から掛けます。この工程はTempering と呼ばれます。ここでも、 を愚直に掛けることはせず、 を掛けることに相当する演算として、パラメータ に対して以下を行います。
例えば、 は に
という形式の行列を右から掛けることに対応します。他の式に関しても、左下の 1 の列が右上に移動したり、この 1 の列に 0 が混じるような列となります[11]。各行列は明らかに rank が で正則なので、 は正則です。こうして得られる を擬似乱数として出力します（内部状態は で持っておくことに注意してください）。
メルセンヌ・ツイスタの仕組みは上に説明した通りです。ここで、
とパラメータを設定すると、周期が で 623-distribution を 32 ビットで満足するようになります。このパラメータが設定されたメルセンヌ・ツイスタは MT19937 と呼ばれます[8][10]。
詳細と性質
さて、 を計算するには、 が必要でした。したがって、 から までを状態として保存しておく必要があります。これは、長さ の配列で、 の位置にヘッドがくるようにすればよいです。
また、初期状態 をどうするかという問題があります。実は、 および の上位 ビットが全て である場合を除けば、どの初期状態からでも 通りの内部状態全てを巡ります[8]。このことは、次のように説明できます。
MT19937 の内部状態のうち以降生成される擬似乱数に関わるビットは と の上位 ビットの合計 ビットです。 ビットの列は 通り存在します。ここで、 なので、 ビットが全て である場合、内部状態はずっと ビットの であることがわかります。したがって、周期 は以降生成される擬似乱数に関わる ビットについて、全て である場合を除いてありうる全ての状態を巡ることを意味しています。
上の議論から周期 というのは、 ビットの内部状態に対して殆ど限界といえる性能をしていることがわかります。このことは623次元に一様分布するという性質についても同様です。
と周期 に関して、 という不等式が成り立ちます[7]。もし だと、 個のベクトルで ビットベクトルの全て（ の登場回数は 1 少ない）をカバーできないからです。 ここで、MT19937について で考えてみると、 なので、 となり、 がわかります。
ここまで説明したように、メルセンヌ・ツイスタは非常に優れた性質を持っています。しかし、連続する 624 個の擬似乱数を取得すれば、内部状態が再現できてしまい以降の擬似乱数が予測できてしまいます[11]。これは、Tempering に用いる が正則であることから内部状態 が取得できるからです。
暗号論的擬似乱数へ
暗号文は決定的に生成されるべきでなく、暗号は乱択アルゴリズムであることが望まれます[5]。
しかし、ここまで紹介してきた擬似乱数生成器のように次の乱数が予測可能であることは、暗号に利用する上で問題となります[1]。例えば、共通鍵の推測による解読や、初期化ベクトルの推測による暗号化プロセスのシミュレートにつながります。
このように、暗号の用途で使われる乱数はシミュレーションなどで使われる擬似乱数とは異な
る性質が要求されるため、通常の擬似乱数と区別するために暗号論的擬似乱数と呼ばれます（暗号に利用する前提なら単に擬似乱数ということもあります）[1]。
次回は、暗号論的擬似乱数に求められる性質やその構成について解説します。
まとめ
- コンピュータで乱数を扱うため、擬似乱数というものを考える
- 擬似乱数を生成するアルゴリズムのなかには、平方採中法、線形合同法、メルセンヌ・ツイスタがある。また、それらの仕組みについて説明した
- 擬似乱数にとって、周期と分布は重要な性質である
- ここで紹介した擬似乱数は暗号利用されるべきでなく、暗号利用するには暗号論的擬似乱数を用いるべき
参考文献
- 光成滋生, “クラウドを支えるこれからの暗号技術”, 2021. https://herumi.github.io/ango/
- Rabin, Michael O. "Probabilistic algorithm for testing primality." Journal of number theory 12.1 (1980): 128-138. https://www.sciencedirect.com/science/article/pii/0022314X80900840
- J. Von Neumann, “Various techniques used in connection with random digits,” Applied Math Series, Notes by G. E. Forsythe, in National Bureau of Standards, Vol. 12, pp. 36-38, 1951.
- 二乗中抜き法の欠点とTAOCP https://qiita.com/saka1_p/items/30ad1770d26e60e630b4
- Park, Stephen K., and Keith W. Miller. "Random number generators: good ones are hard to find." Communications of the ACM 31.10 (1988): 1192-1201. https://dl.acm.org/doi/pdf/10.1145/63039.63042
- 線形合同法（擬似乱数生成法）の周期, tsujimotterのノートブック. https://tsujimotter.hatenablog.com/entry/period-of-linear-congruential-generators
- 松本眞, “あなたの使っている乱数、大丈夫？-危ない標準乱数と、メルセンヌ・ツイスター開発秘話-,” 第50回市村学術賞記念 先端技術講演会(2014) http://www.math.sci.hiroshima-u.ac.jp/m-mat/TEACH/ichimura-sho-koen.pdf
- Matsumoto Makoto, and Takuji Nishimura. "Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator." ACM Transactions on Modeling and Computer Simulation (TOMACS) 8.1 (1998): 3-30. https://dl.acm.org/doi/pdf/10.1145/272991.272995
- Mersenne Twister with improved initialization (2002) http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/MT2002/mt19937ar.html
- メルセンヌ・ツイスタをわかった気になる https://6715.jp/posts/5/
- メルセンヌ・ツイスタを倒す https://6715.jp/posts/6/