Approximations and Bounds for the Ratio
Gamma(x+1) / Gamma(x+1/2)

 
   Let G(x) = Gamma(x+1) / Gamma(x+1/2)

Simple Bounds (based on sqrt)

  
  [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 )

Sharp Bounds (based on exp)

   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 )

New Bounds? (Aug 20, 2006)

   
   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 ) 
   

Which bounds to choose?

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 )

Numerics: What perfomance can we expect?

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

Notes

(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)))))))

Numerical Recipe

   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   

Links

---

Back to the Homepage of Factorial Algorithms.

~ FIN ~