1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
|
# this code is a version of the mersenne twister random number generator which
# is supposed to be used from RPython without the Python interpreter wrapping
# machinery etc.
# this is stolen from CPython's _randommodule.c
from rpython.rlib.rarithmetic import r_uint, intmask
N = 624
M = 397
MATRIX_A = r_uint(0x9908b0df) # constant vector a
UPPER_MASK = r_uint(0x80000000) # most significant w-r bits
LOWER_MASK = r_uint(0x7fffffff) # least significant r bits
MASK_32 = r_uint(0xffffffff)
TEMPERING_MASK_A = r_uint(0x9d2c5680)
TEMPERING_MASK_B = r_uint(0xefc60000)
MAGIC_CONSTANT_A = r_uint(1812433253)
MAGIC_CONSTANT_B = r_uint(19650218)
MAGIC_CONSTANT_C = r_uint(1664525)
MAGIC_CONSTANT_D = r_uint(1566083941)
class Random(object):
def __init__(self, seed=r_uint(0)):
self.state = [r_uint(0)] * N
self.index = 0
self.init_genrand(seed)
def init_genrand(self, s):
mt = self.state
mt[0]= s & MASK_32
for mti in range(1, N):
mt[mti] = (MAGIC_CONSTANT_A *
(mt[mti - 1] ^ (mt[mti - 1] >> 30)) + r_uint(mti))
# See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
# In the previous versions, MSBs of the seed affect
# only MSBs of the array mt[].
# for >32 bit machines
mt[mti] &= MASK_32
self.index = N
def init_by_array(self, init_key):
key_length = len(init_key)
mt = self.state
self.init_genrand(MAGIC_CONSTANT_B)
i = 1
j = 0
if N > key_length:
max_k = N
else:
max_k = key_length
for k in range(max_k, 0, -1):
mt[i] = ((mt[i] ^
((mt[i - 1] ^ (mt[i - 1] >> 30)) * MAGIC_CONSTANT_C))
+ init_key[j] + r_uint(j)) # non linear
mt[i] &= MASK_32 # for WORDSIZE > 32 machines
i += 1
j += 1
if i >= N:
mt[0] = mt[N - 1]
i = 1
if j >= key_length:
j = 0
for k in range(N - 1, 0, -1):
mt[i] = ((mt[i] ^
((mt[i - 1] ^ (mt[i - 1] >> 30)) * MAGIC_CONSTANT_D))
- i) # non linear
mt[i] &= MASK_32 # for WORDSIZE > 32 machines
i += 1
if (i>=N):
mt[0] = mt[N - 1]
i = 1
mt[0] = UPPER_MASK
def _conditionally_apply(self, val, y):
if y & r_uint(1):
return val ^ MATRIX_A
return val
def genrand32(self):
mt = self.state
if self.index >= N:
for kk in range(N - M):
y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK)
mt[kk] = self._conditionally_apply(mt[kk+M] ^ (y >> 1), y)
for kk in range(N - M, N - 1):
y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK)
mt[kk] = self._conditionally_apply(mt[kk + (M - N)] ^ (y >> 1),
y)
y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK)
mt[N - 1] = self._conditionally_apply(mt[M - 1] ^ (y >> 1), y)
self.index = 0
y = mt[self.index]
self.index += 1
y ^= y >> 11
y ^= (y << 7) & TEMPERING_MASK_A
y ^= (y << 15) & TEMPERING_MASK_B
y ^= (y >> 18)
return y
def random(self):
a = intmask(self.genrand32() >> 5)
b = intmask(self.genrand32() >> 6)
return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0)
def jumpahead(self, n):
mt = self.state
for i in range(N - 1, 1, -1):
j = n % i
mt[i], mt[j] = mt[j], mt[i]
nonzero = False
for i in range(1, N):
mt[i] += r_uint(i + 1)
mt[i] &= r_uint(0xffffffff)
nonzero |= bool(mt[i])
# Ensure the state is nonzero: in the unlikely event that mt[1] through
# mt[N-1] are all zero, set the MSB of mt[0] (see issue #14591). In the
# normal case, we fall back to the pre-issue 14591 behaviour for mt[0].
if nonzero:
mt[0] += r_uint(1)
mt[0] &= r_uint(0xffffffff)
else:
mt[0] = r_uint(0x80000000)
self.index = N
|
