# Generation of uniform random numbers

# linear congruential generators

# setting the parameters

m <- 256 # small example 
a <- 137
const <- 1
d <- 2

m <- 2^29 # generator Randu
a <- 2^16+3
const <- 0
d <- 3

m <- 2^29 # generator which was used in an Unix version
a <- 62605
const <- 113218009
d <- 3

# generating nrep elements of the sequence
nrep <- min(m,2^16)
x <- rep(1,nrep+d-1) 
for (i in (2:(nrep+d-1)))  x[i] <- (a*x[i-1]+const)%%m
hist(x/m,breaks=20)
n1 <- min(m,5000)
plot(x[1:n1],x[2:(n1+1)],cex=0.5)

#Visualisation for d>2
library(rggobi) # package with dynamical graphics 
n1 <- min(m,5000)
x3 <- cbind(x[1:n1],x[2:(n1+1)],x[3:(n1+2)])
ggobi(x3)
x5 <- cbind(x3,x[4:(n1+3)],x[5:(n1+4)])
ggobi(x5)


# Illustration of part 3 of the Theorem in the course.
m <- 64
a <- 11
x <- rep(1,m)
for (i in (2:m))  x[i] <- (a*x[i-1])%%m
x
a <- 13
x <- rep(1,m)
for (i in (2:m))  x[i] <- (a*x[i-1])%%m
x
(x[1:16]-1)/4
m2 <- m/4
const <- (a-1)/4
y <- rep(0,m2)
for (i in (2:m2))  y[i] <- (a*y[i-1]+const)%%m2
y

# Behaviour of low order bits
par(mfrow=c(4,1))
plot(x[1:64]%%2,type="l")
plot(x[1:64]%%4,type="l")
plot(x[1:64]%%8,type="l")
plot(x[1:64]%%16,type="l")
## x%%2^k is again a linear congruential generator with a'= a%%2^k


# Analysis of the lattice structure of the d-tupels
par(mfrow=c(1,1))


# setting up the matrix of d-tupels (shifted so that 0 is in the lattice)
L <- x[1:nrep] 
c1 <- 1
for (i in (1:(d-1))) {
  L <- cbind(L,x[(1+i):(nrep+i)]-c1*const)
  c1 <- a*c1+1
}


# generating system for the dual matrix
# (columns of the matrix Gd without the factor 1/m)
Gd <- diag(d) 
Gd[1,1] <- m
for (i in (2:d)) Gd[1,i] <- - ((a^(i-1))%%m)
Gd


# Visualisation for d=2: parallel equispaced lines which contain the
# lattice points 
plot(L[,1],L[,2],xlab="x[i]",ylab="x[i+1]",cex=0.2) # Plot of the pairs
abline(v=(10:20),col=2) # first set of lines, parallel to y axis
a1 <- m/Gd[2,i] 
b1 <- -Gd[1,i]/Gd[2,i]
k <- round(0.5*(Gd[1,i]+Gd[2,i]))-5 # selecting suitable intercepts
for (j in (k:(k+10))) abline(j*a1,b1,col=3) # second set of lines

#Reduction: project a vector onto another one and take a lattice point nearby
s1 <- sum(Gd[,2]*Gd[,2])
s2 <- sum(Gd[,1]*Gd[,2])
s2/s1
Gd[,1] <- Gd[,1]-round(s2/s1)*Gd[,2] 
Gd
a1 <- m/Gd[2,1] # set of lines after reduction
b1 <- -Gd[1,1]/Gd[2,1]
k <- round(0.5*(Gd[1,1]+Gd[2,1]))-5 
for (j in (k:(k+10))) abline(j*a1,b1,col=4) 
# exchanging the role of the two vectors 
s1 <- sum(Gd[,1]*Gd[,1])
s2 <- sum(Gd[,1]*Gd[,2])
s2/s1
Gd[,2] <- Gd[,2]-round(s2/s1)*Gd[,1]
Gd
a1 <- m/Gd[2,2] # set of lines after the reduction
b1 <- -Gd[1,2]/Gd[2,2]
k <- round(0.5*(Gd[1,2]+Gd[2,2]))-5 
for (j in (k:(k+10))) abline(j*a1,b1,col=5) 

  

# Simplifying a generating system 
Gd <- f.gitter(Gd) 
Gd
apply(Gd,2,function(x) sqrt(sum(x*x))) # lengths of the generating vectors

f.gitter <- function(L)
{
  ## Purpose: Modifying a generating system of a lattice such that the
  ##          vectors are shorter 
  ## -------------------------------------------------------------------------
  ## Arguments: L= matrix with columns = vectors of the generating system
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date:  1 Dec 2003, 17:50
  d <- nrow(L)
  change <- T 
  while (change) {
    change <- F
    for (j in (1:d)){
      s1 <- sum(L[,j]*L[,j])
      for (i in ((1:d)[-j])){
        s2 <- sum(L[,i]*L[,j])
        if (2*abs(s2)>s1){
          change <- T
          L[,i] <- L[,i]-round(s2/s1)*L[,j]
        }}}}
  L
}

# Studying the quality of generators for m=256 and arbitrary a.
m <- 256 
result <- 0
for (a in seq(1,m-1,by=2)) {
  Gd <- rbind(c(m,-a),c(0,1)) 
  Gd.neu <- f.gitter(Gd)
  while (any(Gd.neu!=Gd)) {
    Gd <- Gd.neu
    Gd.neu <- f.gitter(Gd.neu)
  }
  result <- c(result,min(Gd[1,1]^2+Gd[2,1]^2,Gd[1,2]^2+Gd[2,2]^2))
}
result <- cbind(seq(1,m-1,by=2),round(1/sqrt(result[-1]),3))
result


#Testing generators


#Test for overlapping trippels of 3 bits
n <- 2^16
x <- runif(n) # Testing the generator in R
x <- x[1:n]/m # Testing our own random numbers
x <- floor(8*x) #Bits 1,2,3 
x <- floor(16*x)%%8 #Bits 2,3,4 
x <- floor(32*x)%%8 #Bits 3,4,5 
x <- floor(1024*x)%%8 #Bits 8,9,10 
x <- c(x,x[1:2]) # periodic extension
y <- x[1:n]+8*x[2:(n+1)] # values of overlapping pairs
obs2 <- table(factor(y,0:63)) # counting the frequencies
q2 <- sum((obs2-2^10)^2/2^10) # Chisquare statistics
y <- y+64*x[3:(n+2)] # values for overlapping trippels 
obs3 <- table(factor(y,0:511)) # counting the frequencies
q3 <- sum((obs3-2^7)^2/2^7) # Chisquare statistics
1-pchisq(q3-q2,512-64) # correct chisquare test in case of overlap

#Overlapping pairs of 10 bits. We generate 2^21+1 values, in groups of
#length 4096. For each such group we count how many of the 2^20 pairs
#do not occur. r is the last random number in a such group which has
#to be carried over to the next group.
empty <- rep(TRUE,2^20)
r <- floor(1024*runif(1))
for (i in (1:512)){
  x <- c(r,floor(1024*runif(4096)))
  r <- x[4097]
  y <- 1+x[-4097]+1024*x[-1]
  empty[y] <- FALSE
}
(sum(empty)-141909)/290.26 #Standardisation following Marsaglia

#Overlapping 10-tupels of 2 bits. We generate 2^21+9 values, in groups of
#length 4096. For each such group we count how many of the 2^20 10-tupels
#do not occur. r contains the last 9 random number in a such group which has
#to be carried over to the next group.
empty <- rep(TRUE,2^20)
r <- floor(4*runif(9))
for (i in (1:512)){
  x <- c(r,floor(4*runif(4096)))
  r <- x[4097:4105]
  y <- 1+x[1:4096]+4*x[2:4097]+16*x[3:4098]+64*x[4:4099]+256*x[5:4100]+
       1024*x[6:4101]+4096*x[7:4102]+16384*x[8:4103]+65536*x[9:4104]+
       262144*x[10:4105]
  empty[y] <- FALSE
}
(sum(empty)-141911)/339 #Standardisation following Marsaglia

#Same thing for a linear congruential generator
m <- 2^29 # generator Randu
a <- 2^16+3
const <- 0
x0 <- 6899342 #starting value
x <- rep(x0,4105)
for (i in (2:9)) x[i] <- (a*x[i-1]+const)%%m
empty <- rep(TRUE,2^20)
for (j in (1:512)){
  for (i in (10:4105)) x[i] <- (a*x[i-1]+const)%%m
  z <- floor(4*x/m)
  x[1:9] <- x[4097:4105]
  y <- 1+z[1:4096]+4*z[2:4097]+16*z[3:4098]+64*z[4:4099]+256*z[5:4100]+
       1024*z[6:4101]+4096*z[7:4102]+16384*z[8:4103]+65536*z[9:4104]+
       262144*z[10:4105]
  empty[y] <- FALSE
}
(sum(empty)-141911)/339 #Standardisation following Marsaglia

# Information on random number generators in R
help(.Random.seed)


# Figures for a book

ps.latex(file="zufallszahl.eps",height=5)
par(mfrow=c(1,3),pty="s")
const <- 1
x <- rep(1,m+1)
a <- 165 #gute Wahl
for (i in (2:(m+1)))  x[i] <- (a*x[i-1]+const)%%m
L <- x[1:m] 
c1 <- 1
for (i in (1:(d-1))) {
  L <- cbind(L,x[(1+i):(nrep+i)]-c1*const)
  c1 <- a*c1+1
}
plot(L[,1]/m,L[,2]/m,xlab="u[i]",ylab="u[i+1]",xlim=c(0,1),ylim=c(0,1),pty="s")
a <- 125 # mttelgute Wahl
for (i in (2:(m+1)))  x[i] <- (a*x[i-1]+const)%%m
L <- x[1:m] 
c1 <- 1
for (i in (1:(d-1))) {
  L <- cbind(L,x[(1+i):(nrep+i)]-c1*const)
  c1 <- a*c1+1
}
plot(L[,1]/m,L[,2]/m,xlab="u[i]",ylab="u[i+1]",xlim=c(0,1),ylim=c(0,1),pty="s")
a <- 85  # schlechte Wahl
for (i in (2:(m+1)))  x[i] <- (a*x[i-1]+const)%%m
L <- x[1:m] 
c1 <- 1
for (i in (1:(d-1))) {
  L <- cbind(L,x[(1+i):(nrep+i)]-c1*const)
  c1 <- a*c1+1
}
plot(L[,1]/m,L[,2]/m,xlab="u[i]",ylab="u[i+1]",xlim=c(0,1),ylim=c(0,1),pty="s")
ps.end()