A Note on a Method of Computing the Gamma Function

Abstract

Numerous formulas are available for the computation of the Gamma function [1, 2]. The purpose of this note is to indicate the value of a well-known method that is easily extended for higher accuracy requirements. Using the recursion formula for the Gamma function, Γ( x + 1) = x Γ( x ), (1) and Stirling's asymptotic expansion for ln Γ( x ) [3], we have ln Γ( x ) ∼ ( x - 1/2) ln x - x + 1/2 ln 2π + ∑ N r =1 C r / x 2 r -1 . (2) It follows that, if k and N are appropriately selected positive integers, Γ( x + 1) can be represented by Γ( x + 1) ∼ √2π exp ( x + k - 1/2) ln ( x + k ) - ( x + k ) exp ∑ N r =1 C r /( x + k ) 2 r -1 /( x + 1)( x + 2) … ( x + k - 1) (3) where C r = (- 1) r -1 B r /(2 r - 1)(2 r ), B r being the Bernoulli numbers [4]. These coefficients have been published by Uhler [5]. Requiring the range 0 ≦ x ≦ 1 is no restriction since, if necessary, Γ( x + 1) can be generated for other arguments using (1). For a given N , the error in (2) can be estimated from |ε| < | C N +1 |/ x 2 N +1 . (4) The curves of Figure 1 show contours of constant error bound as a function of N and x . These curves represent single and double-precision floating-arithmetic requirements of ε < 5·10 -9 and ε < 5·10 -17 . For a given N , k is defined as the minimum integral x greater than or equal to those on the curves. Then N and k can be chosen to minimize round-off and computing time. For N and k equal to 4, formula (3) yields Γ( x + 1) ∼ √2π exp ( x + 4 - 1/2) ln ( x + 4) - ( x + 4) exp ∑ 4 r =1 C r /( x + 4) 2 r -1 /( x + 1)( x + 2)( x + 3). (5) A similar expression suitable for double precision results for N = 8 and k = 9. The exponents in (5) are split to reduce roundoff. Various algebraic manipulations might result in a further reduction of roundoff.