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()))
--EDIT--
mt[] values are not getting initialized correctly. The first few values in C after calling init_by_array are:
-2147483648
1827812183
1371430253
-735590895
In python the first few are:
2147483648
1827812183
1371430253
3559376401
So it seems the issue is with positive/negative hexadecimals. I modified the init_by_array here the code below but I'm guessing genrand_int32 has similar issues (I've tried adding similar code at the end and throughout without fixing the problem):
for i,num in enumerate(self.mt):
if num>=2147483648: # (2**32)/2
self.mt[i]-=4294967296 # 2**32
