# CS 590M
# Peter J. Haas

# This file contains some R code to demonstrate the pitfalls 
# of linear congruential generators
#============================================================

# Code for the various random number generators

# These implementations are NOT efficient in general, they just illustrate the pitfalls

# generate vector of random numbers from Lewis-Learmonth LCG
LLgen = function(seed,num) {
  vec = c(seed)
  for (i in 2:num) {vec = append(vec,(16807*vec[i-1]) %% 2147483647)}
  vec/2147483647
}

# generate vector of random numbers for Box-Muller demo
# (not an LCG that you would use in practice, but it generates a clear plot)

BMgen = function(seed,num) {
  vec = c(seed)
  for (i in 2:num) {vec = append(vec,(7*vec[i-1]) %% 2^31)}
  vec/2^31
}

# generate vector of random numbers from RANDU

RANDUgen = function(seed,num) {
  vec = c(seed)
  a = 2^16 + 3
  for (i in 2:num) {vec = append(vec,(a*vec[i-1]) %% 2^31)}
  vec/2^31
}

# Set up Mersenne Twister

mtplot = function(x,y) plot(x,y, main = "(x1,x2) pairs from Mersenne Twister",
  xlab="x1",ylab="x2",xlim=c(0,1),ylim=c(0,1),pin=c(1.0,1.0))

mttest1 = function(jump) {
z = runif(jump)
u = runif(100)
x = u[seq(1,100,2)]
y = u[seq(2,100,2)]
par(pty="s")
mtplot(x,y)
}

mttest2 = function(jump) {
.Random.seed <<- c(403,10,1:624)
.Random.seed <<- as.integer(.Random.seed)
z = runif(jump)
u = runif(100)
x = u[seq(1,100,2)]
y = u[seq(2,100,2)]
par(pty="s")
mtplot(x,y)
}

#============================================================
#============================================================
#============================================================

# Pitfalls of choosing initial seed for LL generator

LLinit = function(numpts) {
  x <- LLgen(1,numpts) # Note: seeds 1 and 2 are separated by over 1.3 billion steps
  y <- LLgen(2,numpts)
  par(ask=TRUE)
  plot(x,y)
  hist(x+y)
}

LLinit(10000)

#============================================================
# Pitfalls of using Box-Muller method with an LCG

par(ask=F)
boxmuller = function(numpts) {
  unifs = BMgen(2^31-1,2*numpts)
  u1 = unifs[seq(1,2*numpts,2)]
  u2 = unifs[seq(2,2*numpts,2)]
  x = sqrt(-2*log(u1))*cos(2*pi*u2)
  y = sqrt(-2*log(u1))*sin(2*pi*u2)
  plot(x,y)
}

boxmuller(5000)

#============================================================
# Pitfalls of initializing Mersenne twister

mttest1(0) # plot 50 pairs from beginning of default stream
mttest1(1000) # jump by 1000 steps and plot next 50 pairs

mttest2(0) # same as above, but initial seed is (1,2,...,623)
mttest2(100)
mttest2(1000)
mttest2(10000)

#============================================================

# Pitfalls of RANDU LCG

# The following package supports rotatable plots
# (may need to install from CRAN repository if not present)
library(rgl)

randu2 = function(numpts) {
  unifs = RANDUgen(1,2*numpts)
  x = unifs[seq(1,2*numpts,2)]
  y = unifs[seq(2,2*numpts,2)]
  par(pty="s")
  plot(x,y)
}

randu2(500)

randu3 = function(numpts) {
  unifs = RANDUgen(1,3*numpts)
  x = unifs[seq(1,3*numpts,3)]
  y = unifs[seq(2,3*numpts,3)]
  z = unifs[seq(3,3*numpts,3)]
  plot3d(x,y,z)
}

randu3(5000)