m <- 2^48
a <- 25214903917
c <- 11

#m <- 10000
#a <- 17
#c <- 1

X0 <- 0

n <- 10000
Xi <- X0
X <- numeric(n)
for (i in 1:n) {
  Xi <- (a*Xi + c) %% m
  # cat("X[", i, "] =", Xi, "\n")
  X[i] <- Xi
}
#hist(X, breaks = seq(-0.5, m-0.5, by = 1))

U <- (X + 1) / (m + 1)
hist(U)
plot(U[-n], U[-1])

U <- runif(n)
plot(U[-n], U[-1])