Let G(x) = Gamma(x+1) / Gamma(x+1/2)
[S0] (x + 1/2)/sqrt(x + 1) < G(x) < sqrt(x + 1/2) ( x > 0 ) [S1] sqrt(x + 1/4) < G(x) < (x + 1/2)/sqrt(x + 3/4) ( x > 0 ) [S2] (x + 1/2)/sqrt(x + 3/4 + 1/(32x + 24)) < G(x) ( x > -1/2 ) [S3] G(x) < sqrt(x + 1/4 + 1/(32x + 8)) ( x > 0 )
Let B(f,x) = (x+1/4)^(1/2)*exp(f(x)) u(x) = (8x+2)^(-2) l(x) = (8x+2)^(-4)(64x^2+32x-6) w(x) = (8x+2)^(-6)(4096x^4+4096x^3+896x^2-64x+904/3) [LU] B(l,x) < G(x) < B(u,x) ( x > 0 ) [LW] B(l,x) < G(x) < B(w,x) ( x > 0 )
The following bounds were inspired by Windschitl's approximation of the Gamma function. [For the Windschitl formula see also the approximation page and the references given there.] [Q] Q(z) = sqrt(z*e) (1/2 + exp(-1/z)/2)^z [P] P(z) = sqrt(u/e) (1/2 + exp(-1/u)/2)^(-u) (u = z + 1/2) = (z + 1/2) / Q(z + 1/2) G(z) = Q(z)(1 + O(z^(-5)) and [QP] Q(x) < G(x) < P(x) ( x > 0 )
If you need simple bounds, choose [S1] sqrt(x + 1/4) < G(x) < (x + 1/2)/sqrt(x + 3/4) ( x > 0 ) If you need simple and sharp bounds, choose [S2S3] (x + 1/2)/sqrt(x + 3/4 + 1/(32x + 24)) < G(x) < sqrt(x + 1/4 + 1/(32x + 8)) ( x > 0 ) If you need sharper bounds, choose (with u = x + 1/2) [QP] sqrt(x*e)(1/2 + exp(-1/x)/2)^x < G(x) < sqrt(u/e)(1/2 + exp(-1/u)/2)^(-u) ( x > 0 ) If you need very sharp bounds, choose [LW] B(l,x) < G(x) < B(w,x) ( x > 0 )
The tables display the number of exact significant decimal digits of the formulas for some values of x. (The sign indicates 'lower' or 'upper' bound.) The bracketing bounds [S1]: x = 5, 10, 15, 20, 25, 30, 35, 40, 45, 50. ------------------------------------------------------------- -3.2, -3.8, -4.2, -4.4, -4.6, -4.8, -4.9, -5.0, -5.1, -5.2 3.3, 3.9, 4.2, 4.4, 4.6, 4.8, 4.9, 5.0, 5.1, 5.2 The bracketing bounds [S2, S3]: x = 5, 10, 15, 20, 25, 30, 35, 40, 45, 50. ------------------------------------------------------------- -5.7, -6.8, -7.4, -7.9, -8.3, -8.6, -8.9, -9.1, -9.3, -9.5 5.5, 6.7, 7.4, 7.9, 8.3, 8.6, 8.8, 9.1, 9.3, 9.5 The bracketing bounds [Q, P]: x = 5, 10, 15, 20, 25, 30, 35, 40, 45, 50. ----------------------------------------------------------------- -6.4, -7.9, -8.8, -9.4, -9.9, -10.3, -10.6, -10.9, -11.2, -11.4 6.6, 8.0, 8.9, 9.5, 9.9, 10.3, 10.7, 10.9, 11.2, 11.4 The pairwise bracketing bounds [B(u,x), B(l,x), B(w,x)]: x = 5, 10, 15, 20, 25, 30, 35, 40, 45, 50. ---------------------------------------------------------------------- 5.5, 6.7, 7.3, 7.8, 8.2, 8.5, 8.8, 9.0, 9.2, 9.4 -7.24, -8.97, -10.0, -10.7, -11.3, -11.8, -12.2, -12.5, -12.8, -13.1 8.67, 11.0, 12.3, 13.3, 14.1, 14.7, 15.3, 15.7, 16.1, 16.5
(1) Observe that G(x) < Sqrt(x+1/2) follows from G(x) < B(u,x) because exp(1/(2(2x+1)^2) < (2x+2)/(2x+1). Several proofs for this inequality were given in the thread "Eine Ungleichung mit Exp". (2) Note how S3(x) and S2(x) are related to each other: S2(x) = (x+1/2)/S3(x+1/2) By a theorem proved by Alexandru Lupas in the thread "Eine Ungleichung mit Exp" we know that S2(x) is a lower bound for G(x) as soon as we know that S3(x) is an upper bound for G(x). (3) If we expand exp(z) in [Q] up to O(z^(-6)) we get simplified forms of [Q] and [P] without losing numerical performance. / 1 1 1 1 5 1 21 1 3619 1 \ [Q*] sqrt(z) | 1 + - - + --- --- - ---- --- - ----- --- + -------- --- | \ 8 z 128 z^2 1024 z^3 32768 z^4 11796480 z^5 / / 1 1 1 1 5 1 21 1 32291 1 \ [P*] sqrt(z) | 1 + - - + --- --- - ---- --- - ----- --- + -------- --- | \ 8 z 128 z^2 1024 z^3 32768 z^4 11796480 z^5 / (4) Averaging [Q*] and [P*] and smoothing the last term we get the following simple and remarkable efficient approximation to G(z) [A] sqrt(z)(1+1/(8*z)*(1+1/(16*z)*(1-5/(8*z)*(1+21/(160*z)*(1-(50/(21*z)))))))
Say we want to approximate G(z) for z = 32 with 10 exact decimal digits. Then we use formula A(z). Say we want to compute for z = 32 a lower and upper bound for G(z) with 10 exact decimal digits. Then we choose formulas [Q*] and [P*] and proceed as follows: a(z) = (1+1/(8*z)*(1+1/(16*z)*(1-5/(8*z)*(1+21/(160*z))))) s(z) = sqrt(z) p(z) = 11796480*z^5 LB(z) = s(z)*(a(z) + 3619/p(z)) # low bound HB(z) = s(z)*(a(z) + 32291/p(z)) # high bound ErrEst = (HB(z)+LB(z))/(HB(z)-LB(z)) # error estimate
Back to the Homepage of Factorial Algorithms.