#Date: Tue, 27 Apr 1999 07:53:33 +0100 (BST)
#From: Prof Brian D Ripley <ripley@stats.ox.ac.uk>
#To: Paul Gilbert <pgilbert@bank-banque-canada.ca>
#Subject: Re: [R] random sequence
#
#On Mon, 26 Apr 1999, Paul Gilbert wrote:
#
#> Brian
#> 
#> >What exactly do you want?  My class exercise is to implement the
#> >Wichmann-Hill algorithm in C/Fortran and use it with S-PLUS! But
#> >while one is at it, why not use a decent generator with both?
#> 
#> It would certainly be nice to have all the R generators in S, but for now I
#> would be happy to have just one. What I would like is to be able to run the
#> programs which test my code, and get the same random sequences so that I 
#> don't need different check values for R and Splus.
#
#>From my files (from before R was thought of):
#

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
}

#> assign(".WHseed",  1:3, frame=0)
#> WichHill(10)
# [1] 0.03381877 0.77754189 0.05273525 0.74462407 0.49036219 0.98285437
# [7] 0.80915099 0.71338138 0.80102091 0.98958603
#> WichHill(10)
# [1] 0.9168563 0.4666467 0.6701995 0.9274966 0.8431496 0.7882993 0.7696486
# [8] 0.2939858 0.6187752 0.9892436
#
#in R:
#
#> .Random.seed <- c(0, 1:3)  # select Wichmann-Hill, seed = 1:3
#> runif(10)
# [1] 0.03381877 0.77754189 0.05273525 0.74462407 0.49036219 0.98285437
# [7] 0.80915099 0.71338138 0.80102091 0.98958603
#> runif(10)
# [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.