Computation of the Gamma, Digamma, and Trigamma Functions

Abstract

This paper gives an approximation for the gamma function that, while different, has the same form as one by Lanczos [SIAM J. Numer. Anal., B1 (1964), pp. 86–96]. Both approximations correct Stirling’s approximation with contributions from the gamma function’s poles, and both require $O( - \log \epsilon )$ time independent of z to calculate $z!$ with a relative error $\epsilon $. At comparable accuracies, this paper’s approximation requires slightly more computation than Lanczos’ but is clearly superior in two ways; its coefficients are given by simple formulas and the error estimates apply not only to the gamma function, but also to its derivatives. Thus approximations for the digamma and trigamma functions are also given. Let $F_{a,{1}/ {2}} (z) = z!(z + a)^{ - z - ({1 / 2})} e^{z + a} (2\pi )^{ - ({1 / 2})} $ and $f_a (z) = \ln F_{a,{1 / 2}} (z)$, with z complex and a real. Several lemmas in the paper require “Stirling inequalities,” i.e., bounds on $| {F_{a,{1 / 2}} (z)} |$. To this end, the classical Binet integral for $f_0 (z)$ is generalized from $a = 0$ to all real a. Since the generalized Binet integral is a Laplace transform, the requisite Stirling inequalities follow from the complete monotonicity of such transforms. The proof of the main theorem uses the calculus of residues and presents a lemma generalizing Plana’s theorem, which, like the Euler–Maclaurin sum formula, evaluates the difference between a sum and the corresponding integral. This paper also presents a factorial approximation that is both simpler and more accurate than Stirling’s approximation.