Implement Park-Miller Random Number Generation Algorithm in C++

The Park-Miller algorithm is a type of linear congruential generator (LCG) that is used to generate pseudo-random numbers. It is also called the minimal standard random number generator. The general formula of a random number generator (RNG) of this type is: X_{k+1} = g X(k) mod n where the modulus n is a prime number or a power of a prime number, the multiplier g is an element of high multiplicative order modulo n, and the seed X0 is co-prime to n.

In this article, we will be using Park-Miller algorithm to generate 5 pseudo-random numbers using C++. Here is the formula of the Park-Miller algorithm:

X(n+1) = (a * Xn) mod m

where,
Xn is the current number (seed)
X(n+1) is the next random number
a = 16807, prime multiplier
m = 2^31 - 1 = 2147483647, large prime modulus

Steps to Implement Park-Miller Algorithm

We have implemented the Park-Miller algorithm to generate random numbers using the following steps:

• At the beginning, we have defined various variables (a , m , q , and r) that is required in the formula of the park-miler algorithm.
• The variable a represents the multiplier, m is the modulus, q is m / a, and r is m % a. The q and r are used to avoid the overflow.
• We provide an initial seed value for generating the sequence of random numbers.
• We have used the func() function to generate a new random number using the LCG formula.
• Inside func(), we calculate the value of hi and lo to get the upper and lower bound that is later used to find the value of t.
• We calculate t to avoid overflow. We use the if/else statement to check if t is positive or negative.
• If t is positive, it becomes the new seed, otherwise, m is added to keep it within bounds. Then the random number is returned in the range of [0, 1) by dividing the seed by m.

C++ Program to Implement Park-Miller Algorithm

Below is the C++ program to implement the Park-Miller random number generation algorithm:

#include <iostream>
using namespace std;

const long a = 16807;
const long m = 2147483647;
const long q = 127773;     // m / a
const long r = 2836;       // m % a
long seed = 12345678L;

double func() {
    long hi = seed / q;
    long lo = seed % q;
    long t = a * lo - r * hi;
    if (t > 0)
        seed = t;
    else
        seed = t + m;
    return static_cast<double>(seed) / m;
}

int main() {
    cout << "Random numbers are:\n";
    for (int i = 0; i < 5; i++)
        cout << func() << endl;
    return 0;
}

The output of the above code is:

Random numbers are:
0.621835
0.177248
0.00290613
0.843384
0.76204