﻿ A Ratio of Gamma Values

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
```

Discussions on sci.math and de.sci.mathematik.
* Alexandru Lupas: "G(x+1)/G(x+0.5),G=Gamma" [Aug 14, 2006]
* Peter Luschny: "Eine Ungleichung mit Exp"  [Aug 15, 2006]
* Peter Luschny: "Eine Inklusion mit Gamma"  [Aug 17, 2006]
* Peter Luschny: "Eine neue Naeherung an Gamma(x+1)/Gamma(x+1/2)" [Aug 20, 2006]

[S0] and [S1] were mentioned by Alexandru Lupas,
[S2] and B(u,x) by David W. Cantrell,
[S3], B(l,x) and B(w,x), Q(x) and P(x), and A(x) by Peter Luschny.

See also a bibliography on inequalities for the gamma function compiled by Alexandru Lupas.

---

Back to the Homepage of Factorial Algorithms.