Mersenne Twister

Mersenne Twister

Documentation

#include <cryptopp/mersenne.h>

A Mersenne Twister is a pseudo random number generator designed by Makoto Matsumoto and Takuji Nishimura suitable for Monte-Carlo simulations. The generator is not cryptographically secure.

The Twister is provided for interoperability and to demonstrate the creation of a generator.

A related wiki page is RandomNumberGenerator, which discusses the RandomNumberGenerator interface and generating random numbers.

Requirements

The generator is somewhat tricky to implement because it is word oriented, and not byte oriented. Because the generator is word oriented, there are two implications for an implementation. First, the result of GenerateWord32 must be consistent with the result of calling GenerateBlock with 1, 2, 3 and 4 byte arrays on both big endian and little endian machines. For example, if GenerateWord32 returns 0xD091BB5C, then GenerateBlock must return 0xD0 0x91 0xBB 0x5C for 1, 2, 3 and 4 byte arrays. Second, Discard rounds up to a multiple of a word size, and then discards the required number of words (and not bytes).

The Mersenne Twister also needs to parameterize some magic constants. The magic constants are K used in the twist; M used as the period parameter; N used as the state size; F used as a multiplier; and S used as a default seed.

Implementation

The class declaration, with its magic constants is as follows.

template <unsigned int K, unsigned int M, unsigned int N, unsigned int F, unsigned long S>
class MersenneTwister : public RandomNumberGenerator
{
  public:
    MersenneTwister(unsigned long seed = S);
    virtual ~MersenneTwister();

    virtual void GenerateBlock(byte *output, size_t size);
    virtual word32 GenerateWord32(word32 min=0, word32 max=0xffffffffL);
    virtual void DiscardBytes(size_t n);

  protected:
    word32 NextMersenneWord();
    void Twist();

  private:
    FixedSizeSecBlock<word32, N> m_state;
    unsigned int m_seed, m_idx;
};

The constructor sets the initial state. In Matsumoto and Nishimura's paper, this is the code that executes when index == N+1.

MersenneTwister(unsigned long seed = S) : m_seed(seed), m_idx(N)
{
    m_state[0] = seed;
    for (unsigned int i = 1; i < N; i++)
        m_state[i] = word32(F * (m_state[i-1] ^ (m_state[i-1] >> 30)) + i);
}

NextMersenneWord

The implementation of the Mersenne Twister is centered around NextMersenneWord. Both GenerateWord32 and GenerateBlock call NextMersenneWord. NextMersenneWord returns the next word from the state table. If the index exceeds the state size, then Twist is performed and the index is set to 0.

word32 NextMersenneWord()
{
    if (m_idx >= N) { Twist(); }
		
    word32 temp = m_state[m_idx++];
    temp ^= (temp >> 11);
    temp ^= (temp <<  7) & 0x9D2C5680; // 2636928640
    temp ^= (temp << 15) & 0xEFC60000; // 4022730752
		
    return temp ^ (temp >> 18);
}

void Twist()
{			
    static const unsigned long magic[2]={0x0UL, K};
    word32 kk, temp;

    for (kk=0;kk<N-M;kk++)
    {
        temp = (m_state[kk] & 0x80000000)|(m_state[kk+1] & 0x7FFFFFFF);
        m_state[kk] = m_state[kk+M] ^ (temp >> 1) ^ magic[temp & 0x1UL];
    }
		
    for (;kk<N-1;kk++)
    {
        temp = (m_state[kk] & 0x80000000)|(m_state[kk+1] & 0x7FFFFFFF);
        m_state[kk] = m_state[kk+(M-N)] ^ (temp >> 1) ^ magic[temp & 0x1UL];
    }
		
    temp = (m_state[N-1] & 0x80000000)|(m_state[0] & 0x7FFFFFFF);
    m_state[N-1] = m_state[M-1] ^ (temp >> 1) ^ magic[temp & 0x1UL];
		
    // Reset index
    m_idx = 0;
}

GenerateWord32

GenerateWord32 simply returns the next word from the state array, or the next word that meets the requirements of min and max:

word32 GenerateWord32(word32 min=0, word32 max=0xffffffffL)
{
    const word32 range = max-min;
    if (range == 0xffffffffL)
        return NextMersenneWord();
			
    const int maxBits = BitPrecision(range);
    word32 value;

    do {
        value = Crop(NextMersenneWord(), maxBits);
    } while (value > range);

    return value+min;
}

GenerateBlock

GenerateBlock is more interesting because of the big endian/little endian/array consistency requirements. In the code below, ByteReverse and CRYPTOPP_GET_BYTE_AS_BYTE are provided by the library in misc.h.

void GenerateBlock(byte *output, size_t size)
{
    // Handle word32 size blocks
    word32 temp;
    for (size_t i=0; i < size/4; i++, output += 4)
    {
#if defined(CRYPTOPP_ALLOW_UNALIGNED_DATA_ACCESS) && defined(IS_LITTLE_ENDIAN)
        *((word32*)output) = ByteReverse(NextMersenneWord());
#elif defined(CRYPTOPP_ALLOW_UNALIGNED_DATA_ACCESS)
        *((word32*)output) = NextMersenneWord();
#else
        temp = NextMersenneWord();
        output[3] = CRYPTOPP_GET_BYTE_AS_BYTE(temp, 0);
        output[2] = CRYPTOPP_GET_BYTE_AS_BYTE(temp, 1);
        output[1] = CRYPTOPP_GET_BYTE_AS_BYTE(temp, 2);
        output[0] = CRYPTOPP_GET_BYTE_AS_BYTE(temp, 3);
#endif
    }
		
    // No tail bytes
    if (size%4 == 0)
        return;
		
    // Handle tail bytes
    temp = NextMersenneWord();
    switch (size%4)
    {
        case 3: output[2] = CRYPTOPP_GET_BYTE_AS_BYTE(temp, 1); /* fall through */
        case 2: output[1] = CRYPTOPP_GET_BYTE_AS_BYTE(temp, 2); /* fall through */
        case 1: output[0] = CRYPTOPP_GET_BYTE_AS_BYTE(temp, 3); break;
        default: ;;
    }
}

The gyrations around CRYPTOPP_ALLOW_UNALIGNED_DATA_ACCESS, IS_BIG_ENDIAN and IS_LITTLE_ENDIAN can be hidden behind GetWord, if desired. Its an ordinary named function, but a lot goes on behind the scenes with its use of ConditionalByteReverse. See misc.h for details.

DiscardBytes

DiscardBytes rounds up to a word size, and then discards the requisite number of words.

void DiscardBytes(size_t n)
{
    for(size_t i=0; i < (n+3)/4; i++)
        NextMersenneWord();
}

Downloads

mersenne.zip - Mersenne Twister implementation for Crypto++.