\ 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 Macfs 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 initfed 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 wonft 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 machinesc 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 [itfs supposed to be non-linear] $ FFFFFFFF and \ For 32-bit machinesc 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 0cNN-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