# randmt: MT19937 pseudorandom number generator

## Pascal Getreuer

## 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) {
"%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()`

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;
42);
init_randmt("Ten numbers:\n");
printf(for (i = 0; i < 10; ++i) {
"%f\n", rand_unif());
printf(
}
42);
init_randmt("The same ten numbers:\n");
printf(for (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 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())) {
1);
exit(
}
init_randmt_auto_r(generator);
for (i = 0; i < 10; ++i) {
"%f\n", rand_unif_r(generator));
printf(
}
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

```
16];
randmt_t *generators[int i;
for (i = 0; i < 16; ++i) {
if(!(generators[i] = new_randmt())) {
1);
exit(
}
}
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^{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 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.