Fork me on GitHub
Source file: random-numbers.fut

# Random numbers

Generating random numbers in a purely functional language requires us to manually track the state of the random number generator (RNG). The following code shows the idea. Note that the random numbers we will generate are of mediocre quality at best, but they’re good enough for small visualisations and simulations and such. Random number generation is a very complex subject, but most of the complexities are the same in any language. Here we focus only on the ones that are particular to Futhark.

First we define a module type describing the interface of a random number generator.

``module type rand = {``

It consists of a random number generator state:

``  type rng``

These states can be initialised from a seed:

``  val init : i32 -> rng``

To generate a uniformly distributed random integer, we pass in an RNG state, which will give us back an integer and a new state:

``  val rand : rng -> (rng, i32)``

It is the user’s responsibility not to use the same state more than once (unless that is what is desired).

We will sometimes need to use one seed to generate multiple random streams. The `split` function splits one RNG state in twain:

``  val split : rng -> (rng, rng)``

And `split_n` splits it arbitrary many ways, which is useful when we wish to compute many parallel streams in parallel:

``````  val split_n : (n: i64) -> rng -> [n]rng
}``````

We can implement the `rand` module type using various different algorithms. For this example, we’ll go with a simple linear congruential generator (LCG):

``module lcg : rand = {``

The RNG state is just 32 bits.

``  type rng = u32``

An LCG is defined by the three constants a, c, and m. We use the following values, which are the same that are used by glibc:

``````  let a : u32 = 1103515245
let c : u32 = 12345
let m : u32 = 1<<31``````

Note that with these constants, we will never generate a negative random integer.

``````  let rand rng =
let rng' = (a * rng + c) % m
in (rng', i32.u32 rng')``````

Initialisation consists of a few rounds of a hash function found on StackOverflow, as is tradition.

``````  let init (x: i32): u32 =
let x = u32.i32 x
let x = ((x >> 16) ^ x) * 0x45d9f3b
let x = ((x >> 16) ^ x) * 0x45d9f3b
let x = ((x >> 16) ^ x)
in x``````

The main weakness of this initialisation is that a zero seed will result in a zero state, which will be “stuck”, producing only zeroes. So don’t do that.

The splitting functions creates RNG states by slightly twiddling the states.

``````  let split (rng: rng) =
(init (i32.u32 (rand rng).0),
init (i32.u32 rng))

let split_n n rng =
tabulate n (\i -> init (i32.u32 rng ^ i32.i64 i))
}``````

This is how we might use it:

``````let test seed =
let rng = lcg.init seed
let (rng, x) = lcg.rand rng
let (_, y) = lcg.rand rng
in (x, y)``````

Usually we don’t need integers in the full `i32` range, so let’s define a handy function for generating integers in the range `[0, bound-1]`:

``````let rand_i32 (rng: lcg.rng) (bound: i32) =
let (rng, x) = lcg.rand rng
in (rng, x % bound)``````

Note that the resulting integers are not fully uniformly distributed unless the span of the RNG is divisible by `bound`. In practice, as long as `bound` isn’t too large, it should not matter much.

We can also define a function that gives us a random `f64` in the interval `[0,1]`:

``````let rand_f64 (rng: lcg.rng) =
let (rng, x) = lcg.rand rng
in (rng, f64.i32 x / f64.i32 i32.highest)``````

For a real random number package, we’d then go on defining various other function for generating normally distributed numbers. We would also write these functions as part of a parametric module, rather than hard-coding the use of `lcg`. All this is done in the cpprandom package, which you should use for real Futhark programs.