Random Number Generation via Linear Congruential Generators in C++

Random Number Generation via Linear Congruential Generators in C++

In this article we are going to construct classes to help us encapsulate the generation of random numbers. Random number generators (RNG) are an essential tool in quantitative finance as they are necessary for Monte Carlo simulations that power numerical option pricing techniques. Other articles on QuantStart have so far used RNGs in a procedural manner. In particular, we have utilised the Box-Muller technique to generate one or more random variables distributed as a standard Gaussian.

We will now show how to construct a random number generator class hierarchy. This allows us to separate the generation of random numbers from the Monte Carlo solvers that make use of them. It helps us reduce the amount of code we will need to write in the future, increases extensibility by allowing easy creation of additional random number generators and makes the code more maintainable.

There are further reasons to write our own random number generators:

  • It allows us to make use of pseudo-random numbers. These are sequences of numbers that possess the correct statistical properties to "emulate" random numbers in order to improve the convergence rates of Monte Carlo simulations. The interface for random numbers and pseudo-random numbers is identical and we can hide away the details in the specific classes. In particular we can implement low-discrepancy numbers and anti-thetic sampling in this manner.
  • Relying on the rand function provided with the C++ standard is unreliable. Not only is rand implementation specific, because it varies across multiple vendor compilers, but we are unaware of the efficiency of each implementation. This leads to difficulties in cross-platform testing as we cannot guarantee reproducibility.
  • We are able to provide multiple separate streams of random numbers for different parts of our running program. The seed for the rand function, srand, is a global variable and hence will affect all components of our program, which is usually unwanted behaviour. By implementing our own RNG we avoid this issue.

Random Number Generator Abstract Base Class

Our random number generators will be formed from an inheritance hierarchy. To form the hierarchy we will create an abstract base class that specifies the interface to the random number generator. All subsequent generators will inherit the interface from this class.

The primary considerations of this interface are as follows:

  • Quantity or dimension of the generator: Many of the options pricers we have already created require more than a single random number in order to be accurately priced. This is the case for path-dependent options such as Asians, Barriers and Lookbacks. Thus our first consideration is to make sure that the generator provides a vector of random numbers, with the dimension specified at the creation of the instance.
  • The supported statistical distributions from which to draw random variables: For options pricing, the two main statistical distributions of interest will be the uniform distribution and the standard normal distribution (i.e. the "Gaussian" distribution). Gaussian random draws are calculated from uniform random draws. We can use the statistical classes we created in the previous article in order to obtain random draws from any particular distribution we wish, without modifying the RNG.
  • We will need methods to support obtaining and setting the random seed, so that we can control which random numbers are generated and to ensure reproducibility across separate runs and platforms.

With those considerations in mind, let's create a simple abstract base class for our random number generator, in the file random.h:

#ifndef __RANDOM_H
#define __RANDOM_H

#include <vector>

class RandomNumberGenerator {
 protected:
  unsigned long init_seed;  // Initial random seed value
  unsigned long cur_seed;   // Current random seed value
  unsigned long num_draws;  // Dimensionality of the RNG
  
 public:
  RandomNumberGenerator(unsigned long _num_draws, unsigned long _init_seed) 
    : num_draws(_num_draws), init_seed(_init_seed), cur_seed(_init_seed) {};
  virtual ~RandomNumberGenerator() {};

  virtual unsigned long get_random_seed() const { return cur_seed; }
  virtual void set_random_seed(unsigned long _seed) { cur_seed = _seed; }
  virtual void reset_random_seed() { cur_seed = init_seed; }
  virtual void set_num_draws(unsigned long _num_draws) { num_draws = _num_draws; }

  // Obtain a random integer (needed for creating random uniforms)
  virtual unsigned long get_random_integer() = 0;
  
  // Fills a vector with uniform random variables on the open interval (0,1)
  virtual void get_uniform_draws(std::vector<double>& draws) = 0; 
};

#endif

Let's run through the code. Firstly, note that we have three protected member variables (which are all large unsigned long integers). cur_seed is the RNG current seed value. init_seed is the initial seed value, which does not change once the RNG has been instantiated. The current seed can only be reset to the initial seed. num_draws represents the dimensionality of the random number generator (i.e. how many random draws to create):

protected:
 unsigned long init_seed;  // Initial random seed value
 unsigned long cur_seed;   // Current random seed value
 unsigned long num_draws;  // Dimensionality of the RNG

Since we're creating an abstract base class is it a good idea to use protected data?

This is actually a contentious issue. Sometimes protected variables are frowned upon. Instead, it is argued that all data should be private and that accessor methods should be used. However, inherited classes -are- clients of the base class, just as "public" clients of the classes are. The alternative argument is that it is extremely convenient to use protected member data because it reduces the amount of cluttered accessor and modifier methods. For brevity I have used protected member data here.

Although the class will never be instantiated directly, it still has a constructor which must be called to populate the protected members. We use a member initialisation list to carry this out. We also create an empty method implementation for the constructor ({}), avoiding the need to create a random.cpp source file. Notice that we're setting the current seed to the initial seed as well.

RandomNumberGenerator(unsigned long _num_draws, unsigned long _init_seed) 
    : num_draws(_num_draws), init_seed(_init_seed), cur_seed(_init_seed) {};

We then have four separate access and reset methods (all virtual), which get, set and reset the random seed and another which resets the number of random draws. They are all directly implemented in the header file, once again stopping us from needing to create a random.cpp source file:

virtual unsigned long get_random_seed() const { return cur_seed; }
virtual void set_random_seed(unsigned long _seed) { cur_seed = _seed; }
virtual void reset_random_seed() { cur_seed = init_seed; }
virtual void set_num_draws(unsigned long _num_draws) { num_draws = _num_draws; }

We now need a method to create a random integer. This is because subsequent random number generators will rely on transforming random unsigned longs into uniform variables on the open interval $(0,1)$. The method is declared pure virtual as different RNGs will implement this differently. We don't want to "force" an approach on future clients of our code:

// Obtain a random integer (needed for creating random uniforms)
virtual unsigned long get_random_integer() = 0;

Finally we fill a supplied vector with uniform random draws. This vector will then be passed to a statistical distribution class in order to generate random draws from any chosen distribution that we implement. In this way we are completely separating the generation of the uniform random variables (on the open interval $(0,1)$) and the draws from various statistical distributions. This maximises code re-use and aids testing:

// Fills a vector with uniform random variables on the open interval (0,1)
virtual void get_uniform_draws(std::vector<double>& draws) = 0;

Our next task is to implement a linear congruential generator algorithm as a means for creating our uniform random draws.

Linear Congruential Generator Implementation

Linear congruential generators (LCG) are a form of random number generator based on the following general recurrence relation:

\begin{eqnarray*} x_{k+1} = g \cdot x_k \mod n \end{eqnarray*}

Where $n$ is a prime number (or power of a prime number), $g$ has high multiplicative order modulo $n$ and $x_0$ (the initial seed) is co-prime to $n$. Essentially, if $g$ is chosen correctly, all integers from 1 to $n-1$ will eventually appear in a periodic fashion. This is why LCGs are termed pseudo-random. Although they possess "enough" randomness for our needs (as $n$ can be large), they are far from truly random. We won't dwell on the details of the mathematics behind LCGs, as we will not be making strong use of them going forward in our studies. However, most system-supplied RNGs make use of LCGs, so it is worth being aware of the algorithm. The listing below (lin_con_gen.cpp) contains the implementation of the algorithm. If you want to learn more about how LCGs work, take a look at Numerical Recipes[1].

With the mathematical algorithm described, it is straightforward to create the header file listing (lin_con_gen.h) for the Linear Congruential Generator. The LCG simply inherits from the RNG abstract base class, adds a private member variable called max_multiplier (used for pre-computing a specific ratio required in the uniform draw implementation) and implements the two pure virtual methods that were part of the RNG abstract base class:

#ifndef __LINEAR_CONGRUENTIAL_GENERATOR_H
#define __LINEAR_CONGRUENTIAL_GENERATOR_H

#include "random.h"

class LinearCongruentialGenerator : public RandomNumberGenerator {
 private:
  double max_multiplier;
  
 public:
  LinearCongruentialGenerator(unsigned long _num_draws, unsigned long _init_seed = 1);
  virtual ~LinearCongruentialGenerator() {};
  
  virtual unsigned long get_random_integer();
  virtual void get_uniform_draws(std::vector<double> draws); 
};

#endif

The source file (lin_con_gen.cpp) contains the implementation of the linear congruential generator algorithm. We make heavy use of Numerical Recipes in C[1], the famed numerical algorithms cookbook. The book itself is freely available online. I would strongly suggest reading the chapter on random number generator (Chapter 7) as it describes many of the pitfalls with using a basic linear congruential generator, which I do not have time to elucidate on in this article. Here is the listing in full:

#ifndef __LINEAR_CONGRUENTIAL_GENERATOR_CPP
#define __LINEAR_CONGRUENTIAL_GENERATOR_CPP

#include "lin_con_gen.h"

// This uses the Park & Miller algorithm found in "Numerical Recipes"
// See reference [1] for more details

// Define the constants for the Park & Miller algorithm

const unsigned long a = 16807;       // 7^5
const unsigned long m = 2147483647;  // 2^32 - 1 (and thus prime)

// Schrage's algorithm constants

const unsigned long q = 127773;
const unsigned long r = 2836;

// Parameter constructor
LinearCongruentialGenerator::LinearCongruentialGenerator(unsigned long _num_draws, 
                                                         unsigned long _init_seed) :
  RandomNumberGenerator(_num_draws, _init_seed) {

  if (_init_seed == 0) {
    init_seed = 1;
    cur_seed = 1;
  }

  max_multiplier = 1.0 / (1.0 + (m-1));
}

// Obtains a random unsigned long integer
unsigned long LinearCongruentialGenerator::get_random_integer() {
  unsigned long k = 0;
  k = cur_seed / q;
  cur_seed = a * (cur_seed - k * q) - r * k;
  
  if (cur_seed < 0) {
    cur_seed += m;
  }

  return cur_seed;
}

// Create a vector of uniform draws between (0,1)
void LinearCongruentialGenerator::get_uniform_draws(std::vector<double>& draws) {
  for (unsigned long i=0; i<num_draws; i++) {
    draws[i] = get_random_integer() * max_multiplier;
  }
}

#endif

Firstly, we set all of the necessary constants (see [1] for the explanation of the chosen values). Note that if we created another LCG we could inherit from the RNG base class and use different constants:

// Define the constants for the Park & Miller algorithm

const unsigned long a = 16807;       // 7^5
const unsigned long m = 2147483647;  // 2^32 - 1

// Schrage's algorithm constants

const unsigned long q = 127773;
const unsigned long r = 2836;

Secondly we implement the constructor for the LCG. If the seed is set to zero by the client, we set it to unity, as the LCG algorithm does not work with a seed of zero. The max_mutliplier is a pre-computed scaling factor necessary for converting a random unsigned long into a uniform value on on the open interval $(0,1) \subset \mathbb{R}$:

// Parameter constructor
LinearCongruentialGenerator::LinearCongruentialGenerator(unsigned long _num_draws, 
                                                         unsigned long _init_seed) :
  RandomNumberGenerator(_num_draws, _init_seed) {

  if (_init_seed == 0) {
    init_seed = 1;
    cur_seed = 1;
  }

  max_multiplier = 1.0 / (1.0 + (m-1));
}

We now concretely implement the two pure virtual functions of the RNG base class, namely get_random_integer and get_uniform_draws. get_random_integer applies the LCG modulus algorithm and modifies the current seed (as described in the algorithm above):

// Obtains a random unsigned long integer
unsigned long LinearCongruentialGenerator::get_random_integer() {
  unsigned long k = 0;
  k = cur_seed / q;
  cur_seed = a * (cur_seed - k * q) - r * k;
  
  if (cur_seed < 0) {
    cur_seed += m;
  }

  return cur_seed;
}

get_uniform_draws takes in a vector of the correct length (num_draws) and loops over it converting random integers generated by the LCG into uniform random variables on the interval $(0,1)$:

// Create a vector of uniform draws between (0,1)
void LinearCongruentialGenerator::get_uniform_draws(std::vector<double>& draws) {
  for (unsigned long i=0; i<num_draws; i++) {
    draws[i] = get_random_integer() * max_multiplier;
  }
}

The concludes the implementation of the linear congruential generator. The final component is to tie it all together with a main.cpp program.

Implementation of the Main Program

Because we have already carried out most of the handwork in random.h, lin_con_gen.h, lin_con_gen.cpp, the main implementation (main.cpp) is straightforward:

#include <iostream>
#include "lin_con_gen.h"

int main(int argc, char **argv) {
  // Set the initial seed and the dimensionality of the RNG
  unsigned long init_seed = 1;
  unsigned long num_draws = 20;
  std::vector<double> random_draws(num_draws, 0.0);

  // Create the LCG object and create the random uniform draws
  // on the open interval (0,1) 
  LinearCongruentialGenerator lcg(num_draws, init_seed);
  lcg.get_uniform_draws(random_draws);

  // Output the random draws to the console/stdout
  for (unsigned long i=0; i<num_draws; i++) {
    std::cout << random_draws[i] << std::endl;
  }

  return 0;
}

Firstly, we set up the initial seed and the dimensionality of the random number generator. Then we pre-initialise the vector, which will ultimately contain the uniform draws. Then we instantiate the linear congruential generator and pass the random_draws vector into the get_uniform_draws method. Finally we output the uniform variables. The output of the code is as follows:

7.82637e-06
0.131538
0.755605
0.45865
0.532767
0.218959
0.0470446
0.678865
0.679296
0.934693
0.383502
0.519416
0.830965
0.0345721
0.0534616
0.5297
0.671149
0.00769819
0.383416
0.0668422

As can be seen, all of the values lie between $(0,1)$. We are now in a position to combine our statistical classes with the uniform random number generator. We can extend both of these class hierarchies to include random draws from other distributions, as well as create more efficient generators, such as a Mersenne Twister random number generator. Those tasks will be the subject of later articles.

References