# 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) {
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

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

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

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

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

### 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

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

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

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

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

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

frees the memory associated with `generator`

.

### init_randmt_auto_r

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

initializes `generator`

with seed value `seed`

.

### rand_unif_r

uses `generator`

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

### rand_uint32_r

uses `generator`

to generate a random 32-bit unsigned integer uniformly between 0 and 2^{32} − 1.

### rand_normal_r

uses `generator`

to generate a standard normal distributed random number.

### rand_exp_r

uses `generator`

to generate an exponentially-distributed number with mean `mu`

.

### rand_gamma_r

uses `generator`

to generate a Gamma-distributed number with shape parameter `a`

and scale parameter `b`

.

### rand_geometric_r

uses `generator`

to generate a geometrically-distributed number where `p`

is the probability of success.

### rand_poisson_r

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.