\ Mersenne Twister
\ Random number generator
\ Translated into Mops on the Mac by Michael Kuyumcu
\ July/August 2005 - v.0.5
\ freeware
\ Original disclaimer to be included:
(*
A C-program for MT19937, with initialization improved 2002/1/26.
Coded by Takuji Nishimura and Makoto Matsumoto.
Before using, initialize the state by using init_genrand(seed)
or init_by_array(init_key, key_length).
Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above
copyright
notice, this list of conditions and the following disclaimer in
the
documentation and/or other materials provided with the
distribution.
3. The names of its contributors may not be used to endorse or
promote
products derived from this software without specific prior
written
permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
COPYRIGHT OWNER OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
Any feedback is very welcome.
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
*)
decimal
need zFloating point
:class MT19937 super{ Object }
record{
var N \ number of 4-byte integers used for the state vector
var M \ ?
var matrix_a \ constant âvectorá a
var upper_mask \ most significant w-r bits
var lower_mask \ least significant r bits
var mti
624 array mt \ If the algorithm changes, the 624 might be changed.
\ If I only knew how to write the code so that the 624
\ has to appear (and be changed) only once...
var ii \ local variables of manifold use...
var jj
var kk
var y \ used to calculate and hold integer random numbers
var seedsize \ The array length to be used.
var arraybase \ The location of the array in the Macfs memory.
\ Not the address of the 8-byte descriptor, but the start
\ of the actual array data
}
\ INIT: is the first thing to be called after the instantiation of one
\ MT19937 Mersenne Twister random number generator object.
\ Code example: mt19937 twister init: twister
:m INIT:
624 ?N \ For revised editions of the algorithm, these two
397 ?M \ numbers might be changed - here.
$ 9908B0DF ?matrix_a
$ 80000000 ?upper_mask
$ 7FFFFFFF ?lower_mask
N? 1+ ?mti \ this means that the generator has not been seeded
yet.
;m
\ Here the random number generator is seeded with an user-defined
number (seed)
\ passed on the data stack. Seeding the generator with a 4-byte
integer scalar
\ like here is only ONE way of initializing it. It can be initfed
alternatively
\ by a whole array of seed values. See the method below, INIT_BY_ARRAY
\ Usage: mt19937 twister init: twister
\ init_genrand: twister
:m INIT_GENRAND: ( seed -- )
$ FFFFFFFF AND \ this AND is for machines with register widths
0 to: mt \ greater than 32 bits. The algorithm is defined
\ only for register widths of 32 bits.
N? 1 DO \ After the first state vector element, we calculate
\ the remaining elements in a loop from 1 to 623 (N-1).
i 1- at: mt \ Retrieve the state vector element calculated
before.
dup \ Put another copy of it on the data stack
30 rshift \ Shift logically 30 bits to the right
xor \ XOR with the state vector value calculated before
1812433253 * \ Multiply with this âmagicá number determined
\ mathematically by the authors of the algorithm
i + \ and add the loop index variable.
$ FFFFFFFF and \ this AND is for machines with register widths
\ greater than 32 bits (see above explanation).
i to: mt \ Last, store the current calculated state vector
value.
LOOP
N? ?mti \ Signalling: The GENRAND routine has been called
\ meaning: the generator has been initialized
\ by the user (so the program wonft do it again).
;m
\ Here the random number generator is seeded with an user-defined
ARRAY of
\ seed (key) values. The array must be created by the user. Its
address must
\ be on the data stack, and above that must be the number of cells to
be used
\ (which will normally be its size, the number of cells it holds).
\ Usage: mt19937 twister init: twister
\ init_by_array: twister
:m INIT_BY_ARRAY: ( seed_array_addr seed_array_size -- )
?seedsize \ Store the values on the stack for later
4+ 4+ ?arraybase \ reference for calculations
19650218 INIT_GENRAND: self \ Initialize random generator engine
1 ?ii
0 ?jj
seedsize? N? max ?kk \ If N is greater than the seed_array_size,
\ take N as the loop limit, otherwise the size.
\ First mixup phase
\ ." First phase mixup." kk? . cr \ key drop \ debugging
1 kk? 1- DO
." " \ This is the strangest line of the code:
\ if one omits it, PowerMops hangs.
ii? 1- at: mt \ Retrieve state vector element mt[i-1]
dup \ Duplicate it for the subsequent XOR-ing
30 rshift \ Shift the uppermost copy 30 bits right
xor \ XOR the result with the very vector element
1664525 * \ Multiply by this magic number
ii? at: mt \ Retrieve the recent vector element
xor \ XOR them together
\ key drop \ check the stack, debugging
arraybase? \ Access to the seed key array
jj? 4 * + \ jj * 4 [jj is the index of the array el.]
@ + \ This gives the value of the key array element
jj? +
$ FFFFFFFF and \ For 32-bit machinesc
ii? to: mt \ Now store the calculated vector element
ii? 1+ ?ii \ Increment ii and jj by one
jj? 1+ ?jj
ii? N? >= IF
N? 1- at: mt
0 to: mt \ If j>=N, then mt[0|:=mt[N-1]
1 ?ii \ c and ii:=1
THEN
jj? seedsize? >= IF \ If jj >= seedsize then set jj:=0
0 ?jj
THEN
\ key \ for debugging
\ 65 = IF UNLOOP EXIT THEN \ for debugging
-1 +LOOP
\ Second mixup phase
\ ." Second phase mixup" cr key drop \ debugging
1 N? 1- DO
." " \ Strange that this is necessary --- ???
ii? 1- at: mt \ Retrieve state vector element mt[i-1]
dup \ Duplicate it for the subsequent XOR-ing
30 rshift \ Shift the uppermost copy 30 bits right
xor \ XOR the result with the very vector element
1566083941 * \ Multiply by this magic number
ii? at: mt \ Retrieve the recent vector element
xor \ XOR them together
ii? - \ and subtract ii [itfs supposed to be non-linear]
$ FFFFFFFF and \ For 32-bit machinesc
ii? to: mt \ Now store the mixed-up state vector
ii? N? >= IF
N? 1- at: mt
0 to: mt \ If j>=N, then mt[0|:=mt[N-1]
1 ?ii \ c and ii:=1
THEN
-1 +LOOP
\ Last phase: assign mt[0] a new contents (non-zero it must be)
$ 80000000 0 to: mt
;m
\
\ Now here is the main routine for generating integer pseudo-random
\ numbers. It gives 32-bit integer numbers. There are other routines
\ that return 31-bit integers, and ârealá numbers for different
intervals
\ namely a ârealá out of [0,1], [0,1), & ârealá numbers with 53-bit
resolution
\ out of the interval [0,1). Of course, these are no real ârealá
numbers,
\ but a subset of the rational numbers of the internal floating-point
\ unit. One calls all of these routines without passing anything on
the stack.
\ They leave the desired number (of integer or ârealá format on the
stack).
\
:m GENRAND_INT32: ( -- pseudo random number on [0,$FFFFFFFF])
mti? N? >= \ generate N 32-bit integers at a time
IF
mti? N? 1+ = \ If mti = N+1 (meaning: the generator
IF \ has not been initialized by the user
5489 INIT_GENRAND: self \ so we initialize it here.
THEN
N? M? - 0 \ 1st part of the treatment
DO \ for elements mt[0] to mt[N-M-1]
i at: mt
upper_mask? and \ upper_mask and lower_mask are
i 1+ at: mt \ constant masks that blend out
lower_mask? and \ lower and upper 16-bit of a value
or \ both results are OR-ed together
dup \ we must shift this preliminary result y
$ 1 and \ The result of this AND is 0 or 1
matrix_a? * \ So the result of this is 0 or matrix_a
swap \ y should be still on top of the stack
1 rshift \ and to be shifted 1 bit to the right
xor
M? i + at: mt \ get random number mt[i+M]
xor \ Connect the result with it logically
i to: mt \ Our new random vector is stored.
LOOP
N? 1- N? M? - \ 2nd part of the treatment
DO \ for elements mt[N-M] to mt[N-2]
i at: mt \ The same as in the loop above
upper_mask? and \ upper_mask and lower_mask are
i 1+ at: mt \ constant masks that blend out
lower_mask? and \ lower and upper 16-bit of a value
or \ both results are OR-ed together
dup \ we must shift this preliminary result y
$ 1 and \ The result of this AND is 0 or 1
matrix_a? * \ So the result of this is 0 or matrix_a
swap \ y is to be still on top of the stack
1 rshift \ and to be shifted 1 bit to the right
xor
M? N? - i + at: mt \ Getting random number mt[i+(M-N)]
xor
i to: mt \ The new random vector is calculated.
LOOP
N? 1- at: mt \ 3rd and last part of the treatment
upper_mask? and \ only for element mt[N-1]
0 at: mt
lower_mask? and
or \ Now y, preliminary result, on stack
dup
$ 1 and
matrix_a? *
swap \ We need y on top of stack
1 rshift
xor
M? 1- at: mt
xor
N? 1- to: mt
0 ?mti \ reset mti for the next round of N el.
THEN
mti? at: mt \ Get current random number on stack, =y
mti? 1+ ?mti \ Increment the random number counter
\ Start of the final non-linear tempering
\ On the stack is the current result y
dup \ First step: duplicate y for the XORing
11 rshift xor \ y ^= (y >> 11) [in C notation]
dup \ Second step of the tempering
7 lshift \ y^= (y << 7) & 0x9D2C5680
$ 9D2C5680 and \ AND it with this constant
xor
dup \ Third tempering step
15 lshift \ y^= (y << 15) & 0xEFC60000
$ EFC60000 and \ AND it with this constant
xor
dup \ Fourth and last step of tempering
18 rshift \ y^= (y >> 18)
xor \ Now y, the result, is on the stack;
\ y is the asked-for ârandomá number.
;m
\ Here follow a few helping routines that all call the genrand_init32
\ routine above. They map the obtained pseudo-random numbers onto
different
\ intervals.
:m GENRAND_INT31: ( -- 31-bit pseudo-random-number )
genrand_int32: self \ Get a 32-bit pseudo-random number
1 rshift \ divide it by 2 to get into [0,$7FFFFFFF]
;m
:m GENRAND_[0,1]: ( -- ârealá number from [0,1], 0 and 1 included )
genrand_int32: self \ Get a 32-bit integer pseudo-random number
s>f \ Convert this to a float
4294967295.0 f/ \ Divide by (2^32-1)
\ Now the number is in [-0.5,0.5]
0.5 f+
;m
:m GENRAND_[0,1): ( -- ârealá number from [0,1[ )
BEGIN
genrand_int32: self \ Get a 32-bit integer pseudo-random number
s>f \ Convert this to a float
4294967295.0 f/ \ Divide by (2^32-1)
\ Now the number is in [-0.5,0.5]
0.5 f+
fdup
1.0 f<> \ Produce new until <> 1.0
UNTIL
;m
:m GENRAND_(0,1): ( -- ârealá number from ]0,1[, 0 and 1 excluded )
genrand_int32: self \ Get a 32-bit integer pseudo-random number
s>f \ Convert this to a float
4294967296.0 f/ \ Divide by (2^32)
\ Now the number is in ]-0.5,0.5[
0.4999999999999999 f+
;m
:m GENRAND_(-0.5,0.5): ( -- ârealá number from ]-0.5,0.5[ )
genrand_int32: self \ Get a 32-bit integer pseudo-random number
s>f \ Convert this to a float
4294967296.0 f/ \ Divide by (2^32)
;m
:m GENRAND_[-0.5,0.5]: ( -- ârealá number from [-0.5,0.5] )
genrand_int32: self \ Get a 32-bit integer pseudo-random number
s>f \ Convert this to a float
4294967295.0 f/ \ Divide by (2^32-1)
;m
\ Now this is the most typical case: the user wants some integer number
\ between 0 (inclusive) and NN (exclusive).
:m GENRAND_[0,N): ( NN -- integer number from 0cNN-1 )
s>f \ Convert the max border to a float
genrand_[0,1): self \ Get a float pseudo-random number
f* \ Now the float stack holds some value
\ between 0 (included) and NN (excluded)
f>s \ Retranslate into an integer
;m
;class