Updated 2008-01-04 22:43:23 by dkf

The Mersenne Twister is a (pseudo-)Random Number generator. The designers of the Mersenne Twister algorithm have their homepage at http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html

The random package also includes an implementation of the Mersenne Twister algorithm.

The following code implements the Mersenne Twister algorithm. It provides three public commands in the "mt" namespace:
srand seed
Initializes the algorithm with a random (binary) string as a seed value. The seed string can have any length (only the first 2496 bytes are significant)
int32
Returns a 32 bit integer pseudo-random number.
rand
Returns a floating-point random number in the interval [0..1) (i.e., the this command may return 0, but never returns 1).

 #
 # ----------------------------------------------------------------------
 # Mersenne Twister Random Number Generator
 #
 # Derived from the source code for MT19937 by Takuji Nishimura
 # and Makoto Matsumoto, which is available from their homepage
 # at http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
 #
 # Written by Frank Pilhofer. Released under BSD license.
 # ----------------------------------------------------------------------
 #

 namespace eval mt {
    variable N 624
    variable M 397
    variable MATRIX_A 0x9908b0df
    variable UPPER_MASK 0x80000000
    variable LOWER_MASK 0x7fffffff
    variable mag010 0
    variable mag011 $MATRIX_A
    variable mt
    variable mti
 }

 #
 # Initializes with a seed
 #

 proc mt::seed {s} {
    variable N
    variable mt
    variable mti

    set mt [list [expr {$s & 0xffffffff}]]
    set mtimm $mt

    for {set mti 1} {$mti < $N} {incr mti} {
	set t1 [expr {$mtimm ^ ($mtimm >> 30)}]
	set t2 [expr {1812433253 * $t1 + $mti}]
	set mtimm [expr {$t2 & 0xffffffff}]
	lappend mt $mtimm
    }
 }

 #
 # Initialize from a (binary) seed string
 #

 proc mt::init {s} {
    variable N
    variable mt

    seed 19650218

    set i 1
    set j 0

    #
    # The algorithm wants a list of 32 bit integers for the key
    #

    set slen [string length $s]

    if {($slen % 4) != 0} {
	append s [string repeat "\0" [expr {4-($slen%4)}]]
    }

    binary scan $s i* key

    if {$N > [llength $key]} {
	set k $N
    } else {
	set k [llength $key]
    }

    set mtimm [lindex $mt 0]

    for {} {$k} {incr k -1} {
	set keyj [lindex $key $j]
	set mti [lindex $mt $i]
	set t1 [expr {$mtimm ^ ($mtimm >> 30)}]
	set t2 [expr {$mti ^ ($t1 * 1664525)}]
	set t3 [expr {$t2 + $keyj + $j}]
	set mtimm [expr {$t3 & 0xffffffff}]
	set mt [lreplace $mt $i $i $mtimm]
	incr i
	incr j
	if {$i >= $N} {
	    set mt [lreplace $mt 0 0 $mtimm]
	    set i 1
	}
	if {$j >= [llength $key]} {
	    set j 0
	}
    }

    for {set k [expr {$N-1}]} {$k} {incr k -1} {
	set mti [lindex $mt $i]
	set t1 [expr {$mtimm ^ ($mtimm >> 30)}]
	set t2 [expr {$mti ^ ($t1 * 1566083941)}]
	set t3 [expr {$t2 - $i}]
	set mtimm [expr {$t3 & 0xffffffff}]
	set mt [lreplace $mt $i $i $mtimm]
	incr i
	if {$i >= $N} {
	    set mt [lreplace $mt 0 0 $mtimm]
	    set i 1
	}
    }

    set mt [lreplace $mt 0 0 0x80000000]
 }

 #
 # Produce some more random numbers
 #

 proc mt::more {} {
    variable N
    variable M
    variable mt
    variable mti
    variable MATRIX_A
    variable UPPER_MASK
    variable LOWER_MASK
    variable mag010
    variable mag011

    if {$mti == [expr {$N+1}]} {
	seed 5489
    }

    set newmt [list]

    for {set kk 0} {$kk<[expr {$N-$M}]} {incr kk} {
	set mtkk [lindex $mt $kk]
	set mtkkpp [lindex $mt [expr {$kk+1}]]
	set mtkkpm [lindex $mt [expr {$kk+$M}]]
	set y [expr {($mtkk & $UPPER_MASK) | ($mtkkpp & $LOWER_MASK)}]
	if {($y & 1) == 0} {
	    set mag01 $mag010
	} else {
	    set mag01 $mag011
	}
	set mtkk [expr {$mtkkpm ^ ($y >> 1) ^ $mag01}]
	lappend newmt $mtkk
    }
    for {} {$kk<[expr {$N-1}]} {incr kk} {
	set mtkk [lindex $mt $kk]
	set mtkkpp [lindex $mt [expr {$kk+1}]]
	set mtkkpm [lindex $newmt [expr {$kk+$M-$N}]]
	set y [expr {($mtkk & $UPPER_MASK) | ($mtkkpp & $LOWER_MASK)}]
	if {($y & 1) == 0} {
	    set mag01 $mag010
	} else {
	    set mag01 $mag011
	}
	set mtkk [expr {$mtkkpm ^ ($y >> 1) ^ $mag01}]
	lappend newmt $mtkk
    }
    set mtnm1 [lindex $mt [expr {$N-1}]]
    set mt0 [lindex $newmt 0]
    set mtmmm [lindex $newmt [expr {$M-1}]]
    set y [expr {($mtnm1 & $UPPER_MASK) | ($mt0 & $LOWER_MASK)}]
    if {($y & 1) == 0} {
	set mag01 $mag010
    } else {
	set mag01 $mag011
    }
    set mtkk [expr {$mtmmm ^ ($y >> 1) ^ $mag01}]
    lappend newmt $mtkk

    set mti 0
    set mt $newmt
 }

 #
 # ----------------------------------------------------------------------
 # Public interface
 # ----------------------------------------------------------------------
 #

 #
 # Initialize with a random (binary string) seed
 #

 proc mt::srand {seed} {
    init $seed
 }

 #
 # Generates an integer random number in the [0,0xffffffff] interval
 #

 proc mt::int32 {} {
    variable N
    variable mt
    variable mti

    if {$mti >= $N} {
	more
    }

    set y [lindex $mt $mti]
    incr mti

    set y [expr {$y ^  ($y >> 11)}]
    set y [expr {$y ^ (($y <<  7) & 0x9d2c5680)}]
    set y [expr {$y ^ (($y << 15) & 0xefc60000)}]
    set y [expr {$y ^  ($y >> 18)}]

    return [expr {$y & 0xffffffff}]
 }

 #
 # Generates a floating-point random number in the [0,1) interval
 #

 proc mt::rand {} {
    set i [int32]
    return [expr {double($i) / 4294967296.0}]
 }

 #
 # ----------------------------------------------------------------------
 # Print test vectors, for comparison with the original code
 # ----------------------------------------------------------------------
 #

 proc mt::test {} {
    srand [binary format i4 [list 0x123 0x234 0x345 0x456]]
    puts "1000 outputs of int32"
    for {set i 0} {$i < 1000} {incr i} {
	puts -nonewline [format "%10u " [int32]]
	if {($i % 5) == 4} {
	    puts ""
	}
    }
    puts "1000 outputs of real"
    for {set i 0} {$i < 1000} {incr i} {
	puts -nonewline [format "%10.8f " [rand]]
	if {($i % 5) == 4} {
	    puts ""
	}
    }
 }

Initially written by FPX.