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