Park-Miller乱数生成アルゴリズムをC++で実装する方法

Park-Miller乱数生成アルゴリズムは、擬似乱数列を生成する古典的な手法の一つです。乗算と剰余演算のみを用いて乱数を生成できるシンプルな構造が特徴で、1988年にStephen ParkとKeith Millerによって発表されたことから「Lehmer型乱数生成器」としても知られています。

Park-Millerアルゴリズムの基本式

このタイプの乱数生成器(RNG)の一般式は次のように表されます。

X_k+1 = g × X(k) mod n

この式が正しく機能するためには、各パラメータが以下の条件を満たす必要があります。

• 法(modulus)n:素数、または素数のべき乗であること
• 乗数(multiplier)g:nを法として高い乗法位数を持つ要素であること
• 初期値(シード)X₀:nと互いに素であること

アルゴリズムの手順

以下は、Park-Millerアルゴリズムの処理手順を示した疑似コードです。

開始
　　変数 n、a、b、c、seed を宣言する
　　変数 n、a、b、c、seed の値を読み込む
　　Uniform() を呼び出す
　　変数 hi、lo、t を宣言する
　　hi = seed ÷ b
　　lo = seed − b × hi
　　t = a × lo − c × hi
　　もし(t > 0)ならば
　　　　seed = t
　　そうでなければ
　　　　seed = t + n
　　seed を返す
完了
i = 0 から n まで繰り返す
　　関数 random を呼び出す
完了
終了

C++による実装例

次のコードは、上記のアルゴリズムをC++で実装したものです。uniform() 関数が乱数の生成を担い、main() 関数内で10個の乱数を配列に格納して表示します。

#include <iostream>
using namespace std;

const long n = 2145678965L;
const long a = 763214L;
const long b = 88844L;
const long c = 7766L;
static long seed = 12345678L;

double uniform() {
    long hi = seed / b;
    long lo = seed - b * hi;
    long t = a * lo - c * hi;
    if (t > 0)
        seed = t;
    else
        seed = t + n;
    return seed;
}

int main(int argc, char **argv) {
    double A[10];
    for (int i = 0; i < 10; i++)
        A[i] = uniform();

    cout << "Random numbers are:\n";
    for (int i = 0; i < 10; i++)
        cout << A[i] << endl;
}

実行結果

このプログラムを実行すると、以下のような乱数列が出力されます。

Random numbers are:
6.50293e+10
4.27187e+10
2.1539e+10
4.62058e+10
1.70792e+10
8.24569e+09
5.93381e+10
3.63839e+10
4.81931e+10
8.91007e+09

まとめ

Park-Millerアルゴリズムは、法・乗数・シードの条件さえ満たせば、非常に少ない計算資源で再現性のある乱数列を生成できる点が大きな利点です。ただし、単純な線形合同法と同様に周期性や統計的性質の限界があるため、暗号用途ではなくシミュレーションや学習目的での利用に適しています。実際の標準的なパラメータとしては、n = 2³¹ − 1、g = 16807(MINSTD)が広く知られています。