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:

  def a : u32 = 1103515245
  def c : u32 = 12345
  def m : u32 = 1<<31

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

  def 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.

  def 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.

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

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

This is how we might use it:

def 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]:

def 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]:

def 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.

Random numbers

See also