mtrand.cpp

// mtrand.cpp, see include file mtrand.h for information

#include "mtrand.h"

namespace Util 
{

   /*
   * Default constructor.
   *
   * Uses default seed only if this is the first instance
   */
   MTRand_int32::MTRand_int32() 
    : p(0),
      init(false)
   { 
      for (int i=0; i< n; ++i) state[i] = 0x0UL;
      if (!init) seed(5489UL); 
      init = true; 
   }
 
   /*
   * Constructor with 32 bit int as seed.
   */
   MTRand_int32::MTRand_int32(unsigned long s) 
    : p(0),
      init(false)
   { 
      for (int i=0; i< n; ++i) state[i] = 0x0UL;
      seed(s); 
      init = true; 
   }
 
   /*
   * Constructor with array of size 32 bit ints as seed
   */
   MTRand_int32::MTRand_int32(const unsigned long* array, int size) 
    : p(0),
      init(false)
   { 
      for (int i=0; i< n; ++i) state[i] = 0x0UL;
      seed(array, size); 
      init = true; 
   }

   /*
   * Copy constructor.
   */
   MTRand_int32::MTRand_int32(const MTRand_int32& other) 
    : p(other.p),
      init(other.init)
   { 
      for (int i=0; i< n; ++i) {
         state[i] = other.state[i];
      }
   }
 
   /*
   * Assignment constructor.
   */
   void MTRand_int32::operator = (const MTRand_int32& other) 
   { 
      p = other.p;
      init = other.init;
      for (int i=0; i< n; ++i) {
         state[i] = other.state[i];
      }
   }

   /* 
   * Generate new state vector.
   */
   void MTRand_int32::gen_state() 
   { 
     for (int i = 0; i < (n - m); ++i)
       state[i] = state[i + m] ^ twiddle(state[i], state[i + 1]);
     for (int i = n - m; i < (n - 1); ++i)
       state[i] = state[i + m - n] ^ twiddle(state[i], state[i + 1]);
     state[n - 1] = state[m - 1] ^ twiddle(state[n - 1], state[0]);
     p = 0; // reset position
   }

   /*
   * Init by 32 bit seed
   */
   void MTRand_int32::seed(unsigned long s) {  
     state[0] = s & 0xFFFFFFFFUL; // for > 32 bit machines
     for (int i = 1; i < n; ++i) {
       state[i] = 1812433253UL * (state[i - 1] ^ (state[i - 1] >> 30)) + i;
   // see Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier
   // in the previous versions, MSBs of the seed affect only MSBs of the array state
   // 2002/01/09 modified by Makoto Matsumoto
       state[i] &= 0xFFFFFFFFUL; // for > 32 bit machines
     }
     p = n; // force gen_state() to be called for next random number
   }
   
   void MTRand_int32::seed(const unsigned long* array, int size) { // init by array
     seed(19650218UL);
     int i = 1, j = 0;
     for (int k = ((n > size) ? n : size); k; --k) {
       state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525UL))
         + array[j] + j; // non linear
       state[i] &= 0xFFFFFFFFUL; // for > 32 bit machines
       ++j; j %= size;
       if ((++i) == n) { state[0] = state[n - 1]; i = 1; }
     }
     for (int k = n - 1; k; --k) {
       state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941UL)) - i;
       state[i] &= 0xFFFFFFFFUL; // for > 32 bit machines
       if ((++i) == n) { state[0] = state[n - 1]; i = 1; }
     }
     state[0] = 0x80000000UL; // MSB is 1; assuring non-zero initial array
     p = n; // force gen_state() to be called for next random number
   }

}