## George Marsaglia's PRNG

AM (21 june 2007) I came across a relatively simple PRNG recommended by George Marsaglia, arguably an authority on the subject:

``` # marsaglia.tcl --
#     Implementation of a PRNG according to George Marsaglia
#
namespace eval ::PRNG {
variable mod [expr {wide(256)*wide(256)*wide(256)*wide(256)-5}]
variable fac [expr {wide(256)*wide(32)}]
variable x1 [expr {wide(\$mod*rand())}]
variable x2 [expr {wide(\$mod*rand())}]
variable x3 [expr {wide(\$mod*rand())}]

puts \$mod
}

proc ::PRNG::marsaglia {} {
variable mod
variable fac
variable x1
variable x2
variable x3

set xn [expr {(\$fac*(\$x3+\$x2+\$x1))%\$mod}]

set x1 \$x2
set x2 \$x3
set x3 \$xn

return [expr {\$xn/double(\$mod)}]
}

# main --
#     A small test
#
package require Tk
package require Plotchart

pack [canvas .c]
set p [::Plotchart::createXYPlot .c {0 1 0.2} {0 1 0.2}]
set x [::PRNG::marsaglia]
\$p dataconfig series -type symbol -symbol dot
for {set i 0} {\$i < 1000} {incr i} {
set y [::PRNG::marsaglia]
\$p plot series \$x \$y
set x \$y
}```

The implementation has a few disadvantages: using the builtin random number generator for seeding it means you only get 2**31 different series, but the run-length is about 2**92 for each series.

AM (22 june 2007) I am thinking about a small package/module to make Monte Carlo simulations easy to do.

Here is my first attempt at such simulations: A simple Monte Carlo simulation

And here is a first attempt at a package that could be useful:

``` # random.tcl --
#     Create procedures that return various types of pseudo-random
#     number generators (PRNGs)
#

# math::random --
#     Create the namespace
#
namespace eval ::math::random {
variable count 0
variable pi    [expr {4.0*atan(1.0)}]
}

# prng_Binomial --
#     Create a PRNG with a binomial distribution
#
# Arguments:
#     p         Probability that the outcome will be 1
#
# Result:
#     Name of a procedure that returns a binomially distributed
#     random number
#
proc ::math::random::prng_Binomial {p} {
variable count

incr count

set name ::math::random::PRNG_\$count
proc \$name {} [string map [list P \$p] {return [expr {rand()<P? 1 : 0}]}]

return \$name
}

# prng_Uniform --
#     Create a PRNG with a uniform distribution in a given range
#
# Arguments:
#     min       Minimum value
#     max       Maximum value
#
# Result:
#     Name of a procedure that returns a uniformly distributed
#     random number
#
proc ::math::random::prng_Uniform {min max} {
variable count

incr count

set name ::math::random::PRNG_\$count
set range [expr {\$max-\$min}]
proc \$name {} [string map [list MIN \$min RANGE \$range] {return [expr {MIN+RANGE*rand()}]}]

return \$name
}

# prng_Exponential --
#     Create a PRNG with an exponential distribution with given mean
#
# Arguments:
#     mean      Mean value
#
# Result:
#     Name of a procedure that returns a uniformly distributed
#     random number
#
proc ::math::random::prng_Uniform {min max} {
variable count

incr count

set name ::math::random::PRNG_\$count
proc \$name {} [string map [list MEAN \$mean] {return [expr {-MEAN*log(rand())}]}]

return \$name
}

# prng_Discrete --
#     Create a PRNG with a uniform but discrete distribution
#
# Arguments:
#     n         Outcome is an integer between 0 and n-1
#
# Result:
#     Name of a procedure that returns such a random number
#
proc ::math::random::prng_Discrete {n} {
variable count

incr count

set name ::math::random::PRNG_\$count
proc \$name {} [string map [list N \$n] {return [expr {int(N*rand())}]}]

return \$name
}

# prng_Normal --
#     Create a PRNG with a normal distribution
#
# Arguments:
#     mean      Mean of the distribution
#     stdev     Standard deviation of the distribution
#
# Result:
#     Name of a procedure that returns such a random number
#
# Note:
#     Use the Box-Mueller method to generate a normal random number
#
proc ::math::random::prng_Normal {mean stdev} {
variable count

incr count

set name ::math::random::PRNG_\$count
proc \$name {} [string map [list MEAN \$mean STDEV \$stdev] \
{
variable pi
set rad [expr {sqrt(-log(rand()))}]
set phi [expr {2.0*\$pi*rand()}]
set r   [expr {\$rad*cos(\$phi)}]
return [expr {MEAN + STDEV*\$r}]
}]

return \$name
}

# prng_Disk --
#     Create a PRNG with a uniform distribution of points on a disk
#
# Arguments:
#
# Result:
#     Name of a procedure that returns the x- and y-coordinates of
#     such a random point
#
proc ::math::random::prng_Disk {rad} {
variable count

incr count

set name ::math::random::PRNG_\$count
proc \$name {} [string map [list RAD \$rad] \
{
variable pi
set phi [expr {2.0*\$pi*rand()}]
set x   [expr {\$rad*cos(\$phi)}]
set y   [expr {\$rad*sin(\$phi)}]
return [list \$x \$y]
}]

return \$name
}

# main --
#     Test code
#
set bin [::math::random::prng_Binomial 0.2]

set ones  0
set zeros 0
for { set i 0} {\$i < 100000} {incr i} {
if { [\$bin] } {
incr ones
} else {
incr zeros
}
}

puts "Binomial: \$ones - \$zeros"

set discrete [::math::random::prng_Discrete 10]

for { set i 0} {\$i < 100000} {incr i} {
set v [\$discrete]

if { [info exists count(\$v)] } {
incr count(\$v)
} else {
set count(\$v) 1
}
}

puts "Discrete:"
parray count

set rect [::math::random::prng_Rectangle 10 3]

puts "Rectangle:"
for { set i 0} {\$i < 10} {incr i} {
puts [\$rect]
}

set normal [::math::random::prng_Normal 0 1]

puts "Normal:"
for { set i 0} {\$i < 10} {incr i} {
puts [\$normal]
}```

The idea is that each random variable will be associated with its own command. That way there can be no mixup of distributions or parameters, as there is a command specific for that variable.

 Category Mathematics