tests/test-runifInterface-LCG.R

In randtoolbox: Toolbox for Pseudo and Quasi Random Number Generation and Random Generator Tests

library(randtoolbox)

#see e.g. https://en.wikipedia.org/wiki/Linear_congruential_generator

  RNGkind()
  
  #Park Miller congruential generator : mod = 2147483647= 2^31-1, mult=16807, incr=0
  LCG.par <- c("2147483647", "16807", "0")
  set.generator(name="congruRand", mod=LCG.par[1], mult=LCG.par[2], incr=LCG.par[3], seed=12345)
  LCG <- get.description()
  LCG$parameters == LCG.par
  x1 <- runif(5)
  setSeed(12345)
  x2 <- congruRand(5, dim=1, mod=2^31-1, mult=16807, incr=0)
  print(cbind(x1, x2))
  print(sum(abs(x1 - x2)))
  
  RNGkind()
  
  # the Knuth Lewis RNG : mod=4294967296 == 2^32, mult=1664525, incr=1013904223
  LCG.par <- c("4294967296", "1664525", "1013904223")
  set.generator(name="congruRand", mod=LCG.par[1], mult=LCG.par[2], incr=LCG.par[3], seed=1)
  LCG <- get.description()
  LCG$parameters == LCG.par
  x1 <- runif(5)
  setSeed(1)
  x2 <- congruRand(5, dim=1, mod=4294967296, mult=1664525, incr=1013904223)
  print(cbind(x1, x2))
  print(sum(abs(x1 - x2)))
  
if(.Platform$OS.type != "windows")
{
  # the POSIX rand48 : 281474976710656 == 2^48
  LCG.par <- c("281474976710656", "25214903917", "11")
  set.generator(name="congruRand", mod=LCG.par[1], mult=LCG.par[2], incr=LCG.par[3], seed=1)
  LCG <- get.description()
  LCG$parameters == LCG.par
  x1 <- runif(5)
  setSeed(1)
  x2 <- congruRand(5, dim=1, mod=281474976710656, mult=25214903917, incr=11)
  print(cbind(x1, x2))
  print(sum(abs(x1 - x2)))
}
  
if(FALSE) #congruRand() does not handle 2^64 correctly but set.generator() does 
{
  
  # the MMIX RNG by Donald Knuth => produce two different result after the second term
  # 18446744073709551616 == 2^64
  LCG.par <- c("18446744073709551616", "1442695040888963407", "1013904223")
  set.generator(name="congruRand", mod=LCG.par[1], mult=LCG.par[2], incr=LCG.par[3], seed=1)
  LCG <- get.description()
  LCG$parameters == LCG.par
  x1 <- runif(5)
  setSeed(1)
  x2 <- congruRand(5, dim=1, mod=18446744073709551616, mult=1442695040888963407, incr=1013904223)
  print(cbind(x1, x2))
  #only first value is correct
  (1442695040888963407 * 1 + 1013904223) / 2^64
  
  #Haynes RNG
  LCG.par <- c("18446744073709551616", "636412233846793005", "1")
  set.generator(name="congruRand", mod=LCG.par[1], mult=LCG.par[2], incr=LCG.par[3], seed=1)
  LCG <- get.description()
  LCG$parameters == LCG.par
  x1 <- runif(5)
  setSeed(1)
  x2 <- congruRand(5, dim=1, mod=18446744073709551616, mult=636412233846793005, incr=1, echo = TRUE)
  print(cbind(x1, x2))
  
  #only first value is correct
  (636412233846793005 * 1 + 1) / 2^64
}