randmt: MT19937 pseudorandom number generator

Overview

This C/C++ software contains a high-quality pseudorandom number generator, the Mersenne Twister MT19937 by Takuji Nishimura and Makoto Matsumoto, and routines for sampling from several distributions.

Warning

Do not use for cryptographic purposes. Read Internet RFC4086.

License

This software is distributed under the BSD license. See the included file license.txt for details.

Use

To use randmt in a program, only the files randmt.h and randmt.c are needed; all other included files are for documentation and testing purposes. Include randmt.h and call init_randmt_auto() at the beginning of the program to initialize the pseudorandom number generator. The initialization function seeds the generator with the current time so that different numbers are produced on each run of the program.

Uniform random numbers are generated using the rand_unif() function. Samples from several other distributions can be generated as well, see below.

Example:

#include <stdio.h>
#include "randmt.h"

int main(void)
{
  int i;
  init_randmt_auto();
  for (i = 0; i < 20; ++i) {
    printf("%10.8f\n", rand_unif());
  }
  return 0;
}

This program can be compiled with GCC using the included makefile as

make -f makefile.gcc example

or with MSVC,

nmake -f makefile.vc example.exe

Pseudorandom numbers can be generated for several different distributions:

rand_unif(): continuous uniform distribution on the open interval (0,1)
rand_uint32(): 32-bit unsigned integers uniformly in \(\{0,\ldots, 2^{31}-1\}\)
rand_normal(): the standard normal distribution (zero mean, unit variance)
rand_exp(mu): exponential distribution with mean mu
rand_gamma(a, b): Gamma distribution with shape a and scale b
rand_geometric(p): geometric distribution
rand_poisson(mu): Poisson distribution with mean mu

Reentrant Versions

For use in multithreaded applications, reentrant versions of the functions are also included which have the same name suffixed with _r. With these functions the generator state is passed using a randmt_t object, for instance as rand_unif_r(generator) or rand_exp_r(generator, mu).

A randmt_t object represents the state of an MT19937 pseudo-random number generator. Use the following functions to create, initialize, and destroy randmt_t objects:

new_randmt(): create a new randmt_t
free_randmt(generator): free memory associated with a randmt_t
init_randmt_auto_r(generator): initialize randmt_t with time and address
init_randmt_r(generator, seed): initialize randmt_t with seed

History

This code uses the Mersenne Twister MT19937 code by Makoto Matsumoto and Takuji Nishimura (2002/1/26).
Seiji Nishimura (2008) wrote a reentrant modification of MT19937, where a generator struct is passed to sampling functions.
Nicolas Limare wrote a version with Doxygen documentation, simpler interface, and automatic initialization using the time.
In the code documented here, Pascal Getreuer added samplers for normal, exponential, Gamma, geometric, and Poisson distributions and a test program.

Test Program

To verify the distributions of the samplers, a test program is included. Compile the test program using

make -f makefile.gcc test

or with MSVC,

nmake -f makefile.vc test.exe

The test program applies the Kolmogorov–Smirnov and chi-squared tests to verify that the pseudorandom samplers produce the intended distributions. Note that the output of this program is different on each run.

Typical output is shown below:

For each random number generator, we sample N=1000000
values and compare the sample distribution to the theoretical
density function with the Kolmogorov-Smirnov test.

    Sampler              D
    rand_unif()          0.001095
    rand_normal()        0.001077
    rand_exp(1)          0.000811
    rand_gamma(0.2,1)    0.000516
    rand_gamma(  1,1)    0.000906
    rand_gamma(  2,1)    0.001051
    rand_gamma( 20,1)    0.001079

Supposing the distributions are correct, the D values should be
small with high probability:
    D < 0.001358 with probability 0.95
    D < 0.001627 with probability 0.99
    D < 0.001949 with probability 0.999

We apply the chi-squared test to verify the distributions of the
geometric and Poisson generators (the Kolmogorov-Smirnov test
applies only to continuous distributions).

    Sampler              p-value
    rand_geometric(0.1)  0.883482
    rand_geometric(0.5)  0.527853
    rand_geometric(0.9)  0.651309
    rand_poisson(0.2)    0.700759
    rand_poisson(  1)    0.551421
    rand_poisson(  2)    0.619268
    rand_poisson( 20)    0.656092
    rand_poisson(200)    0.422257

Supposing the distributions are correct, the p-values should be
above zero with high probability:
    p-value > 0.05 with probability 0.95
    p-value > 0.01 with probability 0.99
    p-value > 0.001 with probability 0.999

Function Reference

init_randmt_auto

void init_randmt_auto()

seeds the global generator with the current time. It should be called once at the beginning of the program so that different pseudorandom values are produced on different runs.

It is enough to call this function once. Seeding multiple times does not improve statistical quality of the generator.

init_randmt

void init_randmt(unsigned long seed)

seeds the global random number generator with an unsigned 32-bit integer value.

A constant seed can be used to reproduce the same pseudorandom numbers.

int i;

init_randmt(42);
printf("Ten numbers:\n");
for (i = 0; i < 10; ++i) {
  printf("%f\n", rand_unif());
}

init_randmt(42);
printf("The same ten numbers:\n");
for (i = 0; i < 10; +++i) {
  printf("%f\n", rand_unif());
}

rand_unif

double rand_unif()

generates a random number uniformly on the open interval (0,1) with 53-bit resolution. The global generator is used to generate the value.

CDF \| \(\min(\max(x,0),1)\)
Mean	1/2
Variance	1/12

rand_uint32

unsigned long rand_uint32()

a random 32-bit unsigned integer uniformly between 0 and 0xFFFFFFFF. The global generator is used to generate the value.

rand_normal

double rand_normal()

generates a standard normal distributed random number with probability density

\[ f(x) = \tfrac{1}{\sqrt{2\pi}} \mathrm{e}^{-x^2 / 2}, \]

using the Box–Muller transform. The global generator is used to generate the value.

CDF \| \(\tfrac{1}{2}\bigl(1 + \mathrm{erf}(x/\sqrt{2})\bigr)\)
Mean	0
Variance	1

rand_exp

double rand_exp(double mu)

generates an exponentially-distributed number with mean mu with probability density

\[ f(x;\mu) = \begin{cases} \tfrac{1}{\mu} \mathrm{e}^{-x/mu} & x \ge 0, \\ 0 & x < 0, \end{cases} \]

by inversion. The global generator is used to generate the value.

CDF \| \(1 - \exp(-x/\mu)\)
Mean	\(\mu\)
Variance	\(\mu^2\)

rand_gamma

double rand_gamma(double a, double b)

generates a Gamma-distributed number with shape parameter a and scale parameter b with probability density

\[ f(x;a,b) = x^{a-1} \frac{\exp(-x/b)}{\Gamma(a)b^k}, \quad x \ge 0, \]

using the method of Marsaglia and Tsang, 2000. The global generator is used to generate the value.

CDF \| \(\gamma(a,x/b) / \Gamma(a)\)
Mean	\(ab\)
Variance	\(ab^2\)

where \(\gamma\) is the lower incomplete Gamma function.

rand_geometric

double rand_geometric(double p)

generates a geometrically-distributed number where p is the probability of success with

\[ P(X = k) = (1 - p)^{k-1} p, \quad k = 1, 2, \ldots \]

The global generator is used to generate the value.

CDF \| \(1 - (1 - p)^{\lfloor k \rfloor}\)
Mean	\(1/p\)
Variance	\((1-p)/p^2\)

rand_poisson

double rand_poisson(double mu)

generates a Poisson-distributed number with mean mu with

\[ P(X = k) = \frac{\mu^k \mathrm{e}^{-\mu}}{k!}, \quad k = 0, 1, \ldots \]

using a simple direct algorithm for mu < 10 and the “PTRS” algorithm of Hormann for larger mu. The global generator is used to generate the value.

CDF \| \(1 - \gamma(\lfloor k + 1\rfloor, \mu) / \lfloor k\rfloor !\)
Mean	μ
Variance	μ

new_randmt

randmt_t *new_randmt()

creates a new randmt_t, or returns NULL on memory allocation failure. The newly-created randmt_t should be initialized using init_randmt_auto_r or init_randmt_r. After use, call free_randmt to free memory associated with the randmt_t.

randmt_t *generator = NULL;
int i;

if (!(generator = new_randmt())) {
  exit(1);
}

init_randmt_auto_r(generator);

for (i = 0; i < 10; ++i) {
  printf("%f\n", rand_unif_r(generator));
}

free_randmt(generator);

free_randmt

void free_randmt(randmt_t *generator)

frees the memory associated with generator.

init_randmt_auto_r

void init_randmt_auto_r(randmt_t *generator)

initializes generator with the current time and memory address. The generator is seeded using the current time added to the memory address of the generator. The memory address is included so that different generators are initialized with different seeds. An array of generators can be created initialized as

randmt_t *generators[16];
int i;

for (i = 0; i < 16; ++i) {
  if(!(generators[i] = new_randmt())) {
    exit(1);
  }
}

for (i = 0; i < 16; i++) {
  init_randmt_auto_r(generators[i]);
}

init_randmt_r

void init_randmt_r(randmt_t *generator, unsigned long seed)

initializes generator with seed value seed.

rand_unif_r

double rand_unif_r(randmt_t *generator)

uses generator to generate a random number uniformly on the open interval (0,1) with 53-bit resolution.

rand_uint32_r

unsigned long rand_uint32_r(randmt_t *generator)

uses generator to generate a random 32-bit unsigned integer uniformly between 0 and 2³² − 1.

rand_normal_r

double rand_normal_r(randmt_t *generator)

uses generator to generate a standard normal distributed random number.

rand_exp_r

double rand_exp(randmt_t *generator, double mu)

uses generator to generate an exponentially-distributed number with mean mu.

rand_gamma_r

double rand_gamma(randmt_t *generator, double a, double b)

uses generator to generate a Gamma-distributed number with shape parameter a and scale parameter b.

rand_geometric_r

double rand_geometric(randmt_t *generator, double p)

uses generator to generate a geometrically-distributed number where p is the probability of success.

rand_poisson_r

double rand_poisson(randmt_t *generator, double mu)

uses generator to generate a Poisson-distributed number with mean mu.

This material is based upon work supported by the National Science Foundation under Award No. DMS-1004694. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.