mardi 28 mars 2017

Mersenne Twister in python

In an attempt to better understand the Mersenne Twister (MT) I decided to convert it into python. Here you can run the code and see the commented C source code at the bottom.

The code runs but I'm not getting the expected output. Based on the C output here The first numbers I should be getting are 1067595299 955945823 477289528 4107218783 4228976476 but I'm getting 2594155707 3648022958 667752806 3967333545 3535325493.


Side note:

I found another implementation on SO saying they got the same output as C here but when I run it in python2 or python3 (xrange->range and no L in hexadecimal) the output also doesn't match Cs expected, it gives 0.817330 0.999061 0.510354 0.131533 0.035416 instead of the expected 0.76275443 0.99000644 0.98670464 0.10143112 0.27933125. ... If you're wondering changing my code to give a result in [0,1) interval gives 0.6039989429991692 0.8493715333752334 0.15547331562265754 0.9237168228719383 0.823132110061124 which is also different from both of these results.


Why does this python code not give the same output as the C code? Did I incorrectly implement something?

Code:

class RandomGenerators:
  # period parameters
  N = 624
  M = 394
  MATRIX_A = 0x9908b0df    # constant vector a
  UPPER_MASK = 0x80000000  # most significant w-r bits
  LOWER_MASK = 0x7fffffff  # least significant r bits

  def __init__(self):
    self.mt = [None]
    self.mti = self.N + 1

  '''
  initialize by an array with array-length
  init_key is the array for initializing keys
  key_length is its length
  slight change for C++, 2004/2/26
  '''
  def init_by_array(self, init_key, key_length):
    #int i, j, k;
    self.init_genrand(19650218)
    i,j = 1,0
    k = self.N if self.N > key_length else key_length
    #k = (self.N > key_length ? self.N : key_length)
    #for (; k; k--) {
    while k > 0:
      self.mt[i] = (self.mt[i] ^ ((self.mt[i-1] ^ (self.mt[i-1] >> 30)) * 1664525)) + init_key[j] + j # non linear 
      self.mt[i] &= 0xffffffff # for WORDSIZE > 32 machines 
      i+=1
      j+=1
      if i >= self.N:
        self.mt[0] = self.mt[self.N-1]
        i = 1
      if j >= key_length:
        j = 0
      k-=1

    k = self.N-1
    #for (k=N-1; k; k--) {
    while k > 0:
      self.mt[i] = (self.mt[i] ^ ((self.mt[i-1] ^ (self.mt[i-1] >> 30)) * 1566083941)) - i # non linear 
      self.mt[i] &= 0xffffffff # for WORDSIZE > 32 machines 
      i += 1
      if i >= self.N:
        self.mt[0] = self.mt[self.N-1]
        i=1
      k-=1

    self.mt[0] = 0x80000000 # MSB is 1; assuring non-zero initial array

  # initializes mt[N] with a seed
  def init_genrand(self,s):
    self.mt[0] = s & 0xffffffff
    self.mti = 1
    while self.mti < self.N:
      self.mt.append(1812433253 * (self.mt[self.mti-1] ^ (self.mt[self.mti - 1] >> 30)) + self.mti)
      self.mt[self.mti] &= 0xffffffff
      self.mti += 1

  # generates a random number on [0,0xffffffff]-interval
  def genrand_int32(self):
    y = 0
    mag01=[0x0, self.MATRIX_A]
    if self.mti >= self.N: # generate N words at one time
      kk = 0
      if self.mti == self.N + 1: # if init_genrand() has not been called,
        self.init_genrand(5489) # a default initial seed is used

      while kk < self.N - self.M:
        y = (self.mt[kk] & self.UPPER_MASK) | (self.mt[kk+1] & self.LOWER_MASK)
        self.mt[kk] = self.mt[kk+self.M] ^ (y >> 1) ^ mag01[y & 0x1]
        kk += 1

      while kk < self.N - 1:
        y = (self.mt[kk] & self.UPPER_MASK) | (self.mt[kk+1] & self.LOWER_MASK)
        self.mt[kk] = self.mt[kk+(self.M-self.N)] ^ (y >> 1) ^ mag01[y & 0x1]
        kk += 1

      y = (self.mt[self.N-1] & self.UPPER_MASK) | (self.mt[0] & self.LOWER_MASK)
      self.mt[self.N-1] = self.mt[self.M-1] ^ (y >> 1) ^ mag01[y & 0x1]

      self.mti = 0

    y = self.mt[self.mti]
    self.mti += 1

    # Tempering 
    y ^= (y >> 11)
    y ^= (y << 7) & 0x9d2c5680
    y ^= (y << 15) & 0xefc60000
    y ^= (y >> 18)

    return y

  # generates a random number on [0,0x7fffffff]-interval 
  def genrand_int31(self):
    #return (long)(genrand_int32() >> 1)
    return (self.genrand_int32() >> 1)

  # generates a random number on [0,1]-real-interval 
  def genrand_real1(self):
    return self.genrand_int32() * (1.0 / 4294967295.0)
    # divided by 2^32-1 

  # generates a random number on [0,1)-real-interval 
  def genrand_real2(self):
    return self.genrand_int32() * (1.0 / 4294967296.0)
    # divided by 2^32 

  # generates a random number on (0,1)-real-interval 
  def genrand_real3(self):
    return ((self.genrand_int32()) + 0.5) * (1.0 / 4294967296.0)
    # divided by 2^32 

  # generates a random number on [0,1) with 53-bit resolution
  def genrand_res53(self):
    a, b = self.genrand_int32()>>5, self.genrand_int32()>>6 
    return(a * 67108864.0 + b) * (1.0 / 9007199254740992.0)

r = RandomGenerators()
init = [0x123, 0x234, 0x345, 0x456]
length = 4
r.init_by_array(init, length)
for i in range(10):
  print('mt: '+str(r.genrand_int32()))




Aucun commentaire:

Enregistrer un commentaire