Mersenne Twister

// Mersenne Twister.
// This code is placed into the public domain.

#light

open UInt32

type state = { mt : uint32 array; mutable mti : int; }

let n = 624
let m = 397
let upper_mask = 0x80000000u
let lower_mask = 0x7fffffffu
let mag01 = [| 0u; 0x9908b0dfu |]
let mask_b = 0x9d2c5680u
let mask_c = 0xefc60000u

let init_state () = { mt = Array.create n 0u; mti = n+1; }

let state = init_state ()

let init_genrand seed =
  let ml = 1812433253u
  let mfill = 0xffffffffu
  state.mt.[0] <- (seed &&& mfill)
  for i = 1 to n - 1 do
    let mtprev = state.mt.[i - 1]
    state.mt.[i] <- (ml * (mtprev ^^^ (mtprev >>> 30)) + (of_int i))

let init_by_array key =
  init_genrand 19650218u
  let i = ref 1
  let j = ref 0
  let k = ref (max n (Array.length key))
  let mval1 = 1664525u
  let mval2 = 1566083941u

  while !k > 0 do
    let mtprev = state.mt.[!i - 1]
    state.mt.[!i] <- (state.mt.[!i] ^^^ ((mtprev ^^^ (mtprev >>> 30)) * mval1)) + key.[!j] + (of_int !j)
    incr i
    incr j
    if !i >= n then
      state.mt.[0] <- state.mt.[n - 1]
      i := 1
    if !j >= Array.length key then j := 0
    decr k

  k := n - 1

  while !k > 0 do
    let mtprev = state.mt.[!i - 1]
    state.mt.[!i] <- (state.mt.[!i] ^^^ ((mtprev ^^^ (mtprev >>> 30)) * mval2)) - (of_int !i)
    incr i
    if !i >= n then
      state.mt.[0] <- state.mt.[n - 1]
      i := 1
    decr k
    state.mt.[0] <- 0x80000000u

let genrand_int32 () =
    if (state.mti >= n) then
        if (state.mti = n+1) then init_genrand 5489u
        let y = ref 0u
        for i = 0 to n - m - 1 do
          y := ((state.mt.[i] &&& upper_mask) ||| (state.mt.[i + 1] &&& lower_mask))
          state.mt.[i] <- state.mt.[i + m] ^^^ (!y >>> 1) ^^^ (mag01.[to_int (!y &&& one)])

        for i = n - m to n - 2 do
          y := (state.mt.[i] &&& upper_mask) ||| (state.mt.[i + 1] &&& lower_mask)
          state.mt.[i] <- state.mt.[i + m - n] ^^^ (!y >>> 1) ^^^ (mag01.[to_int (!y &&& one)])

        y := (state.mt.[n - 1] &&& upper_mask) ||| (state.mt.[0] &&& lower_mask)
        state.mt.[n - 1] <- state.mt.[m - 1] ^^^ (!y >>> 1) ^^^ (mag01.[to_int (!y &&& one)])
        state.mti <- 0

    let y = ref (state.mt.[state.mti])
    y := !y ^^^ (!y >>> 11)
    y := !y ^^^ ((!y <<< 7) &&& 0x9d2c5680u)
    y := !y ^^^ ((!y <<< 15) &&& 0xefc60000u)
    y := !y ^^^ (!y >>> 18)
    state.mti <- state.mti + 1
    !y

let genrand_int31 () = (to_int32 (genrand_int32 () >>> 1))
let genrand_real1 () = (System.Convert.ToDouble (genrand_int32 ())) * (1.0 / 4294967295.0)
let genrand_real2 () = (System.Convert.ToDouble (genrand_int32 ())) * (1.0 / 4294967296.0)
let genrand_real3 () = ((System.Convert.ToDouble (genrand_int32 ())) + 0.5) * (1.0 / 4294967296.0)

// 1000 outputs of genrand_int32
let test0 () =
  for i = 0 to 999 do
    Printf.printf "%10A " (genrand_int32 ())
    if (i % 5 = 4) then Printf.printf "\n"

// 1000 outputs of genrand_int31
let test1 () =
  for i = 0 to 999 do
    Printf.printf "%10d " (genrand_int31 ())
    if (i % 5 = 4) then Printf.printf "\n"

// 1000 output of genrand_real2
let test2 () =
  for i = 0 to 999 do
    Printf.printf "%10.8f " (genrand_real2 ())
    if (i % 5 = 4) then Printf.printf "\n"