1. Mersenne Twister Random Generator
Hi,
here is an implementation of the Mersenne Twister (MT) algorithm. I'd
appreciate some people validating it. I converted it from a C program. The
tricky parts in the conversion were the bitwise operations and the
zero-based array references.
It got the C code from http://www.math.keio.ac.jp/~matumoto/emt.html and
there are other varieties there as well.
----
Derek.
--------------
-- A Euphoria program for MT19937: Real number version
-- rand_mt() generates one pseudorandom real number (double)
-- which is uniformly distributed on [0,1]-interval, for each
-- call. seedrand_mt(seed) set initial values to the working area
-- of 624 words. Before rand_mt(), seedrand_mt(seed) must be
-- called once. 'seed' is any 32-bit integer, however if 0 is supplied
-- the value of time(0) is used instead.
-- Integer generator is obtained by modifying two lines.
-- This Euphoria version is based on the 'C' version coded by
-- Takuji Nishimura, considering the suggestions by
-- Topher Cooper and Marc Rieffel in July-Aug. 1997.
-- This library is free software; you can redistribute it and/or
-- modify it under the terms of the GNU Library General Public
-- License as published by the Free Software Foundation; either
-- version 2 of the License, or (at your option) any later
-- version.
-- This library is distributed in the hope that it will be useful,
-- but WITHOUT ANY WARRANTY; without even the implied warranty of
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-- See the GNU Library General Public License for more details.
-- You should have received a copy of the GNU Library General
-- Public License along with this library; if not, write to the
-- Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
-- 02111-1307 USA
-- Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
-- Any feedback is very welcome. For any question, comments,
-- see http://www.math.keio.ac.jp/matumoto/emt.html or email
-- matumoto at math.keio.ac.jp
include machine.e
-- Period parameters
constant N = 624
constant M = 397
constant MATRIX_A = #9908B0DF -- constant vector a
constant UPPER_MASK = #80000000 -- most significant w-r bits
constant LOWER_MASK = #7FFFFFFF -- least significant r bits
-- Tempering parameters
constant TEMPERING_MASK_B = #9D2C5680
constant TEMPERING_MASK_C = #EFC60000
function shift_bits(atom a, integer distance)
sequence abits
abits = int_to_bits(a,32)
if distance < 0 then
distance = -distance
abits[distance+1 .. 32] = abits[1 .. 32-distance]
abits[1 .. distance] = repeat(0, distance)
else
abits[1 .. 32-distance] = abits[distance+1 .. 32]
abits[32 - distance + 1 .. 32] = repeat(0, distance)
end if
return bits_to_int(abits)
end function
sequence mt
mt = repeat(0,N) -- the array for the state vector
integer mti
mti = N + 2 -- mti=N+2 means mt[N] is not initialized
-- initializing the array with a NONZERO seed
global procedure seedrand_mt(atom seed)
atom temp
-- setting initial seeds to mt[N] using
-- the generator Line 25 of Table 1 in
-- [KNUTH 1981, The Art of Computer Programming
-- Vol. 2 (2nd Ed.), pp102]
if seed = 0 then
seed = floor(time())
end if
mt[1]= remainder(seed , #FFFFFFFF)
mti = 1
while mti < N do
mt[mti+1] = remainder((69069 * mt[mti]),#FFFFFFFF)
mti += 1
end while
end procedure
sequence mag01
global function rand_mt(integer range)
atom y
integer kk
mag01= {0, MATRIX_A}
if (mti >= N) then -- generate N words at one time
if (mti = N+2) then -- if seedrand_mt() has not been called,
seedrand_mt(4357) -- a default initial seed is used
end if
kk = 0
while kk < N-M do
y = or_bits(and_bits(mt[kk+1],UPPER_MASK),
and_bits(mt[kk+2],LOWER_MASK)
)
mt[kk+1] = xor_bits(xor_bits(mt[kk+M+1] , shift_bits(y , 1)) ,
mag01[and_bits(y,1)+1])
kk += 1
end while
while kk<N-1 do
y = or_bits(and_bits(mt[kk+1],UPPER_MASK),
and_bits(mt[kk+2],LOWER_MASK)
)
mt[kk+1] = xor_bits(xor_bits(mt[kk+(M-N)+1] , shift_bits(y,1)) ,
mag01[and_bits(y,1)+1])
kk += 1
end while
y = or_bits(and_bits(mt[N],UPPER_MASK),
and_bits(mt[1],LOWER_MASK))
mt[N] = xor_bits(xor_bits(mt[M] , shift_bits(y,1)) ,
mag01[and_bits(y,1)+1])
mti = 0
end if
mti += 1
y = mt[mti]
y = xor_bits(y, shift_bits(y, 11))
y = xor_bits(y,and_bits(shift_bits(y,-7), TEMPERING_MASK_B))
y = xor_bits(y,and_bits(shift_bits(y,-15), TEMPERING_MASK_C))
y = xor_bits(y, shift_bits(y, 18))
return floor((y + #80000000) * (range/#FFFFFFFF) ) + 1
end function
----------- END OF SOURCE CODE ------------
