Uniform pseudo-random number generators

Arnab Rakshit

Sudipta Sarkar

2024-04-05

Introductory section

Introduction

Crude idea

“Lest men suspect your tale untrue, Keep probability in view.” — JOHN GAY (1727)

People who think about randomness almost invariably get into philosophical discussions about what the word random means.

What is RANDOM?

In a sense, there is no such thing as a random number; for example, is 2 a random number? Rather, we speak of a sequence of independent random numbers with a specified distribution, and this means loosely that each number was obtained merely by chance, having nothing to do with other numbers of the sequence, and that each number has a specified probability of falling in any given range of values.
There are several ways to formulate decent abstract definitions of randomness, but for the moment, let us content ourselves with an intuitive understanding of the concept.

History

Many years ago, people who needed random numbers in their scientific work would draw balls out of a “well-stirred urn,” or they would roll dice or deal out cards.
A table of over 40,000 random digits, “taken at random from census reports,” was published in 1927 by L. H. C. Tippett. Since then, a number of devices have been built to generate random numbers mechanically.
The first such machine was used in 1939 by M. G. Kendall and B. Babington-Smith to produce a table of 100,000 random digits. The Ferranti Mark I computer, first installed in 1951, had a built-in instruction that put 20 random bits into the accumulator using a resistance noise generator; this feature had been recommended by A. M. Turing.
In 1955, the RAND Corporation published a widely used table of a million random digits obtained with the help of another special device. A famous random-number machine called ERNIE has been used for many years to pick the winning numbers in the British Premium Savings Bonds lottery.

Motivation & purpose

Numbers that are “chosen at random” are useful in many different kinds of applications.
For example:
- Simulation: When a computer is being used to simulate natural phenomena, random numbers are required to make things realistic. Simulation covers many fields, from the study of nuclear physics (where particles are subject to random collisions) to operations research (where people come into, say, an airport at random intervals).
- Sampling: It is often impractical to examine all possible cases, but a random sample will provide insight into what constitutes “typical” behavior. Numerical analysis. Ingenious techniques for solving complicated numerical problems have been devised using random numbers. Several books have been written on this subject.
- Computer programming: Random values make a good source of data for testing the effectiveness of computer algorithms. More importantly, they are crucial to the operation of randomized algorithms, which are often far superior to their deterministic counterparts.
- Decision making: There are reports that many executives make their decisions by flipping a coin or by throwing darts, etc. It is also rumored that some college professors prepare their grades on such a basis. Sometimes it is important to make a completely “unbiased” decision. Randomness is also an essential part of optimal strategies in the theory of matrix games.
- Cryptography: A source of unbiased bits is crucial for many types of secure communications, when data needs to be concealed.

Goal

Uniform(0,1) random numbers are the key to random variate generation in simulation. Then necessary transformations can be done to get other RVs. So here we specifically reduce our goal to Give an algorithm that produces a sequence of pseudo-random numbers (PRNs) R1, R2, … that “appear” to be iid Unif(0,1).
Some desired properties of the proposed algorithms should be:
- Output appears to be iid Unif(0,1)
- Correctness
- Works fast
- Ability to reproduce any sequence it generates

Development of Pseudo-random number generators

Some naive generators(which we wouldn’t use)

Random devices

Though these holds randomness properties, but stroing Unif(0,1) sequence being mot so easy, it’s tough to repeat the experiment.
Examples:
- Flipping a coin
- Particle count by Geiger counter
- Least significant digits of atomic clock

Random number tables

This refers to the list of digits supplied in tables. These method is cumbersome and very much time-consuming. Besides, whenever the numbers are tabled, the method no longer remains random. Hence these methods should not be used, in general.

Middle-square method

History

The method was invented by John von Neumann, and was described by him at a conference in 1949.
In the 1949 talk, Von Neumann quipped that “Anyone who considers arithmetical methods of producing random digits is, of course, in a state of sin.” What he meant, he elaborated, was that there were no true “random numbers”, just means to produce them, and “a strict arithmetic procedure”, like the middle-square method, “is not such a method”. Nevertheless, he found these methods hundreds of times faster than reading “truly” random numbers off punch cards, which had practical importance for his ENIAC work. He found the “destruction” of middle-square sequences to be a factor in their favor, because it could be easily detected: “one always fears the appearance of undetected short cycles”. Nicholas Metropolis reported sequences of 750,000 digits before “destruction” by means of using 38-bit numbers with the “middle-square” method.

The method

To generate a sequence of n-digit pseudo-random numbers, an n-digit starting value is created and squared, producing a 2n-digit number.
If the result has fewer than 2n digits, leading zeroes are added to compensate.
The middle n digits of the result would be the next number in the sequence and returned as the result.
This process is then repeated to generate more numbers.

Drawback

n must be even: The value of n must be even in order for the method to work – if the value of n is odd, then there will not necessarily be a uniquely defined “middle n-digits” to select from. Consider the following: If a 3-digit number is squared, it can yield a 6-digit number (e.g. 5402 = 291600). If there were to be middle 3 digits, that would leave 6 − 3 = 3 digits to be distributed to the left and right of the middle. It is impossible to evenly distribute these digits equally on both sides of the middle number, and therefore there are no “middle digits”. It is acceptable to pad the seeds with zeros to the left in order to create an even valued n-digit number (e.g. 540 → 0540).
Repetition: For a generator of n-digit numbers, the period can be no longer than 8n. If the middle n digits are all zeroes, the generator then outputs zeroes forever. If the first half of a number in the sequence is zeroes, the subsequent numbers will be decreasing to zero. While these runs of zero are easy to detect, they occur too frequently for this method to be of practical use. The middle-squared method can also get stuck on a number other than zero. For n = 4, this occurs with the values 0100, 2500, 3792, and 7600. Other seed values form very short repeating cycles, e.g., 0540 → 2916 → 5030 → 3009. These phenomena are even more obvious when n = 2, as none of the 100 possible seeds generates more than 14 iterations without reverting to 0, 10, 50, 60, or a 24 ↔︎ 57 loop.

Illustration

Here we fix n = 4 and carry out the middle-square algorithm, generate 1000 random numbers and study the behavioral patterns of the algorithm.
Here we fix the initial entry as 6969, then generate random numbers from 0001 to 9999, both bounds being inclusive, then represent them as random numbers from uniform[0,1] distribution, by operating necessary normalization.

Code

# Function to generate random numbers using Middle Square Method
middle_square <- function(seed, digits, n) {
  result <- numeric(n)
  for (i in 1:n) {
    # Square the seed
    squared <- seed^2
    # Extract middle digits
    middle <- as.numeric(substr(as.character(squared), start = (nchar(squared) - digits) / 2 + 1, stop = (nchar(squared) + digits) / 2))
    # Update seed
    seed <- middle 
    # Normalize to [0, 1]
    result[i] <- middle / (10^digits)
  }
  return(result)
}


# Set parameters
initial_seed <- 6969
num_digits <- 4
num_random_numbers <- 1000

# Generate random numbers
set.seed(initial_seed) # Set seed for reproducibility
random_numbers <- middle_square(initial_seed, num_digits, num_random_numbers)

# Print generated random numbers
print(random_numbers)

   [1] 0.5669 0.1375 0.8906 0.3168 0.0362 0.3104 0.6348 0.2971 0.8268 0.3598
  [11] 0.9456 0.4159 0.2972 0.8327 0.3389 0.4853 0.5516 0.4262 0.1646 0.7093
  [21] 0.3106 0.6472 0.8867 0.6236 0.8876 0.7833 0.3558 0.6593 0.4676 0.8649
  [31] 0.8052 0.8347 0.6724 0.2121 0.4986 0.8601 0.9772 0.4919 0.1965 0.8612
  [41] 0.1665 0.7722 0.6292 0.5892 0.7156 0.2083 0.3388 0.4785 0.8962 0.3174
  [51] 0.0742 0.5056 0.5631 0.7081 0.1405 0.9740 0.8676 0.2729 0.4474 0.0166
  [61] 0.2755 0.5900 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
  [71] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
  [81] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
  [91] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [101] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [111] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [121] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [131] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [141] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [151] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [161] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [171] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [181] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [191] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [201] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [211] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [221] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [231] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [241] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [251] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [261] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [271] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [281] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [291] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [301] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [311] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [321] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [331] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [341] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [351] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [361] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [371] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [381] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [391] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [401] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [411] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [421] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [431] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [441] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [451] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [461] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [471] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [481] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [491] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [501] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [511] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [521] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [531] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [541] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [551] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [561] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [571] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [581] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [591] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [601] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [611] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [621] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [631] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [641] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [651] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [661] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [671] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [681] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [691] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [701] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [711] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [721] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [731] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [741] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [751] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [761] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [771] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [781] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [791] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [801] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [811] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [821] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [831] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [841] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [851] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [861] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [871] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [881] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [891] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [901] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [911] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [921] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [931] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [941] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [951] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [961] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [971] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100
 [981] 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100
 [991] 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100 0.2100 0.4100 0.8100 0.6100

Here it can be observed that starting from the 63^rd random number, the algorithm becomes stuck in the cycle, 2100 → 4100 → 8100 → 6100
For convenience, the frequency distribution of the generated random numbers is shown here. Just look at the frequencies of 0.2100, 0.4100, 0.6100 and 0.8100!

Code

table(middle_square(6969,4,1000))


0.0166 0.0362 0.0742 0.1375 0.1405 0.1646 0.1665 0.1965 0.2083   0.21 0.2121 
     1      1      1      1      1      1      1      1      1    234      1 
0.2729 0.2755 0.2971 0.2972 0.3104 0.3106 0.3168 0.3174 0.3388 0.3389 0.3558 
     1      1      1      1      1      1      1      1      1      1      1 
0.3598   0.41 0.4159 0.4262 0.4474 0.4676 0.4785 0.4853 0.4919 0.4986 0.5056 
     1    234      1      1      1      1      1      1      1      1      1 
0.5516 0.5631 0.5669 0.5892   0.59   0.61 0.6236 0.6292 0.6348 0.6472 0.6593 
     1      1      1      1      1    235      1      1      1      1      1 
0.6724 0.7081 0.7093 0.7156 0.7722 0.7833 0.8052   0.81 0.8268 0.8327 0.8347 
     1      1      1      1      1      1      1    235      1      1      1 
0.8601 0.8612 0.8649 0.8676 0.8867 0.8876 0.8906 0.8962 0.9456  0.974 0.9772 
     1      1      1      1      1      1      1      1      1      1      1

ACF Plots

Here we have plotted the ACF plots, taking the initial entries as 1346, 2795, 6969, 9940. It should be mentioned that these numbers are not chosen randomly though.

Code

# ACF plot
random_numbers <- middle_square(initial_seed, num_digits, num_random_numbers)
  par(mfrow=c(2,2))
  acf(middle_square(1346,4,1000))
  acf(middle_square(2795,4,1000))
  acf(middle_square(6969,4,1000))
  acf(middle_square(9940,4,1000))

Histograms

Code

# Histogram plot
par(mfrow=c(2,2))
hist(middle_square(1346,4,1000))
hist(middle_square(2795,4,1000))
hist(middle_square(6969,4,1000))
hist(middle_square(9940,4,1000))