Research‎ > ‎

randmt: MT19937 pseudorandom number generator


Download randmt

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.

Usage

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 usage:

    #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() the continuous uniform distribution on (0,1)
rand_uint32() 32-bit unsigned integers on [0,0xFFFFFFFF]
rand_normal() the standard normal distribution
rand_exp(mu) exponential distribution with mean mu
rand_gamma(a,b) Gamma distribution
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.” For these functions, the generator state is passed using a randmt_t object.

rand_unif_r(generator) the continuous uniform distribution on (0,1)
rand_uint32_r(generator) 32-bit unsigned integers on [0,0xFFFFFFFF]
rand_normal_r(generator) the standard normal distribution
rand_exp_r(generator,mu) exponential distribution with mean mu
rand_gamma_r(generator,a,b) Gamma distribution
rand_geometric_r(generator,p) geometric distribution
rand_poisson_r(generator,mu) Poisson distribution with mean 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 a specified seed value

History

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.

This function only needs to be called once. Seeding multiple times does not improve the 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

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

CDF
Mean 0
Variance   1

rand_exp

double rand_exp(double mu) generates an exponentially-distributed number with mean mu with probability density

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

CDF 1 − exp(−x/μ)
Mean μ
Variance   μ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

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

CDF γ(a,x/b)/Γ(a), where γ is the lower incomplete Gamma function
Mean ab
Variance   ab2

rand_geometric

double rand_geometric(double p) generates a geometrically-distributed number where p is the probability of success with

The global generator is used to generate the value.

CDF 1 − (1 − p)k
Mean 1/p
Variance   (1 − p)/p2

rand_poisson

double rand_poisson(double mu) generates a Poisson-distributed number with mean mu with

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

CDF 1 − γ(⌊k + 1⌋,μ)/⌊k⌋!
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 initiailized 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 0xFFFFFFFF.

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.