rand 関数や java.util.Random と結晶構造

この記事は，SYSKEN Advent Calendar 2018 の23日目の記事です（？）．

代打のふぃらっはです．もうすぐクリスマスですね．自分宛てのプレゼントと思って何冊か本を注文しましたが，ひたすら虚しいだけでした．

今日は，線形合同法（Linear congruential generator）の結晶構造についてちょっと語ろうと思います．

1. そもそも線形合同法って？

線形合同法とは，擬似乱数生成器のひとつです．

(1) $\begin{equation*} X_n = a X_{n-1} + c \pmod{M} \end{equation*}$

なる漸化式によって，順次整数を生成させることができます．ここで $M, a, c$ は正整数のパラメータで，それぞれ法，乗数，加数と呼ばれます．ね，めちゃくちゃ簡単でしょ．皆さんがよくご存知だと思われる C 言語の rand() だったり，Java の java.util.Random でも使われています（ rand() については，厳密にはコンパイラによります）．

2. 高次元疎結晶構造

で，話は急に変わりますが，線形合同法はあまり（というか，かなり）評判がよくありません．その理由の一つとして，「高次元になると一様でもなんでもなくなる」ことが上げられます．
$(X_n, X_{n-1})$ なる点を順次グラフに打っていくと，たとえば以下のようになります．

ぜんぜんランダムじゃないやんけ，と思いますよね．このような点列の構造を $n-$ 次元疎結晶構造と呼びます．線形合同法によって得られる点列が結晶構造をなす理由は以下のとおりです：
漸化式(1)の周期が $k$ だとしましょう． $X_n$ は $0, 1, 2, \cdots , k-1$ をそれぞれ1回ずつたどり，重複しません．したがって，最初の点が $(a, b)$ なる座標に乗ると，以降の点は直線 $x = a, y = b$ の上に乗りません．2点目以降についても同様で，どんどん $X_n$ の取りうる範囲が狭まっていきます． $k$ 個目の点に至っては座標の候補がただ一つに絞り込まれてしまいます．．こうして得られる点列は，いわば $k-$ 飛車問題（参考：n-クイーン問題）の解となる盤面をなします． $k-$ 飛車問題の解となる盤面は必ず結晶構造に帰着します．

これは比較的困った問題ですよね．なお困ったことに，結晶構造は次元が上がるに応じてどんどん疎になってゆきます．次節以降で，この問題を解決する「次善の策」について考えてみます．

3. スペクトル検定

線形合同法を利用する以上，先に示したような疎結晶構造が高次元について出てくるのは避けられません．しかし，結晶は結晶でも，パラメータによって疎な場合と密な場合があります．

このプロット結果は， $M = 256, a = 45, c = 1$ として点列を打ったものです．先に示したプロット結果とは異なり，かなり一様っぽくなってることがわかると思います．

実のところ，高次元でどれほど疎な結晶構造を描くかは，パラメータに依存する部分が大きいです．したがって，あるパラメータを与えたときに，それがどれほど「よい」点列を生成するのかについて検定したいというモチベーションが生まれます．ここで用いるのがスペクトル検定です．詳細は省きますが，スペクトル検定を行うと，結晶構造の「精度」が数値的に得られます．

4. RANDU の置き土産

RANDU という，メインフレーム（部屋中を埋めつくすでっかいコンピュータ）の時代から使われてきたパラメータがあります．具体的には， $M = 2^{31}, a = 65539, c=0$ です．こいつから生成される3次元の点列には，ちょっとびっくりするくらい偏りがあります（下図）．

一目でわかると思いますが，これでは乱数として使えませんね．

「線形合同法」みたいなワードでググって色んなサイトを見ると，RANDU の疎結晶構造がたくさん出てきます．

先に述べたスペクトル検定に RANDU のパラメータをかけてみると，3次元では $\mu = \log_2 \nu_3 \simeq 3.441$ なる数値が出てきます．Knuthによれば，3次元で $\log_2 \nu_3 \ge 10$ 程度の値が一つの基準らしいので，RANDU はかなりヤバい方のパラメータに属すると言えるでしょう．

RANDU は相当長い間使われてきました（今でも使っているシーンが散見されるそうです）が，高次元については実用に堪えないクオリティだったと見なせます．

5. 他のパラメータ

しかし，RANDU が酷いからと言っても，線形合同法自体がまるっきりだめというわけではありません．Borland C Compiler の rand() 関数で使われているパラメータ $M = 2^{32}, a = 22695477 , c = 1013904223$ についてスペクトル検定を行うと， $\mu = \log_2 \nu_3 \simeq 10.088$ という比較的良好な値が出てきます．一般に，現在使われているコンパイラ（BCC, GCC）では，少なくとも3次元までなら実用に堪えうる程度に質のよいパラメータが選ばれています．以下の図は GCC で使われているパラメータをもとに発生させた3次元点列のプロット結果です．RANDU の場合とは全く異なり，かなり一様性が高いのがわかると思います．