R Under development (unstable) (2024-10-17 r87242 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(randtoolbox)
Loading required package: rngWELL
This is randtoolbox. For an overview, type 'help("randtoolbox")'.
> 
> #see e.g. https://en.wikipedia.org/wiki/Linear_congruential_generator
> 
>   RNGkind()
[1] "Mersenne-Twister" "Inversion"        "Rejection"       
>   
>   #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
 mod mult incr 
TRUE TRUE TRUE 
>   x1 <- runif(5)
>   setSeed(12345)
>   x2 <- congruRand(5, dim=1, mod=2^31-1, mult=16807, incr=0)
>   print(cbind(x1, x2))
             x1         x2
[1,] 0.09661653 0.09661653
[2,] 0.83399463 0.83399463
[3,] 0.94770250 0.94770250
[4,] 0.03587859 0.03587859
[5,] 0.01154585 0.01154585
>   print(sum(abs(x1 - x2)))
[1] 0
>   
>   RNGkind()
[1] "user-supplied" "Inversion"     "Rejection"    
>   
>   # 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
 mod mult incr 
TRUE TRUE TRUE 
>   x1 <- runif(5)
>   setSeed(1)
>   x2 <- congruRand(5, dim=1, mod=4294967296, mult=1664525, incr=1013904223)
>   print(cbind(x1, x2))
             x1         x2
[1,] 0.23645553 0.23645553
[2,] 0.36927067 0.36927067
[3,] 0.50424203 0.50424203
[4,] 0.70488326 0.70488326
[5,] 0.05054363 0.05054363
>   print(sum(abs(x1 - x2)))
[1] 0
>   
> 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
+ }
> 
> proc.time()
   user  system elapsed 
   0.20    0.12    0.26