https://doi.org/10.5201/ipol.2012.g-tvi
Tip revision: 2cc1be6ae636163bfaaf300aeb376dbea4f887f9 authored by Software Heritage on 01 July 2012, 00:00:00 UTC
ipol: Deposit 1195 in collection ipol
ipol: Deposit 1195 in collection ipol
Tip revision: 2cc1be6
randmt.h
/**
* \warning Do NOT use for cryptographic purposes. Read Internet RFC4086,
* <http://tools.ietf.org/html/rfc4086>.
*/
/**
* \file randmt.h
* \brief Mersenne Twister MT19937 pseudorandom number generator
* \author Makoto Matsumoto (1997-2002)
* \author Takuji Nishimura (1997-2002)
* \author Seiji Nishimura (2008) <seiji1976@gmail.com>
* \author Nicolas Limare <nicolas.limare@cmla.ens-cachan.fr>
* \author Pascal Getreuer <getreuer@gmail.com>
*
* This software provides a high-quality pseudorandom number generator, the
* Mersenne Twister MT19937.
*
* \par Basic Usage
* To initialize the generator, call init_randmt_auto() at the beginning of the
* program. Example usage:
@code
#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;
}
@endcode
* Pseudorandom numbers can be generated for several different distributions:
* - #rand_unif, the continuous uniform distribution on (0,1)
* - #rand_uint32, 32-bit unsigned integers uniformly on [0,0xFFFFFFFF]
* - #rand_normal, the standard normal distribution
* - #rand_exp, exponential distribution
* - #rand_gamma, Gamma distribution
* - #rand_geometric, geometric distribution
* - #rand_poisson, Poisson distribution
*
* \par 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, the continuous uniform distribution on (0,1)
* - #rand_uint32_r, 32-bit unsigned integers uniformly on [0,0xFFFFFFFF]
* - #rand_normal_r, the standard normal distribution
* - #rand_exp_r, exponential distribution
* - #rand_gamma_r, Gamma distribution
* - #rand_geometric_r, geometric distribution
* - #rand_poisson_r, Poisson distribution
*/
#ifndef _RANDMT_H_
#define _RANDMT_H_
#ifndef M_PI
/** \brief The constant \f$ \pi \f$ */
#define M_PI 3.14159265358979323846264338327950288419176939937510
#endif
/**
* \brief An MT19937 pseudo-random number generator
*
* This 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
* - #init_randmt_auto_r, initialize randmt_t with time and address
* - #init_randmt_r, initialize randmt_t with a specified seed value
* - #free_randmt, free memory associated with a randmt_t
*
* The randmt_t is encapsulated by forward declaration, and defined in
* randmt.c.
*
*
* Copyright (C) 1997-2002, Makoto Matsumoto and Takuji Nishimura \n
* Copyright (C) 2008, Seiji Nishimura <seiji1976@gmail.com> \n
* Copyright (C) 2010-2011 Nicolas Limare <nicolas.limare@cmla.ens-cachan.fr> \n
* Copyright (C) 2011, Pascal Getreuer <getreuer@gmail.com> \n
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
*
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* 3. The names of its contributors may not be used to endorse or promote
* products derived from this software without specific prior written
* permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
* A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
* OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
* EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
* PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
* LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
* NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
typedef struct randmtstruct_t randmt_t;
/** \brief Global randmt_t, used with the global versions of the functions */
extern randmt_t __randmt_global_generator;
/* Samplers using a specified generator */
/**
* \defgroup randmtlocal randmt: Samplers using a specified generator
*
* These functions use a specified #randmt_t to generate pseudorandom numbers.
* \{
*/
/**
* \brief Create a new randmt_t
* \return a randmt_t pointer, or NULL on failure
* \see init_randmt_auto_r()
* \see init_randmt_r()
* \see free_randmt()
*
* A newly-created randmt_t should be seeded with init_randmt_auto_r() or
* init_randmt_r(). After use, call free_randmt() to free memory associated
* with the randmt_t.
@code
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);
@endcode
*/
randmt_t *new_randmt();
/**
* \brief Free a randmt_t
* \param generator the randmt_t to free
* \see new_randmt()
*
* Free memory associated to a randmt_t that was created using new_randmt().
*/
void free_randmt(randmt_t *generator);
/**
* \brief Initialize randmt_t with a seed
* \param generator the randmt_t
* \param seed the seed value
* \see init_randmt()
* \see init_randmt_auto_r()
*/
void init_randmt_r(randmt_t *generator, unsigned long seed);
/**
* \brief Initialize generator with the current time and memory address
* \param generator the randmt_t
* \see init_randmt_auto()
* \see init_randmt_r()
*
* 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 initialized as
@code
randmt_t *generators[16];
int i;
for(i = 0; i < 16; i++)
init_randmt_auto_r(generators[i]);
@endcode
*/
void init_randmt_auto_r(randmt_t *generator);
/**
* \brief Generate a random 32-bit unsigned integer (reentrant)
* \param generator the randmt_t
* \return Pseudo-random integer value uniformly between 0 and 0xFFFFFFFF
* \see rand_uint32()
*
* Generates a pseudorandom 32-bit integer value uniformly between 0 and
* 0xFFFFFFFF using the specified generator.
*/
unsigned long rand_uint32_r(randmt_t *generator);
/**
* \brief Generate a Gamma-distributed number (reentrant)
* \param generator the randmt_t
* \param a shape parameter (positive value)
* \param b scale parameter (positive value)
* \return Pseudo-random Gamma-distibuted value
* \see rand_gamma()
*
* Generates a Gamma-distributed random value using the specified generator.
*/
double rand_gamma_r(randmt_t *generator, double a, double b);
/**
* \brief Generate a Poisson-distributed number (reentrant)
* \param generator the randmt_t
* \param mu the mean parameter of the distribution (positive value)
* \return Pseudo-random Poisson-distibuted value
* \see rand_poisson()
*
* A pseudo-random normally distributed number using the specified generator.
*/
double rand_poisson_r(randmt_t *generator, double mu);
/**
* \brief Generate a uniform random number on (0,1) (reentrant)
* \param generator the randmt_t
* \return Pseudo-random double value in the open interval (0,1)
* \see rand_unif()
*
* This routine generates a random number uniformly on the open interval
* (0,1) with 53-bit resolution using the specified generator. The formula
* used to generate the number is
* \f[ \mathtt{rand\_unif} = (a 2^{26} + b + 0.5 - \epsilon) / 2^{53}, \f]
* where a = integer of 27 random bits, b = integer of 26 random bits.
*/
#define rand_unif_r(generator) ( \
(((rand_uint32_r(generator) >> 5)*67108864.0 \
+ (rand_uint32_r(generator) >> 6)) \
+ 0.4999999999999998) / 9007199254740992.0)
/**
* \brief Generate a standard normal distributed random number (reentrant)
* \param generator the randmt_t
* \return Pseudo-random double value with standard normal distribution
* \see rand_normal()
*
* A pseudo-random normally distributed number using the specified generator.
*/
#define rand_normal_r(generator) ( \
sqrt(-2.0*log(rand_unif_r(generator))) \
* cos(2.0*M_PI*rand_unif_r(generator)))
/**
* \brief Generate an exponentially-distributed number (reentrant)
* \param generator the randmt_t
* \param mu mean parameter of the distribution (positive value)
* \return Pseudo-random exponentially-distibuted value with mean mu
* \see rand_exp()
*
* Generates a pseudo-random exponentially distributed value using the
* specified generator.
*/
#define rand_exp_r(generator, mu) ( \
-(mu)*log(rand_unif_r(generator)))
/**
* \brief Generate a geometrically-distributed number (reentrant)
* \param generator the randmt_t
* \param p probability of success
* \return Pseudo-random geometrically-distibuted value
* \see rand_geometric()
*
* Generates a pseudo-random geometrically-distributed value using the
* specified generator.
*/
#define rand_geometric_r(generator, p) ( \
floor(log(rand_unif_r(generator))/log(1 - p)) + 1)
/** \} */
/* Samplers using the global generator */
/**
* \defgroup randmtglobal randmt: Samplers using the global generator
*
* These functions use the global pseudorandom number generator.
* \{
*/
/**
* \brief Initialize the global generator with a seed
* \param seed the seed value
* \see init_randmt_auto()
* \see init_randmt_r()
*
* This routine 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.
@code
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());
@endcode
*/
#define init_randmt(seed) (init_randmt_r(&__randmt_global_generator, seed))
/**
* \brief Initialize the global generator with the current time
* \see init_randmt()
* \see init_randmt_auto_r()
*
* This function 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.
*/
#define init_randmt_auto() (init_randmt_auto_r(&__randmt_global_generator))
/**
* \brief Generate a random 32-bit unsigned integer (nonreentrant)
* \return Pseudo-random integer value uniformly between 0 and 0xFFFFFFFF
* \see rand_uint32_r()
* \see rand_unif()
*
* The global generator is used to generate the value.
*/
#define rand_uint32() (rand_uint32_r(&__randmt_global_generator))
/**
* \brief Generate a uniform random number on (0,1) (nonreentrant)
* \return Pseudo-random double value in the open interval (0,1)
* \see rand_unif_r()
* \see rand_uint32()
*
* This routine generates a random number uniformly on the open interval
* (0,1) with 53-bit resolution. The formula used to generate the number is
* \f[ \mathtt{rand\_unif} = (a 2^{26} + b + 0.5 - \epsilon) / 2^{53}, \f]
* where a = integer of 27 random bits, b = integer of 26 random bits.
* The global generator is used to generate the value.
*
* Distribution properties:
* - CDF: \f$ \min(\max(x,0),1) \f$
* - Mean: \f$ \frac{1}{2} \f$
* - Variance: \f$ \frac{1}{12} \f$
*/
#define rand_unif() (rand_unif_r(&__randmt_global_generator))
/**
* \brief Generate a standard normal distributed random number (nonreentrant)
* \return Pseudo-random double value with standard normal distribution
* \see rand_normal_r()
*
* A pseudo-random normally distributed number with density
* \f[ f(x) = \frac{1}{\sqrt{2\pi}} \mathrm{e}^{-x^2/2}, \f]
* generated by the Box-Muller transform. The global generator is used
* to generate the value.
*
* Distribution properties:
* - CDF: \f$ \frac{1}{2}\bigl( 1 + \mathrm{erf}(x/\sqrt{2})\bigr) \f$
* - Mean: 0
* - Variance: 1
*/
#define rand_normal() (rand_normal_r(&__randmt_global_generator))
/**
* \brief Generate an exponentially-distributed number (nonreentrant)
* \param mu mean parameter of the distribution (positive value)
* \return Pseudo-random exponentially-distibuted value with mean mu
* \see rand_exp_r()
*
* Generates a pseudo-random value distributed with probability density
* \f[ f(x;\mu) = \begin{cases}
* \frac{1}{\mu} \mathrm{e}^{-x/\mu}, & x \ge 0, \\
* 0, & x < 0,
* \end{cases} \f]
* generated by inversion. The global generator is used to generate
* the value.
*
* Distribution properties:
* - CDF: \f$ 1 - \exp(-x/\mu) \f$
* - Mean: \f$ \mu \f$
* - Variance: \f$ \mu^2 \f$
*/
#define rand_exp(mu) (rand_exp_r(&__randmt_global_generator, mu))
/**
* \brief Generate a Gamma-distributed number (nonreentrant)
* \param a shape parameter (positive value)
* \param b scale parameter (positive value)
* \return Pseudo-random Gamma-distibuted value
* \see rand_gamma_r()
*
* Generates a Gamma-distributed random value with density
* \f[ f(x;a,b) = x^{a-1} \frac{\exp(-x/b)}{\Gamma(a)b^k}, \quad x \ge 0, \f]
* where \f$ \Gamma(a) \f$ is the Gamma function, using the method of
* Marsaglia and Tsang, 2000. The global generator is used to generate
* the value.
*
* Distribution properties:
* - CDF: \f$ \gamma(a,x/b)/\Gamma(a) \f$, where \f$ \gamma \f$ is the
* lower incomplete gamma function
* - Mean: \f$ ab \f$
* - Variance: \f$ ab^2 \f$
*/
#define rand_gamma(a, b) (rand_gamma_r(&__randmt_global_generator, a, b))
/**
* \brief Generate a Poisson-distributed number (nonreentrant)
* \param mu the mean parameter of the distribution (positive value)
* \return Pseudo-random Poisson-distibuted value
* \see rand_poisson_r()
*
* A pseudo-random normally distributed number with probability mass function
* \f[ \mathrm{P}(X = k) = \frac{\mu^k \mathrm{e}^{-\mu}}{k!},
* \quad k = 0, 1,\ldots, \f]
* 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.
*
* Reference:
* Wolfgang Hormann, "The transformed rejection method for generating
* Poisson random variables," Insurance: Mathematics and Economics 12,
* pp. 39-45, 1993.
*
* Distribution properties:
* - CDF: \f$ 1 - \gamma(\lfloor k+1 \rfloor, \mu)/\lfloor k \rfloor ! \f$,
* where \f$ \gamma \f$ is the lower incomplete Gamma function
* - Mean: \f$ \mu \f$
* - Variance: \f$ \mu \f$
*/
#define rand_poisson(mu) (rand_poisson_r(&__randmt_global_generator, mu))
/**
* \brief Generate a geometrically-distributed number (nonreentrant)
* \param p probability of success
* \return Pseudo-random geometrically-distibuted value
* \see rand_geometric_r()
*
* Generates a pseudo-random value distributed with probability mass fuction
* \f[ \mathrm{P}(X = k) = (1 - p)^{k-1} p, \quad k = 1, 2, \ldots, \f]
* where \f$ 0 < p \le 1 \f$, generated by inversion. The global generator
* is used to generate the value.
*
* Distribution properties:
* - CDF: \f$ 1 - (1 - p)^{\lfloor k \rfloor} \f$
* - Mean: \f$ 1/p \f$
* - Variance: \f$ (1 - p)/p^2 \f$
*/
#define rand_geometric(p) (rand_geometric_r(&__randmt_global_generator, p))
/** \} */
#endif /* _RANDMT_H_ */