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) =
0),
(init (i32.u32 (rand rng).
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.
See also
The cpprandom package.