How can a deterministic computer produce genuinely random numbers? And does it even need to? Yes, computers need to have access to random numbers. They are used to encrypt information, deal cards in your game of virtual Solitaire, simulate unknown variables like in weather prediction and airplane scheduling, carry out a computer simulation of physical phenomena, notably simulation of neutron transport in nuclear fission, and so much more. A little bit of randomness helps computers to create this kind of beautiful graphics:

Although random numbers are required in many applications, their generation is often overlooked. We program computers, that is their actions are determined by a prescribed set of instructions. As computers are deterministic, they are not capable of producing truly random numbers. How can it possibly produce a random number?

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.

Computers can generate sequences of truly random numbers, but they depend on some type of external input from a process which we assume to be random. They inherit the randomness of the nature using an extra piece of hardware. For such devices, a physical source of randomness (entropy) can be signals, such as atmospheric noise, thermal noise, radiation, or other external electromagnetic and quantum phenomena. This can also be a good source of entropy for random number generation: the distribution of pepperoni on Aldi pizza.

The picture below illustrates an example of a physical random number generator which exploits an elementary quantum optics process. Photons – light particles – are sent one by one onto a semi-transparent mirror and detected. The exclusive events (reflection-transmission) are associated with « 0 » – « 1 » bit values. This results in producing a sequence of random bits called a bit-stream which is understood as a number by computers.

Hardware random number generators have major disadvantages.

The bit-stream from such devices is prone to be biased. For example, a semi-transparent mirror photon detector in the above picture generates bit streams that consist of mostly “0” with the occasional “1”. Thus, they must be carefully designed to minimize the bias.

The second drawback is that they’re not very fast. It takes a while to harvest all that external data from a natural phenomenon and fill a bit-stream with enough random bits.

Pseudorandom Number Generators

Fortunately, for most applications, pseudorandom numbers are sufficient which has many of the statistical properties associated with truly random numbers. Why are they called Pseudorandom numbers? The reason is these

numbers are generated by systematic methods, that is they follow some human-defined patterns, an algorithm.

To show you what I mean, here is John von Neumann’s Middle-Square Algorithm (don’t confuse it with the Middle-Out algorithm from Silicon Walley ) whose idea was to take the square of the previous random number and to extract the middle digits. If you somehow determine the previous number, you will be able to predict the entire sequence of these “random” numbers. In addition, all sequences of the “middle square” method eventually repeat themselves.

Any one who considers arithmetical methods of producing random digits is, of course, in a state of sin John Von Neumann

Linear Congruential Generator

This is definitely not the only way to create a random-looking sequence of numbers. The simple, yet old and best-known method is the linear congruential generator (LCG). LCG generates a sequence of pseudorandom numbers by successively using the recursively defined function

\begin{equation*}

X_{x+1}=(aX_{n}+c) \mod m\end{equation*}

with the modulus m, the multiplier a, the increment c, and the seed, the sequence’s initial value. Let’s Find the sequence of pseudorandom numbers generated by the linear congruential method with modulus m = 9, multiplier a = 7, increment c = 4, and seed \(x_{0}=3\). We obtain that

\(x_{1}= 7*x_{0} + 4 \mod 9 = 7 * 3 + 4 \mod 9 = 25 \mod 9 = 7\),

\(x_{2}=7*x_{1} + 4 \mod 9 = 7 * 7 + 4\mod 9 = 53 \mod 9 = 8\),

\(x_{3}=7*x_{2} + 4 \mod 9 = 7 * 8 + 4\mod 9 = 60 \mod 9 = 6\),

\(x_{4}= 7*x_{3}+ 4 \mod 9 = 7 * 6 + 4 \mod 9 = 46 \mod 9 = 1\),

\(x_{5}= 7*x_{4}+ 4 \mod 9 = 7 * 1 + 4 \mod 9 = 11 \mod 9 = 2\),

\(x_{6}=7*x_{5}+ 4\mod 9 = 7 * 2 + 4 \mod 9 = 18 \mod 9 = 0\),

\(x_{7}= 7*x_{6} + 4 \mod 9 = 7 * 0 + 4 \mod 9 = 4 \mod 9 = 4\),

\(x_{8}=7*x_{7} + 4 \mod 9 = 7 * 4 + 4 \mod 9 = 32 \mod 9 = 5\),

\(x_{9}= 7*x_{8}+ 4 \mod 9 = 7 * 5 + 4 \mod 9 = 39 \mod 9 = 3\).

Because \(x_{9} = x_{0}\) and because each term depends only on the previous term, we see that the sequence

3, 7, 8, 6, 1, 2, 0, 4, 5, 3, 7, 8, 6, 1, 2, 0, 4, 5, 3, . . .

is generated. This sequence contains nine different numbers before repeating. In other words, it has a period equal to n=9. If we choose a big value for m, like \(2^{32}\) kind of big, the results may be sufficiently complex to make the pattern difficult to identify, rather than the middle-square algorithm which has a very short period. That is why LCG is more preferred pseudorandom number generator.

Here is an implementation of an LCG in Python (Wikipedia):

def lcg(modulus, a, c, seed): while True: seed = (a * seed + c) % modulus yield seed 1 2 3 4 def lcg ( modulus , a , c , seed ) : while True : seed = ( a * seed + c ) % modulus yield seed

Despite their downfalls, the huge benefit of pseudorandom generators over hardware random number ones is you just need to know the seed to generate the identical sequence which is useful for testing and fixing software. Also, there is no need for any kind of hardware.

Pseudorandom Number Sampling

Pseudorandom numbers generated by LCG has the uniform distribution meaning that every number has equal probability of being returned. Suppose you have a fair coin, when you flip it, it has a 50% chance of landing heads up and a 50% chance of landing tails up, right? Similarly, if you roll a die, it comes to rest showing on its upper surface a random integer from one to six, each value being equally likely. Outcomes of a fair coin and die follow the uniform distribution. However, sometimes we need to make a die or a coin loaded so that they will land with a specific side facing upwards more or less often than a fair one would. In other words, we want to generate numbers that come from other distributions.

There exists a couple of methods to generate a random number according to other distribution. These methods involve manipulating a uniform random number in some way.

Inverse Transform Sampling

One method is called Inverse Transform Sampling. The method is very simple, so I’ll try to describe it in simple words. First, we take cumulative distribution function \(F_{X}\) of some distribution that we want to sample from. The function accepts as input some value x and tells you what is the probability of obtaining \(X \leq x\). So

\begin{equation*}

F_{X}(x)=Pr(X \leq x)=p \end{equation*}

the inverse of such function \(F_{X}^{-1}\) would take p as input and return x. Notice that p’s are uniformly distributed – this could be used for sampling from any \(F_{X}\) if you know \(F_{X}^{-1}\). Since it is easy for computers to generate pseudorandom numbers from the uniform distribution, we ask them to give us pseudorandom numbers u ~ U(0,1) between 0 and 1. Then we just pass numbers u through \(F_{X}^{-1}\) to get x‘s

\begin{equation*}

F_{X}^{-1}(u)=x\end{equation*}

To visualize this look at CDF below, generally, we think of distributions in terms of looking at y-axis for probabilities of values from the x-axis. With this sampling method, we do the opposite and start with “probabilities” and use them to pick the values that are related to them on the x-axis.

Unfortunately, this is not always possible since not every function has its inverse, e.g. we cannot use this method with bivariate distributions. Therefore, we need other methods.

Rejection Sampling

The second method is called the Acceptance-Rejection method or just the Rejection Method. The idea is based on “Just throw the bad crap away“. This involves choosing an x and y value and testing whether the function of x is greater than the y value. If it is, the x value is accepted. Otherwise, the x value is rejected and the algorithm tries again.

Here is a great explanation of rejection sampling by Agustinus Kristiadi who describes the method in detail.

Suppose that our task is to produce Gaussian Random values. Of course, there is a bunch of efficient algorithms to generate these values. However, to illustrate rejection sampling, we will use the Marsaglia polar method. One begins by generating a point (x,y) uniformly at random in the unit disc. We then let

\( r= x^2 + y^2 \), \( \mu = \sqrt{\frac{-2 \log 2}{r}} \), \(Z_{1} = \mu x\) and \( Z_{2} = \mu y\)

Then we obtain independent Gaussian random variables \(Z_{1} \) and \( Z_{2} = \). Here is Python program based on this method:

def polar(): while True: x = random.random() * 2 - 1 y = random.random() * 2 - 1 r = x*x + y*y if r < 1: mu = math.sqrt(-2.0 * math.log(r) / r) return mu*x, mu*y 1 2 3 4 5 6 7 8 9 10 def polar ( ) : while True : x = random . random ( ) * 2 - 1 y = random . random ( ) * 2 - 1 r = x * x + y * y if r < 1 : mu = math . sqrt ( - 2.0 * math . log ( r ) / r ) return mu * x , mu * y



As you see, all the methods mentioned above amazingly allows us to transform the sequence of uniform 0, 1 random or pseudorandom numbers into a sequence with a different distribution. Even with a deterministic computer, the truly random and pseudorandom techniques we have explored give us numbers with the uniform 0, 1 distribution, which in turn, gives us sequences with other distributions, like wind speeds or dice rolls.

Mathematical Problems

Now, I want to show you, guys, how to program a computer to solve some math problems for us. They are simple, yet interesting. Also, they may give you an idea of how scientists simulate real-world problems. Again, I will use Python, but you can choose any programming language that you are familiar with. Python uses the Mersenne Twister algorithm for generating pseudorandom numbers.

Relatively Prime Numbers

Let \( p_{n} \) be the probability that two numbers, sampled independently and uniformly from {1,2, . . . ,n} are relatively prime (that is, their greatest common divisor is one). What can we say about \( p_{n} \) as n \to \infty \)?

If we use direct enumeration (using loops) to calculate \( p_{n} \) for various values of n, as n approaches 100,000, the time for the computer to complete the calculation becomes excessive. This motivates us to find another attack. The approach we take is to sample pairs of integers in {1,2, . . . ,n} at random and keep track of how often we find that they are relatively prime.We now present a simple computer program to estimate \( p_{n} \).

import random n = int(input("Enter n (max el't of the set): ")) reps = int(input("Enter the number of pairs to sample: ")) count = 0 for i in range(reps): a = random.randint(1, n) b = random.randint(1, n) if (gcd(a,b) == 1): count += 1 print(count/reps) 1 2 3 4 5 6 7 8 9 10 11 12 13 import random n = int ( input ( "Enter n (max el't of the set): " ) ) reps = int ( input ( "Enter the number of pairs to sample: " ) ) count = 0 for i in range ( reps ) : a = random . randint ( 1 , n ) b = random . randint ( 1 , n ) if ( gcd ( a , b ) == 1 ) : count += 1 print ( count / reps )

NOTE: The program just gives us the estimate \( p_{n} \). The value of \( p_{n} \) we get from our program is correct to some decimal places. To obtain greater accuracy, we need either to increase the number of repetitions or find a better method for calculating \( p_{n} \).

Random Walk

Maybe you have heard about the random walk before. We create a procedure that simulates a random walk on the integers. That is, our procedure returns the position of a particle whose initial position is 0. Each time the function is called, the particle moves one step left or right, each with probability 50%. The return value will be the new position of the particle. Then using additional library matplotlib, we plot down on the graph of this process in 1D space.