let log_gamma =
    let coeffs = [|76.18009172947146; -86.50532032941677; 24.01409824083091;
                   -1.231739572450155; 0.1208650973866179e-2; -0.5395239384953e-5|] in
    fun x ->
      assert (x > 0.);
      let y = ref x in
      let temp = x +. 5.5 in
      let temp = (x +. 0.5) *. log temp -. temp in
      let accu = ref 1.000000000190015 in
      for i = 0 to 5 do
        y := !y +. 1.;
        accu := !accu +. coeffs.(i) /. !y
      done;
      temp +. log (2.5066282746310005 *. !accu /. x)