OpenEuphoria

Mersenne Twister Random Generator

• Posted by Derek Parnell <ddparnell at bigpond.com> Sep 15, 2001
• 586 views

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 ------------