randmt: MT19937 pseudorandom number generator
Pascal Getreuer
Overview
This C/C++ software contains a highquality 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_autofor (i = 0; i < 20; ++i) {
("%10.8f\n", rand_unif());
printf}
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()

32bit 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 scaleb
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
pseudorandom 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
withseed
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 chisquared 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 KolmogorovSmirnov 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 chisquared test to verify the distributions of the
geometric and Poisson generators (the KolmogorovSmirnov test
applies only to continuous distributions).
Sampler pvalue
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 pvalues should be
above zero with high probability:
pvalue > 0.05 with probability 0.95
pvalue > 0.01 with probability 0.99
pvalue > 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 32bit integer value.
A constant seed can be used to reproduce the same pseudorandom numbers.
int i;
(42);
init_randmt("Ten numbers:\n");
printffor (i = 0; i < 10; ++i) {
("%f\n", rand_unif());
printf}
(42);
init_randmt("The same ten numbers:\n");
printffor (i = 0; i < 10; +++i) {
("%f\n", rand_unif());
printf}
rand_unif
double rand_unif()
generates a random number uniformly on the open interval (0,1) with 53bit 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 32bit 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 exponentiallydistributed 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 Gammadistributed number with shape parameter
a
and scale parameter b
with probability
density
\[ f(x;a,b) = x^{a1} \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 geometricallydistributed number where p
is
the probability of success with
\[ P(X = k) = (1  p)^{k1} 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  \((1p)/p^2\) 
rand_poisson
double rand_poisson(double mu)
generates a Poissondistributed 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
*new_randmt() randmt_t
creates a new randmt_t
, or returns NULL
on
memory allocation failure. The newlycreated 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
.
*generator = NULL;
randmt_t int i;
if (!(generator = new_randmt())) {
(1);
exit}
(generator);
init_randmt_auto_r
for (i = 0; i < 10; ++i) {
("%f\n", rand_unif_r(generator));
printf}
(generator); free_randmt
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
*generators[16];
randmt_t int i;
for (i = 0; i < 16; ++i) {
if(!(generators[i] = new_randmt())) {
(1);
exit}
}
for (i = 0; i < 16; i++) {
(generators[i]);
init_randmt_auto_r}
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 53bit resolution.
rand_uint32_r
unsigned long rand_uint32_r(randmt_t *generator)
uses generator
to generate a random 32bit unsigned
integer uniformly between 0 and 2^{32} − 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 exponentiallydistributed
number with mean mu
.
rand_gamma_r
double rand_gamma(randmt_t *generator, double a, double b)
uses generator
to generate a Gammadistributed 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 geometricallydistributed
number where p
is the probability of success.
rand_poisson_r
double rand_poisson(randmt_t *generator, double mu)
uses generator
to generate a Poissondistributed number
with mean mu
.
This material is based upon work supported by the National Science Foundation under Award No. DMS1004694. 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.