Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Which Mersenne Twister does C++11 provide?

I'm having trouble determining which variant of Mersenne Twister C++11 provides. Looking at Matsumoto and Nishimura ACM paper at Mersenne twister: A 623 Dimensionally Equidistributed Uniform Pseudorandom Number Generator, the authors provide the algorithm, an implementation of the algorithm, and call it MT19937.

However, when I test C++11's same-named generator with the small program below, I cannot reproduce the stream created by Matsumoto and Nishimura's MT19937. The streams differ from the very first 32-bit word produced.

Which Mersenne Twister does C++11 provide?


The program below was run on Fedora 22 using GCC, -std=c++11 and GNU's stdlibc++.

std::mt19937 prng(102013);
for (unsigned int i = 0; i <= 625; i++)
{
    cout << std::hex << prng();

    if(i+1 != 625)
        cout << ",";

    if(i && i%8 == 0)
        cout << endl;
}
like image 305
jww Avatar asked Oct 08 '15 16:10

jww


People also ask

What is the Mersenne Twister algorithm?

The Mersenne Twister is a strong pseudo-random number generator. In non-rigorous terms, a strong PRNG has a long period (how many values it generates before repeating itself) and a statistically uniform distribution of values (bits 0 and 1 are equally likely to appear regardless of previous values).

Is Mersenne Twister random?

Implementations generally create random numbers faster than hardware-implemented methods. A study found that the Mersenne Twister creates 64-bit floating point random numbers approximately twenty times faster than the hardware-implemented, processor-based RDRAND instruction set.

What is mt19937 C++?

std::mt19937(since C++11) class is a very efficient pseudo-random number generator and is defined in a random header file. It produces 32-bit pseudo-random numbers using the well-known and popular algorithm named Mersenne twister algorithm.

What is the maximum number that you can get from the Mersenne Twister PRNG?

CUDA's implementation of the Mersenne Twister ( MT ) random number generator is limited to a maximal number of threads/blocks of 256 and 200 blocks/grid, i.e. the maximal number of threads is 51200 .


3 Answers

Looking at the MT19937 from your the linked to paper and the MT19937 defined by the standard it looks like they are the same but an additional layer of tempering was added and an initialization multiplier

If we look at the values defined by [rand.predef] 26.5.5(3) vs the parameters defined by the paper we have

32,624,397,31,0x9908b0df,11,0xffffffff,7,0x9d2c5680,15,0xefc60000,18,1812433253 <- standard
w ,n  ,m  ,r ,a         ,u ,d         ,s,b         ,t ,c         ,l ,f  
32,624,397,31,0x9908b0df,11,          ,7,0x9d2c5680,15,0xefc60000,18,           <- paper

This is where the difference is coming from. Also according to the standard the 10,000th iteration of std::mt19937 is 399268537

like image 91
NathanOliver Avatar answered Oct 27 '22 13:10

NathanOliver


It seems C++11 provides Mersenne Twister with improved initialization

I've just extracted original C implementation, and compared with C++.

#include <iostream>
#include <cstdio>
#include <random>

#define N 624
#define M 397
#define MATRIX_A 0x9908b0dfUL   /* constant vector a */
#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
#define LOWER_MASK 0x7fffffffUL /* least significant r bits */

static unsigned long mt[N]; /* the array for the state vector  */
static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */

void init_genrand(unsigned long s)
{
    mt[0]= s & 0xffffffffUL;
    for (mti=1; mti<N; mti++) {
        mt[mti] =
        (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
        mt[mti] &= 0xffffffffUL;
    }
}

unsigned long genrand_int32()
{
    unsigned long y;
    static unsigned long mag01[2]={0x0UL, MATRIX_A};

    if (mti >= N) { /* generate N words at one time */
        int kk;

        if (mti == N+1)   /* if init_genrand() has not been called, */
            init_genrand(5489UL); /* a default initial seed is used */

        for (kk=0;kk<N-M;kk++) {
            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
            mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
        }
        for (;kk<N-1;kk++) {
            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
            mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
        }
        y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
        mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];

        mti = 0;
    }

    y = mt[mti++];

    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);

    return y;
}

int main()
{
    init_genrand(102013);

    std::mt19937 prng(102013);

    for (size_t i = 0; i < 10000; ++i) {
       if (genrand_int32() != prng()) {
          std::cout << "ERROR" << std::endl;
          return 1;
       }
    }

    std::cout << "OK" << std::endl;
    return 0;
}
like image 40
Stas Avatar answered Oct 27 '22 12:10

Stas


I should point out that C++11 actually provides many Mersenne twisters via the template class:

template <class UIntType,
          size_t word_size,
          size_t state_size,
          size_t shift_size,
          size_t mask_bits,
          UIntType xor_mask,
          size_t tempering_u,
          UIntType tempering_d,
          size_t tempering_s,
          UIntType tempering_b,
          size_t tempering_t,
          UIntType tempering_c,
          size_t tempering_l,
          UIntType initialization_multiplier>
  class mersenne_twister_engine;

If anyone has the bravery to explore these levers and knobs... Of course there are the two standard instantiations of these:

using mt19937
  = mersenne_twister_engine<uint_fast32_t,
                            32,
                            624,
                            397,
                            31,
                            0x9908b0df,
                            11,
                            0xffffffff,
                            7,
                            0x9d2c5680,
                            15,
                            0xefc60000,
                            18,
                            1812433253>;

and the 64-bit version:

using mt19937_64
  = mersenne_twister_engine<uint_fast64_t,
                            64,
                            312,
                            156,
                            31,
                            0xb5026f5aa96619e9,
                            29,
                            0x5555555555555555,
                            17,
                            0x71d67fffeda60000,
                            37,
                            0xfff7eee000000000,
                            43,
                            6364136223846793005>;

I've thought it would be nice to provide a toolbox for checking the quality of RNGs so people could try new instantiations.

Here is a comparison of template parms:

32,624,397,31,        0x9908b0df,11,        0xffffffff,7 ,        0x9d2c5680,15,        0xefc60000,18,1812433253          <- std::mt19937
64,312,156,31,0xb5026f5aa96619e9,29,0x5555555555555555,17,0x71d67fffeda60000,37,0xfff7eee000000000,43,6364136223846793005 <- std::mt19937_64
w ,n  ,m  ,r ,a                 ,u ,d                 ,s ,b                 ,t ,c                 ,l ,f  
32,624,397,31,        0x9908b0df,11,                  ,7 ,        0x9d2c5680,15,        0xefc60000,18,                    <- paper

With thanks to @NathanOliver.

like image 22
emsr Avatar answered Oct 27 '22 12:10

emsr