let cnd x =
(* Modified from C++ code by David Koppstein. Found from
   www.sitmo.com/doc/Calculating_the_Cumulative_Normal_Distribution *)
  let b1,b2,b3,b4,b5,p,c =
    0.319381530, -0.356563782, 1.781477937, -1.821255978,
    1.330274429, 0.2316419, 0.39894228 in
  if x >= 0. then
    let t = 1. /. (1. +. (p *. x)) in
    (1. -. (c *. (exp (-.x *. x /. 2.)) *. t *.
      (t *. (t *. (t *. ((t *. b5) +. b4) +. b3) +. b2) +. b1 )))
    else
      let t = 1. /. (1. -. p *. x) in
      c *. (exp (-.x *. x /. 2.)) *. t *.
        (t *. (t *. (t *. ((t *. b5) +. b4) +. b3) +. b2) +. b1 )