let gamma_x_infty =
    let nb_iterations_max = 1000 in
    let relative_accuracy = 1.e-7 in (* optimizable compared to its usage *)
    let min_float = 1.e-200 in (* current caml implementation is e-308 *)
  fun a x ->
    let b = ref (x +. 1. -. a) in
    let c = ref (1. /. min_float) in
    let d = ref (1. /. !b) in
    let h = ref !d in
    try
      for i=1 to nb_iterations_max do
        let i = float i in
        let an = -. i *. (i -. a) in
        b := !b +. 2.;
        d := an *. !d +. !b;
        if abs_float !d < min_float then d := min_float;
        c := !b +. an /. !c;
        if abs_float !c < min_float then c := min_float;
        d := 1. /. !d;
        let coeff = !d *. !c in
        h := !h *. coeff;
        if abs_float (coeff -. 1.) < relative_accuracy then raise Converged;
      done;
      failwith (Printf.sprintf "gamma_x_infty. Parameters: a=%f x=%f" a x);
    with Converged -> exp (-. x +. a *. log x -. log_gamma a) *. !h