Random Numbers in F#

Many problems in engineering, finance and statistics can not be solved by direct methods, but a great number of them can be solved approximately using randomized algorithms. All those algorithms need flexible and efficient pseudo-random number generators. An effective implementation of PRNG in the F# language is somewhat tricky.

The F# language is based on .Net framework, which already provides a fast pseudo-random generator class System.Random, but it meets only certain statistical requirements for randomness. However, there are two classes of generators that have good performance:

• Combined multiple recursive generators (e.g., L’Ecuyer’s MRG32k3a)
• Twisted general feedback shift register generators (e.g., Matsumoto and Nishimura’s MT19937ar)

Despite the fact that MT19937ar has longer period than MRG32k3a, the state space in MT19937ar is very large and it is somewhat slow. Thus, MRG32k3a becomes a better alternative for many use cases.

Algorithm

MRG32k3a is just one of several good parameters sets that L’Ecuyer obtained for CMRGs ¹. It has two components of order 3. The state at step \(n\) is a pair of vectors $S_{1,n}=(s_{1,n},s_{1,n+1}, s_{1,n+2})$ and $S_{2,n}=(s_{2,n},s_{2,n+1}, s_{2,n+2})$, which change according to the linear recurrences:

The uniformly distributed output \(u_n\) is obtained by

Implementation

Two key requirement for the random number generator are the computational performance and minimal memory consumption. MRG32k3a algorithm requires only a pair of 3 element double valued arrays for storing the current generator state, but F# arrays are heap allocated. Fortunately, the F# language supports structure types, which can be allocated on the stack:

[<Struct>]
type State32k3a =
  val mutable S10 : double
  val mutable S11 : double
  val mutable S12 : double

  val mutable S20 : double
  val mutable S21 : double
  val mutable S22 : double

  new (s10, s11, s12, s20, s21, s22) =
    { S10 = s10; S11 = s11; S12 = s12;
      S20 = s20; S21 = s21; S22 = s22 }

The implementation of MRG32k3a gets a reference to State32k3a structure and updates its fields. In this way no intermediate object allocation and copy happens, and the computational performance is close to the C implementation.

let inline private (%.) p m =
  let k = int32(p / m)
  let r = p - double(k) * m
  if r < 0.0 then r + m else r

let MRG32k3a (state : State32k3a byref) =
  let p1 = (1403580.0 * state.S11 - 810728.0 * state.S10) %. 4294967087.0
  state.S10 <- state.S11; state.S11 <- state.S12; state.S12 <- p1
  let p2 = (527612.0 * state.S21 - 1370589.0 * state.S20) %. 4294944443.0
  state.S20 <- state.S21; state.S21 <- state.S22; state.S22 <- p2
  if p1 <= p2 then p1 - p2 + 4294967087.0 else p1 - p2

Nonetheless, it is possible to create a nice functional wrapper around the low level implementation:

let makeMRG32k2a (seed : int32) =
  let seed = float(seed)
  let state = ref (State32k3a(seed, seed, seed, seed, seed, seed))
  let gen = (fun () -> let mutable s = !state in let r = (MRG32k3a &s) in state := s; r)
  gen() |> ignore; gen() |> ignore; // warm up the generator
  gen

This implementaion returns a closure that can be used for generating a sequence of random numbers:

let gen = MRG.defaultMRG32k2a

for i in 1..10 do
  printfn "%A" (gen())

Though, the price for this nice wrapper is high enough. It involves a heap allocated box for the structure and additional copy operations for reading/writing the structure from the box.

Conslusions

Despite the fact that pure functions and immutable structures are preferable in the most cases because they reduce the number of mistakes, it is important to have clearly designed low level functions for advanced optimizations. Unlike R, Python and other scripting languages, F# has no language border between high level and low level implementations, and, thus, facilitates reuse of low level implementations.