Approximations for Digital Computers
Approximations for Digital Computers
| Hi-index | 0.00 |
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, &Ggr;(x + 1) = x&Ggr;(x), (1) and Stirling's asymptotic expansion for ln &Ggr;(x) [3], we have ln &Ggr;(x) ∼ (x - 1/2) ln x - x + 1/2 ln 2&pgr; + ∑Nr=1 Cr/x2r-1. (2) It follows that, if k and N are appropriately selected positive integers, &Ggr;(x + 1) can be represented by &Ggr;(x + 1) ∼ √2&pgr; exp (x + k - 1/2) ln (x + k) - (x + k) exp ∑Nr=1 Cr/(x + k)2r-1/(x + 1)(x + 2) ··· (x + k - 1) (3) where Cr = (- 1)r-1 Br/(2r - 1)(2r), Br 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, &Ggr;(x + 1) can be generated for other arguments using (1). For a given N, the error in (2) can be estimated from |&egr;| CN+1|/x2N+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 &egr; -9 and &egr; -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 &Ggr;(x + 1) ∼ &radic2&pgr; exp (x + 4 - 1/2) ln (x + 4) - (x + 4) exp ∑4r=1Cr/(x + 4)2r-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.