Pastebin launched a little side project called VERYVIRAL.com, check it out ;-) Want more features on Pastebin? Sign Up, it's FREE!
Guest

Mersenne Twister

By: a guest on Jul 22nd, 2011  |  syntax: F#  |  size: 3.25 KB  |  views: 696  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. // Mersenne Twister.
  2. // This code is placed into the public domain.
  3.  
  4. #light
  5.  
  6. open UInt32
  7.  
  8. type state = { mt : uint32 array; mutable mti : int; }
  9.  
  10. let n = 624
  11. let m = 397
  12. let upper_mask = 0x80000000u
  13. let lower_mask = 0x7fffffffu
  14. let mag01 = [| 0u; 0x9908b0dfu |]
  15. let mask_b = 0x9d2c5680u
  16. let mask_c = 0xefc60000u
  17.    
  18. let init_state () = { mt = Array.create n 0u; mti = n+1; }
  19.  
  20. let state = init_state ()
  21.  
  22. let init_genrand seed =
  23.   let ml = 1812433253u
  24.   let mfill = 0xffffffffu
  25.   state.mt.[0] <- (seed &&& mfill)
  26.   for i = 1 to n - 1 do
  27.     let mtprev = state.mt.[i - 1]
  28.     state.mt.[i] <- (ml * (mtprev ^^^ (mtprev >>> 30)) + (of_int i))
  29.  
  30. let init_by_array key =
  31.   init_genrand 19650218u
  32.   let i = ref 1
  33.   let j = ref 0
  34.   let k = ref (max n (Array.length key))
  35.   let mval1 = 1664525u
  36.   let mval2 = 1566083941u
  37.  
  38.   while !k > 0 do
  39.     let mtprev = state.mt.[!i - 1]
  40.     state.mt.[!i] <- (state.mt.[!i] ^^^ ((mtprev ^^^ (mtprev >>> 30)) * mval1)) + key.[!j] + (of_int !j)
  41.     incr i
  42.     incr j
  43.     if !i >= n then
  44.       state.mt.[0] <- state.mt.[n - 1]
  45.       i := 1
  46.     if !j >= Array.length key then j := 0
  47.     decr k
  48.    
  49.   k := n - 1
  50.  
  51.   while !k > 0 do
  52.     let mtprev = state.mt.[!i - 1]
  53.     state.mt.[!i] <- (state.mt.[!i] ^^^ ((mtprev ^^^ (mtprev >>> 30)) * mval2)) - (of_int !i)
  54.     incr i
  55.     if !i >= n then
  56.       state.mt.[0] <- state.mt.[n - 1]
  57.       i := 1
  58.     decr k
  59.     state.mt.[0] <- 0x80000000u
  60.  
  61. let genrand_int32 () =
  62.     if (state.mti >= n) then
  63.         if (state.mti = n+1) then init_genrand 5489u
  64.         let y = ref 0u
  65.         for i = 0 to n - m - 1 do
  66.           y := ((state.mt.[i] &&& upper_mask) ||| (state.mt.[i + 1] &&& lower_mask))
  67.           state.mt.[i] <- state.mt.[i + m] ^^^ (!y >>> 1) ^^^ (mag01.[to_int (!y &&& one)])
  68.  
  69.         for i = n - m to n - 2 do
  70.           y := (state.mt.[i] &&& upper_mask) ||| (state.mt.[i + 1] &&& lower_mask)
  71.           state.mt.[i] <- state.mt.[i + m - n] ^^^ (!y >>> 1) ^^^ (mag01.[to_int (!y &&& one)])
  72.  
  73.         y := (state.mt.[n - 1] &&& upper_mask) ||| (state.mt.[0] &&& lower_mask)
  74.         state.mt.[n - 1] <- state.mt.[m - 1] ^^^ (!y >>> 1) ^^^ (mag01.[to_int (!y &&& one)])
  75.         state.mti <- 0
  76.  
  77.     let y = ref (state.mt.[state.mti])
  78.     y := !y ^^^ (!y >>> 11)
  79.     y := !y ^^^ ((!y <<< 7) &&& 0x9d2c5680u)
  80.     y := !y ^^^ ((!y <<< 15) &&& 0xefc60000u)
  81.     y := !y ^^^ (!y >>> 18)
  82.     state.mti <- state.mti + 1
  83.     !y
  84.  
  85. let genrand_int31 () = (to_int32 (genrand_int32 () >>> 1))
  86. let genrand_real1 () = (System.Convert.ToDouble (genrand_int32 ())) * (1.0 / 4294967295.0)
  87. let genrand_real2 () = (System.Convert.ToDouble (genrand_int32 ())) * (1.0 / 4294967296.0)
  88. let genrand_real3 () = ((System.Convert.ToDouble (genrand_int32 ())) + 0.5) * (1.0 / 4294967296.0)
  89.  
  90. // 1000 outputs of genrand_int32
  91. let test0 () =
  92.   for i = 0 to 999 do
  93.     Printf.printf "%10A " (genrand_int32 ())
  94.     if (i % 5 = 4) then Printf.printf "\n"
  95.  
  96. // 1000 outputs of genrand_int31
  97. let test1 () =
  98.   for i = 0 to 999 do
  99.     Printf.printf "%10d " (genrand_int31 ())
  100.     if (i % 5 = 4) then Printf.printf "\n"
  101.  
  102. // 1000 output of genrand_real2
  103. let test2 () =
  104.   for i = 0 to 999 do
  105.     Printf.printf "%10.8f " (genrand_real2 ())
  106.     if (i % 5 = 4) then Printf.printf "\n"