Skip to content

random sequence

1 message · Brian Ripley

#
On Mon, 26 Apr 1999, Paul Gilbert wrote:

            
WichHill <- function(n)
{
        output <- numeric(n)
        seed <- get(".WHseed", frame=0)
        x <- seed[1]; y <- seed[2]; z <- seed[3]
        for(i in 1:n){
                x <- (171*x) %% 30269
                y <- (172*y) %% 30307
                z <- (170*z) %% 30323
                output[i] <- (x/30269 + y/30307 + z/30323) %% 1.0
        }
        assign(".WHseed", c(x,y,z), frame=0, immediate=T)
        output
}
[1] 0.03381877 0.77754189 0.05273525 0.74462407 0.49036219 0.98285437
 [7] 0.80915099 0.71338138 0.80102091 0.98958603
[1] 0.9168563 0.4666467 0.6701995 0.9274966 0.8431496 0.7882993 0.7696486
 [8] 0.2939858 0.6187752 0.9892436

in R:
[1] 0.03381877 0.77754189 0.05273525 0.74462407 0.49036219 0.98285437
 [7] 0.80915099 0.71338138 0.80102091 0.98958603
[1] 0.9168563 0.4666467 0.6701995 0.9274966 0.8431496 0.7882993 0.7696486
 [8] 0.2939858 0.6187752 0.9892436


You can do this faster in C/Fortran, of course, but the load/save of the 
seed is slow from compiled code and in pure S I get over 1000 unifs/sec.