Mersenne Twister

Source code for a C++ class implementing the Mersenne Twister MT19937 random number generator.

The original C implementation of the MT19937 by Takuji Nishimura and Makoto Matsumoto and further information about the algorithm can be found here: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html

Downloadable files (right click on the links and select "Save Link As"):

Content of STD_Random_MT19937.h:

//#############################################################################
//##  File Name:        STD_Random_MT19937.h                                 ##
//##  File Version:     1.00 (27 Jan 2008)                                   ##
//##  Author:           Silvestro Fantacci                                   ##
//##  Copyright:        BSD licence - see below                              ##
//#############################################################################
//
// I developed this C++ class by modifying the C souce code for the
// Mersenne Twister MT19937 by Makoto Matsumoto and Takuji Nishimura
// (ref. http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
// file mt19937ar.c - 2002 version), which carries the following copyright
// notice:
//
//   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
//   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.
//   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// I wrote this class with the intention of providing some cute (I hope!)
// souce code, as an alternative for people into C++ to browsing pseudocode
// descriptions of the Mersenne Twister algorithm, and for this reason I
// decided to sacrifice on efficiency with respect to the original code.
//
// I did not port the original genrand_int31, genrand_int32, genrand_real1,
// genrand_real2, genrand_real3 and genrand_res53 functions as I feel that in
// C++ they might be better placed in another more general class encapsulating
// this class and/or other random number generator classes.
//
// The initialisation of the random number generator with a list of seeds
// also differs from the original code in the following points:
// 1. if the seed list is empty, then the random number generator is
//    initialised with the default seed
// 2. if the seed list contains only one seed, then the random number
//    generator is initialised consistently with the single seed
//    initialisation
// 3. a spurious initial re-generation of the state vector is avoided
 
#if !defined (STD_RANDOM_MT19937_INCLUDED)
#define STD_RANDOM_MT19937_INCLUDED
 
#include <vector>
 
//=============================================================================
    class STD_Random_MT19937  
//=============================================================================
// Mersenne Twister MT19937 Random Number Generator
//
// Usage example:
//
//        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//        #include <iostream> // for cout
//        #include "STD_Random_MT19937.h"
//
//        void main ()
//        {
//            STD_Random_MT19937 Random;
//            for (unsigned int i = 1; i <= 10; i++)
//                std::cout << (unsigned int) Random.New_Number () << std::endl;            
//        }
//        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
{
public:
 
    // Constructs a random number generator and initialises it using a default
    // seed
                        STD_Random_MT19937 ();
 
    // Constructs a random number generator and initialises it using the
    // specified seed
                        STD_Random_MT19937 (unsigned __int32 Single_Seed);
 
    // Constructs a random number generator and initialises it using the
    // specified seed list
    // Equivalent to previous constructors for less than 2 seeds
    // Useful for 2 to 624 seeds, but more than 624 seeds can be used
                        STD_Random_MT19937
                            (const std::vector <unsigned __int32> & Seed_List);
 
    // Returns a new random number between 0 and 2^32 - 1
    // where 2^32 - 1 = 4294967295 = Hex FFFFFFFF
    unsigned __int32    New_Number ();
 
    // Re-initialises the random number generator using a default seed
    void                Seed ();
 
    // Re-initialises the random number generator using the specified seed
    void                Seed (unsigned __int32 Single_Seed);
 
    // Re-initialises the random number generator using the specified seed list
    // Equivalent to previous Seed methods for less than 2 seeds
    // Useful for 2 to 624 seeds, but more than 624 seeds can be used
    void                Seed (const std::vector <unsigned __int32> &
                                                            New_Seed_List);
 
    virtual                ~STD_Random_MT19937 ();
 
private:
 
    // Sets the whole MT19937 state vector based on the specified seed
    void                Generate_State (unsigned __int32 Seed);
 
    // Resets the whole MT19937 state vector based on its current state
    void                Regenerate_State ();
 
    // Resets the whole MT19937 state vector based on its current state
    // and the specified seed list
    void                Regenerate_State
                        (const std::vector <unsigned __int32> & Seed_List);
 
    // The MT19937 state vector
    unsigned __int32 *    State;
 
    // An index to the element in the MT19937 state vector which will be used
    // for calculating the number to be returned by New_Number ()
    unsigned int        Next;
};
 
#endif // !defined (STD_RANDOM_MT19937_INCLUDED)

Content of STD_Random_MT19937.cpp:

//#############################################################################
//##  File Name:        STD_Random_MT19937.cpp                               ##
//##  File Version:     1.00 (27 Jan 2008)                                   ##
//##  Author:           Silvestro Fantacci                                   ##
//##  Copyright:        BSD licence - see below                              ##
//#############################################################################
//
// I developed this C++ class by modifying the C souce code for the
// Mersenne Twister MT19937 by Makoto Matsumoto and Takuji Nishimura
// (ref. http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
// file mt19937ar.c - 2002 version), which carries the following copyright
// notice:
//
//   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
//   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 "STD_Random_MT19937.h"
 
// Size of the MT19937 state vector
#define SIZE            624
 
// SIZE_1 and SIZE_2 identify sections of the MT19937 state vector which are
// regenerated in a different way (see the Regenerate_State method).
// Note that SIZE_1 + SIZE_2 == SIZE 
#define SIZE_1            227
#define SIZE_2            397
 
//-----------------------------------------------------------------------------
    STD_Random_MT19937::STD_Random_MT19937 ()
//-----------------------------------------------------------------------------
{
    State = new unsigned __int32 [SIZE];
    Seed ();
}
 
//-----------------------------------------------------------------------------
    STD_Random_MT19937::STD_Random_MT19937 (unsigned __int32 Single_Seed)
//-----------------------------------------------------------------------------
{
    State = new unsigned __int32 [SIZE];
    Seed (Single_Seed);
}
 
//-----------------------------------------------------------------------------
    STD_Random_MT19937::STD_Random_MT19937
                        (const std::vector <unsigned __int32> & Seed_List)
//-----------------------------------------------------------------------------
{
    State = new unsigned __int32 [SIZE];
    Seed (Seed_List);
}
 
//-----------------------------------------------------------------------------
    unsigned __int32 STD_Random_MT19937::New_Number ()
//-----------------------------------------------------------------------------
{
    unsigned __int32 Result;
 
    Result = State [Next];
 
    if (++Next == SIZE)
        Regenerate_State ();
 
    // Tempering
    Result ^=  Result >> 11;
    Result ^= (Result << 7)  & 0x9D2C5680;
    Result ^= (Result << 15) & 0xEFC60000;
    Result ^= (Result >> 18);
 
    return Result;
}
 
//-----------------------------------------------------------------------------
    void STD_Random_MT19937::Seed ()
//-----------------------------------------------------------------------------
{
    Seed (5489);
}
 
//-----------------------------------------------------------------------------
    void STD_Random_MT19937::Seed (unsigned __int32 Single_Seed)
//-----------------------------------------------------------------------------
{
    Generate_State (Single_Seed);
    Regenerate_State ();
}
 
//-----------------------------------------------------------------------------
    void STD_Random_MT19937::Seed
                        (const std::vector <unsigned __int32> & Seed_List)
//-----------------------------------------------------------------------------
{
    if (Seed_List.size () == 0)
 
        Seed ();
 
    else if (Seed_List.size () == 1)
 
        Seed (Seed_List [0]);
    else
    {
        Generate_State (19650218);
        Regenerate_State (Seed_List);
    }
}
 
//-----------------------------------------------------------------------------
    void STD_Random_MT19937::Generate_State (unsigned __int32 Seed)
//-----------------------------------------------------------------------------
{
    State [0] = Seed;
    for (unsigned int i = 1; i < SIZE; i++)
    {
        State [i] = State [i - 1] ^ (State [i - 1] >> 30);
        State [i] = State [i] * 1812433253 + i;
    }
 
    Next = 0;
}
 
//-----------------------------------------------------------------------------
    void STD_Random_MT19937::Regenerate_State ()
//-----------------------------------------------------------------------------
{
 
    #define MSB_ONLY    0x80000000    // Mask for the most significant bit
    #define NO_MSB        0x7FFFFFFF    // Mask for the 31 least significant bits    
    #define CRNFTM        0x9908B0DF    // Coefficient of the rational normal form
                                    // twist matrix 
 
    unsigned __int32    s;
    unsigned int        i;
 
    // Regenerate the first SIZE_1 elements
    for (i = 0; i < SIZE_1; i++)
    {
        s = (State [i] & MSB_ONLY) | (State [i + 1] & NO_MSB);
        State [i] = State [i + SIZE_2] ^ (s >> 1);
        if (s % 2 == 1)
            State [i] ^= CRNFTM;
    }
 
    // Regenerate the next SIZE_2 - 1 elements
    for (i = SIZE_1; i < SIZE - 1; i++)
    {
        s = (State [i] & MSB_ONLY) | (State [i + 1] & NO_MSB);
        State [i] = State [i - SIZE_1] ^ (s >> 1);
        if (s % 2 == 1)
            State [i] ^= CRNFTM;
    }
 
    // Regenerate the last element
    s = (State [SIZE - 1] & MSB_ONLY) | (State [0] & NO_MSB);
    State [SIZE - 1] = State [SIZE_2 - 1] ^ (s >> 1);
    if (s % 2 == 1)
        State [SIZE - 1] ^= CRNFTM;
 
    Next = 0;
}
 
//-----------------------------------------------------------------------------
    void STD_Random_MT19937::Regenerate_State
                            (const std::vector <unsigned __int32> & Seed_List)
//-----------------------------------------------------------------------------
{
    unsigned int i, j, k;
 
    i = 1;
    j = 0;
 
    const unsigned int Max = (SIZE >= Seed_List.size () ?   SIZE
                                                          : Seed_List.size ());
 
    for (k = 1; k <= Max; k++)
    {
        State [i] =   (State [i] ^ (  (State [i - 1] ^ (State [i - 1] >> 30))
                                    * 1664525))
                    + Seed_List [j]
                    + j;
 
        if (++i == SIZE)
        {
            State [0] = State [SIZE - 1];
            i = 1;
        }
 
        if (++j == Seed_List.size ())
            j = 0;
    }
 
    for (k = 1; k < SIZE; k++)
    {
        State [i] =   (State [i] ^ (  (State [i - 1] ^ (State [i - 1] >> 30))
                                    * 1566083941))
                    - i;
 
        if (++i == SIZE)
        {
            State [0] = State [SIZE - 1];
            i = 1;
        }
    }
 
    State [0] = 0x80000000;
 
    Next = 0;
}
 
//-----------------------------------------------------------------------------
    STD_Random_MT19937::~STD_Random_MT19937 ()
//-----------------------------------------------------------------------------
{
    delete [] State;
}
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License