source [file join [info library] init.tcl]
# Here is a complimentary-multiply-with-carry RNG with k=4097 and a near-record period,
# more than 10^33000 times as long as that of the Mersenne Twister. (2^131104 vs. 2^19937)
set maxint [ expr { 2 ** 63 } ]
set rmaxint ${maxint}.0
set millnum 0
proc initmill { { name ""} { seed 0 } } {
if { $name eq "" } { set name mill$::millnum ; incr ::millnum }
upvar $name mill
if { $seed != 0 } { ::mathfunc::srand $seed }
set mill(Q) {}
for { set i 0 } { $i < 4096 } { incr i } {
set init [ expr { wide( rand() * $::rmaxint ) } ]
lappend mill(Q) $init
}
set mill(c) 362436
return $name
}
set prevmill ""
proc spin { { spinner "" } } {
if { $spinner eq "" } {
set spinner $::prevmill
if { $spinner eq "" } {
set spinner [ initmill ]
}
}
set ::prevmill $spinner
upvar $spinner mill
set t 0
set a 18782
set i 4095
set x 0
set r 0xfffffffe
set i [ expr { ($i + 1) & 4095 } ]
set Qsubi [ lindex $mill(Q) $i ]
set t [ expr { wide( $a * $Qsubi + $mill(c) ) } ]
set mill(c) [ expr { wide($t >> 32) } ]
set x [ expr { wide($t + $mill(c)) } ]
if {$x < $mill(c)} { incr x ; incr mill(c) }
lset mill(Q) $i [ expr { wide( $r - $x ) } ]
set result [ lindex $mill(Q) $i ].0
return [ expr { abs ($result / $::rmaxint) } ]
}
# ex: set test [ initmill ] ; set val [spin $test]To use, you need to create a new mill to generate random numbers.
set name [ initmill ]
or
initmill <name> <seed>
A bare call to initmill will return an identifier for the new mill. If you pass a name in as the first argument, the new mill will have that name. If you also include a seed, the new mill will be initialized with that seed.
Once you have the mill you can spin it for the next random number:
set random [ spin $name ]
If you don't specify a mill name, it reuses the most recently used mill.
The actual mill is just an array variable in the local scope, so it will disappear automatically along with any other local variables if your proc exits. Note that the name of the non-existent mill will be saved and (mis)used again if followed by another bare call of spin, so be careful to call spin with the appropriate mill id next time.