def log_factorial(n, a): if n == 0: return 0. return scipy.special.gammaln(n + a) - scipy.special.gammaln(a)