I was obviously not happy when out of all the interesting projects, I was assigned to improve the quality of correlated random numbers. I had started a new job, and this project was seen as relatively simple. I had never thought about that it could be intriguing. Let's talk about some high level concepts and algorithms related random numbers.

The background is, as we understand, we need random numbers to simulate financial markets, solve PDEs, perform encryption, model unknown variables like in weather predictions, in games, and so much more. In the real world, we presume, we have random events or sequences, like flipping a fair coin. The coin flip serves as the starting point for stochastic calculus and probability theory. In case you've ever questioned, computers are deterministic i.e. actions are determined by a prescribed set of instructions, how can it possibly produce a random number? When we say computer-generated random numbers, these are Pseudo Random numbers generated from some specific arithmetic operations. Ultimately, we are interested in sequences or say pseudorandom sequences.

In a simple term, we view a random sequence as a generator that gives a value when asked without us being able to predict that value based on any of the previous values.

A pseudo-random sequence is a deterministically generated sequence of numbers that is indistinguishable from an actual random sequence of numbers.

We need to test out if the pseudo-random sequence looks random enough. There are two popular tests called Chi-square test and spectral test. We will discuss these tests and cryptographically secure random number generator in the part 2 of this post.

Pseudo Random Numbers Generator (PRNG): Generating uniform random numbers

PRNGs are the algorithms that generate sequences of random numbers that approximate the property of randomness. The primary goal of the majority of pseudo-random number generators is to generate uniformly distributed pseudo-random integers throughout a specific range. PRNGs generate numbers fast and they can be reproduced by knowing the initial starting point (seed).

Linear congruential generator

Linear Congruential Generator is the most common and oldest algorithm for generating pseudo-randomized numbers. The generator is defined by the recurrence relation: $$ X_{n+1} = a(X_n + c) \space mod \space m$$

where X is the sequence of pseudo-random values, and

$$ m — the modulus$$

$$ 0<a<m — multiplier $$ $$0<c<m — increment $$ $$X_{0}<m — seed $$

PRNGs are fast, deterministic, and periodic. Periodic is not a desirable property although modern PRNGs have a quite long cycle.

XORshift

This is more robust and fast algorithm however it does not pass all the statistical tests without further improvements

/* Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" */
x = unit64_t
x = x^(x<<13)
x = x^(x>>17)
x = x^(x<<5)

We have family of xorshifts algorithms. A popular random number generator is xorshift128+.

The rust implementation of the high performance xoroshiro128+, xorshift128+, xorshift1024*, and splitmix64 pseudo random number generators can be found in xorshift create. docs.rs/xorshift/latest/xorshift

Mersenne Twister (MT)

This is another standard PRNG, at least for tasks that don't require a cryptographically secure generator. It is one of the fastest pseudo-random number generators that passes almost all statistical tests. However, observing a couple hundred outputs, it is possible to predict all future outputs. The Mersenne Twister has a period of $$2^{19937} − 1$$, a Mersenne prime, hence the “Mersenne” part of the name. The implementation of MT in rust can be found in crate mersenne_twister.

use mersenne_twister::MersenneTwister;
use rand::{Rng, SeedableRng};

fn main() {
    // Get a seed somehow.
    let seed: u64 = 0x123456789abcdef;
    // Create the default RNG.
    let mut rng: MersenneTwister = SeedableRng::from_seed(seed);
}

Generating Uniform Random Numbers in Rust

A group of crates makes up the Rand library. The Rand crate provides the main user interface. The library contains several building blocks: getrandom interfaces with the platform-dependent random number source, rand_core defines the API that generators must implement, and a number of crates like rand_chacha and rand_xoshiro provide pseudo-random generators.

getrandom ┐
          └ rand_core ┐
                      ├ rand_chacha ┐
                      ├ rand_hc     ┤
                      ├ rand_pcg    ┤
                      └─────────────┴ rand ┐
                                           ├ rand_distr
                                           └ statrs
from https://rust-random.github.io/book/crates.html

Deterministic generators in Rust

The following rust crates implement pseudo-random number generators:

1. rand_chacha provides generators using the ChaCha cipher
2. rand_hc implements a generator using the HC-128 cipher
3. rand_isaac implements the ISAAC generators
4. rand_pcg implements a small selection of PCG generators
5. rand_xoshiro implements the SplitMix and Xoshiro generators
6. rand_xorshift implements the basic Xorshift generator

We can use rand create to generate the random numbers in Rust.

    let mut rng = rand_chacha::ChaCha8Rng::seed_from_u64(555);
    println!("{}", rng.gen::<i32>()); // prints 425636761
    let mut rng1 = rand_pcg::Lcg64Xsh32::seed_from_u64(555);
    println!("{}", rng1.gen::<i32>()); // prints-1224090965

The extended documentation for Rust's Random number library can be found in The rust Rand book.

Generating Standard Normal Random variables in Rust

Let's assume that we finally figured out how to generate uniformly distributed random numbers. Now we also need standard normal random numbers for our simulations. There are methods to generate standard normal distributed random numbers from a uniform distribution such as Marsaglia polar method, Box–muller transform, the Ziggurat algorithm, and more. The ziggurat algorithm is the most efficient of these three and is used by Rust to generate standard normal random numbers. Marsaglia polar and Box-Mullar transform methods are outdated its good for our understanding.

Marsaglia polar method

The polar method works by choosing random points (x, y) in the square −1 < x < 1, −1 < y < 1 until

$$ 0<s=x^{2}+y^{2}<1,$$, and then pair of normal random variables are,

$$ x{\sqrt {\frac {-2\ln(s)}{s}}}$$, $$ y{\sqrt {\frac {-2\ln(s)}{s}}}$$

This method uses the fact that if the random variable U is U(0, 1) then the random variable X defined by X = 2U − 1 is U(−1, 1). We now choose two variables defined by $$ X_j = 2 * U_j -1, $$ $$ U_j \sim U(0,1), j = 1,2 $$ then we define, $$ S = X^2 + Y^2 <=1, $$ We try different values till the above inequality is satisfied. Once we find the S, we calculate the pair of normal random numbers as in the above formula.

fn get_normal_marsagliapolar() -> (f64, f64) {
    let mut rng = rand::thread_rng();
    let mut X = 0.0;
    let mut Y = 0.0;
    let mut S = 0.0f64;
    while(true) {
        X = Uniform::new(0.0,1.0).sample(&mut rng)*2.0 -1.0;
        Y = Uniform::new(0.0,1.0).sample(&mut rng)*2.0 -1.0;
        S = X*X + Y*Y;
        if S<1.0f64 && S != 0.0f64 {
            break;
        }
    }
    let I = ((-2.0 * S.ln()) / S).sqrt();
    (I*X,I*Y)
}

Box-Muller method

The Box–Muller transform was developed as a more computationally efficient alternative to the inverse transform sampling method, which we will save for later discussion. Box-muller is similar to Marsaglia polar method, which uses polar coordinates instead of Cartesian coordinates. This method is based on the observation that if r and ϕ are two independent U(0, 1) random variables then the variables.

$$ N1={\sqrt {-2\ln r}}\cos(2\pi ϕ) $$ $$ N2={\sqrt {-2\ln r}}\sin(2\pi ϕ) $$

fn get_normal_boxmullar() -> (f64, f64) {
    let mut rng = rand::thread_rng();
    let r:f64 = Uniform::new(0.0,1.0).sample(&mut rng);
    let p:f64 = Uniform::new(0.0,1.0).sample(&mut rng);
    let tmp:f64 = (-2.0*r.ln()).sqrt();
    (tmp*cos(p*2.0*PI),tmp*sin(p*2.0*PI))
}

Conclusion

"Anyone who considers arithmetic methods of producing random digits is, of course, in a state of sin." John von Neumann

Computer scientists have taken a long journey to generate pseudo-random numbers from Linear congruential generators, XORshift to Chacha, PCG, etc. It is useful to have some background knowledge and the introduction to PRNGs. Finally, we reviewed Box-mullar and Marsaglia-polar methods to generate standard normal random numbers from uniform random numbers. Hopefully this is a good introduction.

Useful references

pcg-random.org

prng.di.unimi.it