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

 # 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:
 #     rad       Radius of the disk
 #
 # 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 rad [expr {RAD*sqrt((rand())}]
         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]
 }