let ltqnorm p =
(*
   Modified from python code by David Koppstein. Original comments follow below.
   First version was written in perl, by Peter J. Acklam, 2000-07-19.
   Second version was ported to python, by Dan Field, 2004-05-03.
*)
    if (p <= 0.) || (p >= 1.) then
      raise (ValueError ("Argument to ltqnorm " ^ (Float.to_string p) ^
        " must be in open interval (0,1)"))
    else
      (* Coefficients in rational approximations. *)
      let a =
        [|-3.969683028665376e+01; 2.209460984245205e+02;
        -2.759285104469687e+02; 1.383577518672690e+02;
        -3.066479806614716e+01; 2.506628277459239e+00|] in
      let b =
        [|-5.447609879822406e+01; 1.615858368580409e+02;
        -1.556989798598866e+02; 6.680131188771972e+01;
        -1.328068155288572e+01|] in
      let c =
        [|-7.784894002430293e-03; -3.223964580411365e-01;
        -2.400758277161838e+00; -2.549732539343734e+00;
        4.374664141464968e+00; 2.938163982698783e+00|] in
      let d =
        [|7.784695709041462e-03; 3.224671290700398e-01;
        2.445134137142996e+00; 3.754408661907416e+00|] in

      (* Define break-points. *)
      let plow  = 0.02425 in
      let phigh = 1. -. plow in
      let f q =
        (((((c.(0)*.q+.c.(1))*.q+.c.(2))*.q+.c.(3))*.q+.c.(4))*.q+.c.(5)) /.
          ((((d.(0)*.q+.d.(1))*.q+.d.(2))*.q+.d.(3))*.q+.1.)
      in

    (* Rational approximation for lower region: *)
    if p < plow then
      let q  = sqrt ((-2.) *. (log p)) in
      f q

    (* Rational approximation for upper region: *)
    else if phigh < p then
      let q = sqrt ((-2.) *. (log (1. -. p))) in
      f q

    (* Rational approximation for central region: *)
    else
      let q = p -. 0.5 in
      let r = q *. q in
    (((((a.(0)*.r+.a.(1))*.r+.a.(2))*.r+.a.(3))*.r+.a.(4))*.r+.a.(5))*.q /.
      (((((b.(0)*.r+.b.(1))*.r+.b.(2))*.r+.b.(3))*.r+.b.(4))*.r+.1.)