module RND = (struct (* Generate Random Normal Deviates *)
  open Random
  let p = ref false
  let v = ref 0.
  let init n = init n; p := false
  let rec rnd _ = if !p then (p := false; !v) else (
    let x = (float 2.) -. 1. and y = (float 2.) -. 1. in
    let r2 = x *. x +. y *. y in
    if r2 < 1. then (
      let q = sqrt (-2. *. (log r2) /. r2) in 
      v := x *. q;
      p := true;
      y*.q)
    else rnd ())
end : sig
 val init : int -> unit
 val rnd : unit -> float
end);;