hubFS: THE place for F#

Just in case anyone has any need for it, for whatever reason or purpose, I'm posting a more or less direct, non-idiomatic F# conversion of a C version of the Mersenne Twister random number generator. Last time I checked, creating a sequence of 1,000,000 random numbers using both the F# and C versions with the same starting seed generated identical sequences.

The original code can be found at:

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c

///
///   mt.fs
///

open Printf

let N = 624
let M = 397
let MATRIX_A = 0x9908b0dful
let UPPER_MASK = 0x80000000ul
let LOWER_MASK = 0x7ffffffful

let mt = Array.create N 0ul
let mti = ref (N+1)

let init_genrand s =
  mt.(0) <- (s &&& 0xfffffffful);
  for i = 1 to N-1 do
    mt.(i) <- (1812433253ul * (mt.(i-1) ^^^ (mt.(i-1) >>> 30)) + (UInt32.of_int i));
    mt.(i) <- mt.(i) &&& 0xfffffffful
  done;
  mti := N

let init_by_array init_key =
  let i = ref 1
  and j = ref 0
  and k = ref (if N > Array.length init_key then N else Array.length init_key) in
    init_genrand 19650218ul;
    while (!k > 0) do
      mt.(!i) <- (mt.(!i) ^^^ ((mt.(!i-1) ^^^ (mt.(!i-1) >>> 30)) * 1664525ul)) +
        init_key.(!j) + (UInt32.of_int !j);
      mt.(!i) <- mt.(!i) &&& 0xfffffffful;
      incr i; incr j;
      if !i >= N then
      begin
        mt.(0) <- mt.(N-1);
        i := 1
      end;
      if !j >= Array.length init_key then j := 0;
      decr k
    done;
    k := N-1;
    while (!k > 0) do
      mt.(!i) <- (mt.(!i) ^^^ ((mt.(!i-1) ^^^ (mt.(!i-1) >>> 30)) * 1566083941ul)) -
        (UInt32.of_int !i);
      mt.(!i) <- mt.(!i) &&& 0xfffffffful;
      incr i;
      if (!i >= N) then
      begin
        mt.(0) <- mt.(N-1);
        i := 1
      end;
      decr k
    done;
    mt.(0) <- 0x80000000ul

let genrand_int32 () =
  let mag01 = [|0x0ul; MATRIX_A|] in
    if !mti >= N then
    begin
      if (!mti = N+1) then init_genrand 5489ul;
   
      for kk = 0 to N-M-1 do
        let y = (mt.(kk) &&& UPPER_MASK) ||| (mt.(kk+1) &&& LOWER_MASK) in
          mt.(kk) <- mt.(kk+M) ^^^ (y >>> 1) ^^^ mag01.(UInt32.to_int (y &&& 0x1ul))
      done;
     
      for kk = N-M to N-2 do
        let y = (mt.(kk) &&& UPPER_MASK) ||| (mt.(kk+1) &&& LOWER_MASK) in
          mt.(kk) <- mt.(kk+(M-N)) ^^^ (y >>> 1) ^^^ mag01.(UInt32.to_int (y &&& 0x1ul))
      done;
     
      let y = (mt.(N-1) &&& UPPER_MASK) ||| (mt.(0) &&& LOWER_MASK) in
        mt.(N-1) <- mt.(M-1) ^^^ (y >>> 1) ^^^ mag01.(UInt32.to_int (y &&& 0x1ul));
        mti := 0
    end;
   
    let y = ref (mt.(!mti)) in
      incr mti;
      y := !y ^^^ (!y >>> 11);
      y := !y ^^^ ((!y <<< 7) &&& 0x9d2c5680ul);
      y := !y ^^^ ((!y <<< 15) &&& 0xefc60000ul);
      y := !y ^^^ (!y >>> 18);
      !y

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

 

Thanks for that, it'll be very usefull. I started to convert an OCAML version a while back but it was a bit much for an functional beginner :)