./articles/ Understanding Random Number Generation

Despite its everyday use, the mechanics of Random Number Generation (RNG) are often misunderstood, particularly in determining how "randomness" is generated and ensuring its reliability across different use cases. In this article, we will explore two common random number generators (RNGs): the Linear Congruential Generator (LCG) and the Permuted Congruential Generator (PCG). But first let's understand what "random" means.

If your question is: are there truly random algorithms? I want to tell you bad news. Algorithms are not truly random because they are inherently deterministic by design. True randomness originates from unpredictable, non-deterministic sources, while algorithms—being sets of defined instructions—cannot achieve this unpredictability.

A fundamental characteristic of an algorithm is that it consistently produces the same output given the same input. This reproducibility allows outcomes to be replicated and even predicted, setting them apart from true randomness. Furthermore, algorithms operate within a finite computational environment, constrained by systems like integers or floating-point numbers with defined limits. These restrictions inevitably lead to repetition in the output, forming what is known as the period of a pseudo-random sequence.

To approximate true randomness, external entropy sources must be incorporated. Cryptographically Secure Pseudo-Random Number Generators (CSPRNGs) employ this strategy by using entropy from unpredictable external sources, such as hardware random number generators (HRNGs). In contrast, Pseudo-Random Number Generators (PRNGs) rely entirely on an initial input or seed, producing outputs that only appear random. Techniques like the LCG and the PCG exemplify these methods, which we will implement later.

True Random Number Generators (TRNGs), however, derive randomness from physical phenomena, such as atmospheric noise, radioactive decay, or thermal fluctuations. These sources are inherently unpredictable, yielding a higher degree of randomness compared to algorithmic approaches.

Therefore, our focus will remain on PRNGs and their deterministic, yet practical, approach to generating sequences that approximate randomness.

1. Linear Congruential Generator (LCG)

The Linear Congruential Generator (LCG) is a widely-used and simple pseudo-random number generator based on a linear recurrence relation:

\[ X_{n+1} = (a \cdot X_n + c) \mod m \]

Where:

• \(a\): The multiplier that determines how much the sequence progresses between consecutive numbers.
• \(c\): The increment that introduces variability, preventing the sequence from falling into zero cycles.
• \(m\): The modulus that sets the range of the generated numbers, often chosen as a power of 2 for computational efficiency.
• \(X_n\): The current state, which can be either the initial seed or the previously generated number.

Let's build a Linear Congruential Generator (LCG) from scratch in Python for simplicity. Our first task is to define a class called LCGPseudoRand, with constructor parameters for a, c, m, and seed. Here's the implementation:

    class LCGPseudoRand:
        def __init__(self, a=48271, c=0, m=2**31-1, seed=1):  # C++11 LCG parameters
            self.a = a
            self.c = c
            self.m = m
            self.x0 = seed
            # Compute the first value based on the seed
            self.prev_x = (a * self.x0 + c) % m
                

In this implementation, I've used the C++11 LCG parameters with a default seed value of 1. The constructor also initializes the first value of the sequence using the LCG formula:

    self.prev_x = (a * self.x0 + c) % m
                

For example, if we create the object lcg = LCGPseudoRand(), the first value is calculated as:

(48271 * 1 + 0) % (2**31-1) = 48271

The starting value of the sequence is 48271. Now, let's add a method to generate the pseudo-random numbers based on the same formula:

    def generate_number(self, min_range=None, max_range=None):
        self.prev_x = (self.a * self.prev_x + self.c) % self.m
        return self.prev_x
                

The method updates the state using the LCG formula and returns the next number in the sequence. Let's generate 5 pseudo-random numbers and print them:

    lcg = LCGPseudoRand()
    for _ in range(5):
        print(lcg.generate_number())
                

The output will look like this:

182605794, 1291394886, 1914720637, 2078669041, 407355683

We have successfully generated 5 pseudo-random numbers.

• The first element of the list is calculated as:

  (48271 * 48271 + 0) % 2147483647 = 182605794

• The second element is calculated as:

  (48271 * 182605794 + 0) % 2147483647 = 1291394886

However, if you re-run the same code without changing the seed, the algorithm will produce the exact same sequence again. This is because LCGs are deterministic. Their output is entirely determined by the initial seed and the parameters a, c, and m.

As we discussed earlier, a more unpredictable starting seed is required to avoid deterministic sequences. Atmospheric noise or radioactive decay might be difficult to integrate. A simple yet effective alternative is to combine the current time with the system's process ID. This method works well and ensures better randomness when initializing the generator's seed in the constructor.

    import os
    import time
    
    class LCGPseudoRand:
        def __init__(self, a=48271, c=0, m=2**31-1, seed=None):
            self.a = a
            self.c = c
            self.m = m
            if seed is None:
                self.x0 = int(os.getpid() + time.time())
            else:
                self.x0 = seed
            self.prev_x = (a * self.x0 + c) % m
                    

In this example, we've added functionality to automatically generate a unique seed if the user doesn't provide one. The line self.x0 = int(os.getpid() + time.time()) combines the system's process ID with the current time to ensure that each execution generates a different seed.

The os.getpid() function fetches the process ID of the current program, which is usually unique per running program. The time.time() function returns the current time in seconds since the epoch, which changes with each program execution. When you add these together and cast them to an integer, the result becomes a fairly unique seed that minimizes the likelihood of repeating values. Notice that each time the program runs, a new sequence is generated, thanks to our dynamic seed generation mechanism.

When generating pseudo-random numbers with a LCG, it is often necessary to scale the output so that it falls within a specified range. For example, you may need a random number between 10 and 50, or between two user-defined values, such as min_range and max_range. By scaling the generated number, you can ensure that it fits within the desired bounds.

The scaling formula is:

\[ \text{output} = \text{min_range} + \left( \frac{X_n}{m} \right) \times (\text{max_range} - \text{min_range}) \]

In the following code snippet, we modify the previously generated pseudo-random number using a scaling formula to fit within the user-defined range:

    def generate_number(self, min_range=None, max_range=None):
        self.prev_x = (self.a * self.prev_x + self.c) % self.m
    
        if min_range is not None and max_range is not None:
            return int((self.prev_x / self.m) * (max_range - min_range) + min_range)
        else:
            return self.prev_x
                    

The generate_number method first calculates the next number in the sequence using the standard LCG formula, and then scales the result to fit within the specified min_range and max_range.

Let’s consider an example where we generate random numbers between 10 and 50:

    lcg = LCGPseudoRand()
    for _ in range(5):
        print(lcg.generate_number(10, 50))
                    

Each time this program runs, the output numbers will be scaled to the specified range (10 to 50).

Initialization Differences Between Programming Languages

The parameters \(a\), \(c\), and \(m\) directly affect the quality of randomness and the period of the LCG. A poorly chosen set of parameters may result in sequences with clear patterns or shorter-than-expected periods. To achieve the maximum period of \(m\), the following conditions must be met:

• The increment \(c\) must be relatively prime to \(m\).
• The multiplier \(a - 1\) must be divisible by all prime factors of \(m\).
• If \(m\) is a multiple of 4, \(a\) and \(m\) should share the same remainder when divided by 4.

Different programming languages and libraries initialize LCGs using unique approaches to determine the initial seed \(X_0\), which impacts reproducibility and sequence variation:

long seed = System.currentTimeMillis();

#include <random>
            std::linear_congruential_engine<unsigned, 16807, 0, 2147483647> engine(12345);

#include <random>
    std::minstd_rand rng(42);

2. Permuted Congruential Generator (PCG)

The Permuted Congruential Generator (PCG) is a modern pseudo-random number generator that improves upon the classic Linear Congruential Generator (LCG). It combines a linear congruential generator with a permutation function to produce higher-quality random numbers with minimal computational overhead.

\[ \text{state}_{n+1} = (\text{state}_n \times \text{multiplier} + \text{increment}) \mod 2^{64} \]

However, unlike the LCG, PCG applies a permutation step to the generated state to enhance randomness:

\[ \text{output} = ((\text{state} >> \text{shift}) \oplus \text{state}) \mod 2^{64} \]

Where:

• state: The internal state of the generator.
• multiplier: A large odd number used in the linear step.
• increment: Another large odd number ensuring a full period.
• shift: A calculated bit-shift for scrambling the state.

The constructor initializes the generator's parameters: multiplier, increment, and modulus. Like previously, if a seed is provided, it is used as the initial state; otherwise, the generator creates a unique seed by combining the process ID and the current time:

    class PCGPseudoRand:
        def __init__(self, multiplier=6364136223846793005, increment=1442695040888963407, seed=None) -> None:
            self.multiplier = multiplier
            self.increment = increment | 1
            self.modulus = 2**64
            if seed is not None:
                self.state = seed
            else:
                self.state = int(os.getpid() + time.time()) & 0xFFFFFFFFFFFFFFFF
    

In this implementation the multiplier, increment and modulus are from Donald Knuth's MMIX.

The expression int(os.getpid() + time.time()) & 0xFFFFFFFFFFFFFFFF generates a pseudo-random 64-bit seed by combining the current process ID and the current time. This approach adds a level of non-determinism since the process ID and time are unlikely to repeat simultaneously (same way as LCG discussed earlier). The bitwise AND operation & 0xFFFFFFFFFFFFFFFF ensures that the resulting value is constrained to a 64-bit integer, even if the sum exceeds the 64-bit range.

The Linear Congruential Step method updates the internal state using the LCG formula:

    def _lcg_step(self):
        self.state = (self.state * self.multiplier + self.increment) & 0xFFFFFFFFFFFFFFFF
    

The new state is calculated by multiplying the current state by the multiplier, adding the increment, and applying a modulo operation to ensure the value stays within a 64-bit range.

Is it similar to an LCG? Yes, to some extent. However, there is an additional permutation step here, which does not exist in a standard LCG. This step involves the following:

1. The first step involves extracting key bits from the uppermost part of the 64-bit integer x to determine how much the number will be shifted. Specifically, shift = (x >> 59) + 1 performs a right shift of x by 59 bits. Since x is a 64-bit integer, this operation isolates the top 5 bits (bits 59-63) after the shift.

   The result is a value between 0 and 31, as the top 5 bits form a 5-bit binary number. Adding 1 ensures that the shift amount is always at least 1, preventing a shift by 0 bits, which would have no effect.

2. Next, the bits of x are shifted to the right by the calculated shift amount. This operation discards the bits on the right and shifts in zeros on the left. The resulting value is then XORed with the original value of x. The XOR operation compares each corresponding bit of the shifted value and x, returning 1 if the bits differ and 0 if they are the same. This step introduces a mix of the bits, ensuring better scrambling of the number.

        
    def _permute(self, x):
        shift = (x >> 59) + 1
        return ((x >> shift) ^ x) & 0xFFFFFFFFFFFFFFFF
    

Let's walk through the calculation step-by-step for the expression shift = (x >> 59) + 1, where the number x = 0xF123456789ABCDEF.

The hexadecimal value x = 0xF123456789ABCDEF is converted to binary:

        x = 0xF123456789ABCDEF

        In binary, x = 11110001001000110100010101100111000100111010111111011111
        

The expression x >> 59 means shifting the bits of x to the right by 59 positions. This operation discards the 59 least significant bits and moves the remaining bits to the right.

        x = 11110001001000110100010101100111000100111010111111011111

        x >> 59 = 00000000000000000000000000000000000000000000000000000111
        

Next, we compute the final value of shift. The top 5 bits of x after shifting right by 59 are 111 in binary, which equals 7 in decimal. Adding 1 ensures a minimum shift value of 1:

        shift = (x >> 59) + 1 = 7 + 1 = 8
    

Now, let’s calculate the expression ((x >> shift) ^ x) & 0xFFFFFFFFFFFFFFFF step by step.

Shifting x to the right by 8 positions discards the least significant 8 bits, shifting zeros into the left:

• The original value of x in binary: 11110001001000110100010101100111000100111010111111011111
• The shifted value (right shifted by 8 bits): 00000000111100010010001101000101011001110001001110101111

Performing the XOR operation between the original value x and the shifted value:

        Original (x):   11110001001000110100010101100111000100111010111111011111
        Shifted:        00000000111100010010001101000101011001110001001110101111
        XOR result:     11110001110100100110011000100010011101001011110001110000
    

In hexadecimal, the result of the XOR operation is: 0xF1D2662274BC70.

Finally, the result is ANDed with 0xFFFFFFFFFFFFFFFF to constrain it to 64 bits. Since the XOR result already fits within 64 bits, the AND operation does not change the value:

        ((x >> shift) ^ x) & 0xFFFFFFFFFFFFFFFF = 0xF1D2662274BC70
        

The final value is 0xF1D2662274BC70 or 68066805493841008 in decimal.

This value represents a 64-bit unsigned integer, making it suitable for scaling to a desired range using the scaling formula:

\[ \text{output} = \text{min_range} + \left( \frac{P(S_n)}{2^{64}} \right) \times (\text{max_range} - \text{min_range}) \]

• \( P(S_n) \): The 64-bit value to be scaled (e.g., 0xF1D2662274BC70).
• \( 2^{64} \): The total range of possible 64-bit unsigned integers.
• \( \text{min_range} \) and \( \text{max_range} \): The target range for scaling.

The fraction \( \frac{P(S_n)}{2^{64}} \) normalizes the 64-bit integer value P(S_n) into the range [0, 1). Multiplying this fraction by \( \text{max_range} - \text{min_range} \) scales it to the size of the target range. Adding \( \text{min_range} \) shifts the result to start at the desired minimum value.

3. Python implementations

    class LCGPseudoRand():
        def __init__(self, a=48271, c=0, m=2**31-1, seed=None) -> None: # C++ version
            self.a = a
            self.c = c
            self.m = m
    
            if seed:
                self.x0 = seed
            else:
                self.x0 = int(os.getpid() + time.time())
            
            self.prev_x = (a * self.x0 + self.c) % self.m
        
        def generate_number(self, min_range=None, max_range=None):
    
            self.prev_x = (self.a * self.prev_x + self.c) % self.m
    
            if min_range and max_range:
                return int((self.prev_x / self.m) * (max_range - min_range) + min_range)
            else:
                return self.prev_x

    class PCGPseudoRand:
        def __init__(self, multiplier=6364136223846793005, increment=1442695040888963407, seed=None) -> None:
            self.multiplier = multiplier  
            self.increment = increment | 1 
            self.modulus = 2**64 
    
            
            if seed is not None:
                self.state = seed
            else:
                self.state = int(os.getpid() + time.time()) & 0xFFFFFFFFFFFFFFFF
        
        def _lcg_step(self):
            self.state = (self.state * self.multiplier + self.increment) & 0xFFFFFFFFFFFFFFFF
    
        def _permute(self, x):
            shift = (x >> 59) + 1
            return ((x >> shift) ^ x) & 0xFFFFFFFFFFFFFFFF

        def generate_number(self, min_range=None, max_range=None):
            self._lcg_step()
            
            permuted_value = self._permute(self.state)
            print(self.state, permuted_value)
    
            if min_range is not None and max_range is not None:
                return int((permuted_value / self.modulus) * (max_range - min_range) + min_range)
            else:
                return permuted_value
                      
            

4. Considerations

While the Linear Congruential Generator (LCG) is faster and simpler to implement, it often suffers from lower randomness quality. LCGs generate sequences with a predictable structure.

The Permuted Congruential Generator (PCG), on the other hand, introduces a permutation step on top of the basic linear congruential approach. This permutation step greatly improves randomness quality by better scrambling the bits of the generated numbers. Although PCG incurs a slight computational overhead compared to LCG, it is negligible for most use cases and worth the tradeoff for significantly better randomness properties.

In the future, I hope to discuss the Mersenne Twister, another widely used pseudorandom number generator. Known for its extremely long period (2¹⁹⁹³⁷ − 1) and excellent randomness properties, the Mersenne Twister addresses many shortcomings of earlier generators like LCG. Interestingly, the Mersenne Twister is the default pseudorandom number generator in Python's random module. Its widespread use stems from its efficiency, ability to produce high-quality randomness, and suitability for general-purpose random number generation.