let gamma_0_x =
    let nb_iterations_max = 1000 in
    let relative_accuracy = 1.e-7 in (* optimizable compared to its usage *)
  fun a x ->
    assert (x >= 0.);
    let a_plus_n = ref a in
    let coeff = ref (1. /. a) in
    let sum = ref !coeff in
    try
      for i=1 to nb_iterations_max do
        a_plus_n := !a_plus_n +. 1.;
        coeff := !coeff *. x /. !a_plus_n;
        sum := !sum +. !coeff;
        if !coeff < !sum *. relative_accuracy then raise Converged;
      done;
      failwith (Printf.sprintf "gamma_0_x. Parameters: a=%f x=%f" a x);
    with Converged -> !sum *. exp (-.x +. a *. log x -. log_gamma a)