線形合同法とは

線形合同法(linear congruential generator : LCG)とは、最も古く、また最もよく知られた擬似乱数生成アルゴリズムの一種です。擬似乱数生成のアルゴリズムは数学的背景の理解が難しいことが多々ありますが、このアルゴリズムは背景理論の理解及びプログラミング言語による実装が容易です。処理速度が高速なこともあり、多くの場面で利用されていましたが、欠点もあるため注意が必要です。

乱数列の生成は次の漸化式を用います。

Xn+1 = (A * Xn + C) mod M

A, C, M は定数とし、それぞれ M>A, M>B, A>0, B≥0 でなければなりません。 X0 はSeed値として与えられることとします(とりあえず 0 や 1 でも構いません)。乱数の出現範囲は最大で 0～M となります。

一例

具体的に適当な値を当てはめて乱数列を生成してみましょう。ここでは A=3, C=4, M=11, X0=1 とします。では漸化式に当てはめていくつかの乱数を生成します。

  • X1 = (3 * 1 + 4) mod 11 = 7
  • X2 = (3 * 7 + 4) mod 11 = 3
  • X3 = (3 * 3 + 4) mod 11 = 2
  • X4 = (3 * 2 + 4) mod 11 = 10
  • X5 = (3 * 10 + 4) mod 11 = 1
  • X6 = (3 * 1 + 4) mod 11 = 7

X1 が X6 と等しくなっています。したがって X7 以下は周期的に同じ値(3,2,10,1,7 ...)が出現することになります。

周期

線形合同法によって生成される擬似乱数は最大で周期Mを持つことができます。この最大周期を実現するには A, C, M が以下の条件を満たす必要があります。例えば、A=13、B=5、M=24 の組み合わせが当てはまります。

このとき、一定周期で 0～M が生成されます。

  • B と M が互いに疎
  • A-1が、Mの持つ全ての素因数で割りきれる
  • Mが4の倍数の場合、A-1も4の倍数である

パラメータ(A, C, M)の設定

Parameters in common use

設定されるパラメータにはいくつかの代表的な組み合わせがあります。上の英語版Wikipediaの記事を参考に幾つかの組み合わせを以下の表にまとめておきます。

ただし生成されるビット列を必ずすべてを利用するわけではなく、下位ビットはランダム性に欠けるためです。詳しくは上記URLを参照してください。

A C M
1664525 1013904223 232
22695477 1 232
1103515245 12345 231 - 1
1103515245 12345 231
134775813 1 232
33797 1 232
214013 2531011 232

メリット・デメリット

線形合同法は最小限のメモリ領域しか必要としません。さらに実装も容易で動作も非常に高速です(剰余算はビットシフト等の効率的な式に変換できます)。これは組み込みやゲームなど限られたメモリしか使えない環境では大きなメリットになります。

問題はあまりよい乱数が生成されないということです。高品質な乱数が必要な場面(暗号など)では利用するべきではありません。Mの値によっては下位ビットのランダム性が低くなりがちです。特にMが偶数の場合には最下位ビットが完全に "0" "1" が交互(つまり偶数基数が交互)に出現してします。

さらに線形合同法による乱数列から2次元座標として(r1, r2), (r3, r4) ... のように順番に割り当ててプロットすると、一定の間隔の格子点状に点が並んでしまうというわかりやすい欠点もあります。組み合わせでも使用することはお勧めできません。

線形合同法によって生成される擬似乱数には下位ビットのランダム性に問題があることがわかってます。したがって下位ビットは使わないように上位ビットのみを利用すればある程度問題は解消されます。

サンプルコード

除算で剰余を求めるのではなく、マスクすることで剰余を求めています。最上位ビットの符号値が常に0になるので正の整数値のみ生成される事になります。ちなみにこの擬似乱数では偶数と奇数が交互に生成されてしまいます。

C#

class LCG
{
    int X = 1;  // 1をシード値とします
    const int A = 1664525;
    const int C = 1013904223;
    const int M = 2147483647; // 2^31
    public int Rand()
    {
        X = (X * A + C) & M;
        return X;
    }
}