SciMark2 Monte Carlo in OCaml

module Random = struct
  let m = Array.create 17 0
  let i = ref 4
  let j = ref 16
  let mdig = 32
  let m1 = (1 lsl (mdig - 2)) + (1 lsl (mdig - 2)) - 1
  let m2 = 1 lsl (mdig / 2)
  let dm1 = 1.0 /. (2. ** float(mdig-1) -. 1.)

  let nextDouble() =
    let k = m.(!i) - m.(!j) in
    let k = if k < 0 then k + m1 else k in
    m.(!j) <- k;
    i := if !i = 0 then 16 else !i - 1;
    j := if !j = 0 then 16 else !j - 1;
    if k>=0 then dm1 *. float k else 1. +. dm1 *. float k

  let initialize seed =
    let jseed = max (abs seed) m1 in
    let jseed = ref (if jseed land 1 = 0 then jseed - 1 else jseed) in
    let k0 = 9069 mod m2 in
    let k1 = 9069 / m2 in
    let j0 = ref (!jseed mod m2) in
    let j1 = ref (!jseed / m2) in
    for i=0 to 16 do
      jseed := !j0 * k0;
      j1 := (!jseed / m2 + !j0 * k1 + !j1 * k0) mod (m2 lsr 1);
      j0 := !jseed mod m2;
      m.(i) <- !j0 + (m2 * !j1)
    done;
    i := 4;
    j := 16

  let () =
    initialize 101010
end

module MonteCarlo = struct
  let flops n = 4 * n

  let exec n =
    let c = ref 0 in
    let sqr x = x *. x in
    for i=1 to n do
      if sqr(Random.nextDouble()) +. sqr(Random.nextDouble()) <= 1. then
	incr c
    done;
    4. *. float !c /. float n
end

let () =
  Printf.printf "%0.15f\n%!" (MonteCarlo.exec 100000000)