https://doi.org/10.5201/ipol.2012.g-tvi
Raw File
Tip revision: 2cc1be6ae636163bfaaf300aeb376dbea4f887f9 authored by Software Heritage on 01 July 2012, 00:00:00 UTC
ipol: Deposit 1195 in collection ipol
Tip revision: 2cc1be6
randmt.c
/**
 * \warning Do NOT use for cryptographic purposes.  Read Internet RFC4086,
 * <http://tools.ietf.org/html/rfc4086>.
 */

/**
 * \file randmt.c
 * \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>
 * 
 * \par History
 * - This code uses the Mersenne Twister MT19937 code by Makoto Matsumoto and
 *   Takuji Nishimura (2002/1/26), the original is available at
 *   <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html>
 * - Seiji Nishimura (2008) wrote a reentrant modification of MT19937, where
 *   a generator struct is passed to the pseudorandom sampler functions,
 *   <http://www.k3.dion.ne.jp/~sn1976/site/Software.html>
 * - Nicolas Limare wrote a version with Doxygen documentation, simpler
 *   interface, and automatic initialization using the time, 
 *   <http://dev.ipol.im/git/?p=nil/mt.git;a=summary>
 * - Pascal Getreuer added samplers for normal, exponential, Gamma, 
 *   geometric, and Poisson distributions and test program with Kolmogorov-
 *   Smirnov and chi-squared testing.
 * 
 * 
 * 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.
 */


#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <math.h>
#include <time.h>
#include "randmt.h"

/* Get the high-precision timer specified by POSIX.1-2001 if available */
#if defined(USE_GETTIMEOFDAY) || defined(unix) || defined(__unix__) || defined(__unix)
#include <unistd.h>
#if defined(USE_GETTIMEOFDAY) || (defined(_POSIX_VERSION) && (_POSIX_VERSION >= 200112L))
#include <sys/time.h>
#define USE_GETTIMEOFDAY
#endif
#endif

#ifndef M_LOGSQRT2PI
/** \brief The constant \f$ \log \sqrt{2\pi} \f$ */
#define M_LOGSQRT2PI    0.9189385332046727417803297
#endif

/* Period parameters */
#define MT_N           624             /**< length of the state vector    */
#define MT_M           397
#define MT_MATRIX_A    0x9908B0DFUL    /**< constant vector a             */
#define MT_UPPER_MASK  0x80000000UL    /**< most  significant w-r bits    */
#define MT_LOWER_MASK  0x7FFFFFFFUL    /**< least significant   r bits    */


/** \brief An MT19937 pseudo-random number generator */
typedef struct randmtstruct_t
{    
    unsigned long mt[MT_N];    /**< the array for the state vector       */    
    int mti;                   /**< current position in the state vector */
} randmttype_t;

/** \brief Global randmt_t, used with the global versions of the functions */
randmt_t __randmt_global_generator = {{0}, MT_N + 1};


/* Create a new randmt_t */
randmt_t *new_randmt(void)
{
    randmt_t *generator;

    if((generator = (randmt_t *)calloc(1, sizeof(randmt_t))))
    {
        /* mti == MT_N+1 means mt[MT_N] is not initialized */
        generator->mti = MT_N + 1;
    }

    return generator;
}


/* Free a randmt_t */
void free_randmt(randmt_t *generator)
{
    if(generator)
        free(generator);
    return;
}


/* Initialize randmt_t with a seed */
void init_randmt_r(randmt_t *generator, unsigned long seed)
{
    int mti;
    unsigned long *mt;

    if(!generator)
        return;
    
#if ULONG_MAX > 0xFFFFFFFFUL
    seed &= 0xFFFFFFFFUL;           /* for WORDSIZE > 32bit */
#endif
    (mt = generator->mt)[0] = seed;

    for(mti = 1; mti < MT_N; mti++) 
    {
        mt[mti] = (1812433253UL * (mt[mti - 1] ^ (mt[mti - 1] >> 30)) + mti);
        /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
         * In the previous versions, MSBs of the seed affect
         * only MSBs of the array mt[].
         * 2002/01/09 modified by Makoto Matsumoto
         */
#if ULONG_MAX > 0xFFFFFFFFUL
        mt[mti] &= 0xFFFFFFFFUL;    /* for WORDSIZE > 32bit machines */
#endif
    }

    generator->mti = MT_N;
    return;
}


/* Initialize generator with the current time and memory address */
void init_randmt_auto_r(randmt_t *generator)
{
    unsigned long int seed;

    /* The basic (weak) initialization uses the current time, in seconds */
    seed = (unsigned long int) time(NULL);
    /* Add to the generator's memory address so that different generators 
       use different seeds. */
    seed += ((unsigned long int) generator);

#if defined(USE_GETTIMEOFDAY)
    /* gettimeofday() provides microsecond resolution */
    {
        struct timeval tp;
        (void) gettimeofday(&tp, NULL);

        seed *= 1000000;
        seed += tp.tv_usec;
    }
#endif

    init_randmt(seed);
    return;
}


/* Generate a random 32-bit unsigned integer (reentrant) */
unsigned long rand_uint32_r(randmt_t *generator)
{
    /* mag01[x] = x * MT_MATRIX_A  for x=0,1 */
    const unsigned long mag01[2] = { 0x0UL, MT_MATRIX_A };
    unsigned long y, *mt;
    int mti;
    
    if(!generator)
        return 0;

    mt = generator->mt;
    mti = generator->mti;

    if (mti >= MT_N)   /* generate MT_N words at one time */
    {    
        int kk;
        
        /* If init_randmt_r() has not been called, use a default seed. */
        if (mti == MT_N + 1)   
            init_randmt_r(generator, 5489UL);
        
        for(kk = 0; kk < MT_N - MT_M; kk++) 
        {
            y = (mt[kk] & MT_UPPER_MASK) | (mt[kk + 1] & MT_LOWER_MASK);
            mt[kk] = mt[kk + MT_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
        }        
        for(; kk < MT_N - 1; kk++) 
        {
            y = (mt[kk] & MT_UPPER_MASK) | (mt[kk + 1] & MT_LOWER_MASK);
            mt[kk] = mt[kk + (MT_M - MT_N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
        }
        
        y = (mt[MT_N - 1] & MT_UPPER_MASK) | (mt[0] & MT_LOWER_MASK);
        mt[MT_N - 1] = mt[MT_M - 1] ^ (y >> 1) ^ mag01[y & 0x1UL];
        mti = 0;
    }

    y = mt[mti++];
    
    /* Tempering */
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9D2C5680UL;
    y ^= (y << 15) & 0xEFC60000UL;    
    y ^= (y >> 18);

    generator->mti = mti;

    return y;
}


/* Generate a Gamma-distributed number (reentrant) */
double rand_gamma_r(randmt_t *generator, double a, double b)
{
    const double d = a - 1.0/3.0;
    const double c = 1.0/sqrt(9*d);
    double x, v, u;
    
    if(a >= 1)
    {
        do
        {
            do
            {
                x = rand_normal_r(generator);
                v = 1 + c*x;
            }while(v <= 0);
            
            v = v*v*v;
            u = rand_unif_r(generator);
            x = x*x;
        }while(u >= 1 - 0.0331*x*x 
            && log(u) >= x/2 + d*(1 - v + log(v)));
                
        return b*d*v;
    }
    else if(a > 0)
        return rand_gamma_r(generator, 1 + a, b) 
            * pow(rand_unif_r(generator), 1/a);
    else
        return 0;
}


/* Compute the natural logarithm of the factorial (used by rand_poisson_r) */
static double logfactorial(double n)
{
    /* Look-up table for log(n!) for n = 2, ..., 11 */
    static const double Table[10] = {
        0.69314718055994530942, 1.7917594692280550008,
        3.1780538303479456197, 4.7874917427820459943,
        6.5792512120101009951, 8.5251613610654143002,
        10.604602902745250228, 12.801827480081469611,
        15.104412573075515295, 17.502307845873885839};
    /* Asymptotic approximation coefficients */
    static const double C[7] = { -1.910444077728e-03, 
        8.4171387781295e-04, -5.952379913043012e-04, 
        7.93650793500350248e-04, -2.777777777777681622553e-03, 
        8.333333333333333331554247e-02, 5.7083835261e-03};
    double x, xsq, res, corr;
    int i;
    
    if(n > 11)  /* Asymptotic approximation, based on the public domain */
    {           /* NETLIB (Fortran) code by W. J. Cody and L. Stoltz    */
        x = 1 + floor(n + 0.5);
        
        if(x <= 2.25e76)
        {
            res = C[6];
            xsq = x*x;
            
            for(i = 0; i < 6; i++)
                res = res/xsq + C[i];
        }
        else
            res = 0;
        
        corr = log(x);
        return res/x + M_LOGSQRT2PI - corr/2 + x*(corr - 1);
    }
    else        /* Use look-up table */
    {
        i = (int)floor(n + 0.5);
        return (i >= 2) ? Table[i - 2] : 0; 
    }
}


/* Generate a Poisson-distributed number (reentrant) */
double rand_poisson_r(randmt_t *generator, double mu)
{
    double k, b, a, U, V, us, vr, invalpha, logmu;
    
    if(mu < 10)
    {   /* Use simple direct algorthm for small mu */
        vr = exp(-mu);
        k = -1;
        V = 1;
        
        do
        {
            k += 1;
            V *= rand_unif_r(generator);
        }while(V > vr);
    }
    else
    {   /* Use the "PTRS" algorithm of Hormann */
        b = 0.931 + 2.53*sqrt(mu);
        a = -0.059 + 0.02483*b;
        vr = 0.9277 - 3.6224/(b - 2);
        invalpha = 1.1239 + 1.1328/(b - 3.4);
        logmu = log(mu);
        
        do
        {
            U = rand_unif_r(generator) - 0.5;
            V = rand_unif_r(generator);
            us = 0.5 - fabs(U);
            k = floor((2*a/us + b)*U + mu + 0.43);
        }while((us < 0.07 || V > vr)
            && (k < 0 || (us < 0.013 && V > us)
                || log(V*invalpha/(a/(us*us) + b)) 
                        > -mu + k*logmu - logfactorial(k)));
    }
    
    return k;
}
back to top