module RND = (struct (* Generate Normal Random Deviates *)
  let p = ref false
  and v = ref 0.
  and dr _ = (Random.float 2.) -. 1.
  let rec rnd d = if !p then (p := false; !v) else (
    let x = dr () and y = dr () 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 rnd : unit -> float end);;