18.

440, March 9, 2009

Stirlings formula
The factorial function n! is important in evaluating binomial, hypergeometric, and
other probabilities. If n is not too large, n! can be computed directly, by calculators
or computers. For larger n, using there are diculties with overow, as for example
70! > 10
100
, 254! > 10
500
, which overows on one calculator I have, which computes 253!.
Also, direct multiplication of many factors becomes inecient. There is a relation with the
gamma function, n! (n+1), where () =
_
+
0
x
1
e
x
dx. The statistical computing
system R (in the version we have as of this date) can nd 170! = (171)
.
= 7.2574 10
306
but it balks at (172), so it breaks down for smaller n than the calculator does. Of course,
some computer systems can nd n! for very large n. Mathematica gave 1000! exactly,
showing all the many digits, which is not necessarily convenient.
Stirlings formula provides an approximation to n! which is relatively easy to compute
and is sucient for most purposes. Using it, one can evaluate log n! to better and better
accuracy as n becomes large, provided that one can evaluate log n as accurately as needed.
Then to compute b(k, n, p) :=
_
n
k
_
p
k
q
nk
, for example, where 0 < p = 1 q < 1, one can
nd log b(k, n, p) = log n! log k! log (n k)! +k log p +(n k) log q . The probability
b(k, n, p) cannot overow, and in interesting cases it will also not underow (1/b(k, n, p)
will not overow).
Two sequences of numbers, a
n
and b
n
, are said to be asymptotic, written a
n
b
n
,
if lim
n
a
n
/b
n
= 1. This does not imply that lim
n
(a
n
b
n
) = 0: for example,
n
2
+ n n
2
but (n
2
+ n) n
2
tends to with n. But a
n
/b
n
1 is equivalent to
log (a
n
) log (b
n
) = log (a
n
/b
n
) 0.
Theorem 1. Stirlings formula.n!
n
n
e
n

2n = n
(n+1/2)
e
n

2 . Thus,
log (n!)
__
n +
1
2
_
log n n +
1
2
log (2)
_
0 as n .
Proof. The sign := will mean equals by denition. Let
d
n
:= log (n!)
_
n +
1
2
_
log n + n .
Then we need to prove d
n
converges to a constant, [log (2)]/2. First,
d
n
d
n+1
= log (n + 1)
_
n +
1
2
_
log n +
_
n +
3
2
_
log (n + 1) 1
=
_
n +
1
2
_
log
_
n + 1
n
_
1 .
We have the Taylor series
log (1 + t) = t
t
2
2
+
t
3
3

1
for |t| < 1. For t > 0 the terms alternate in sign. A transformation will help to get terms
of the same sign. The trick is to notice that
n + 1
n
=
1 +
1
2n+1
1
1
2n+1
.
Then
log
_
1 + t
1 t
_
= log (1 + t) log (1 t) = 2
_
t +
t
3
3
+
t
5
5
+
_
,
where now all terms are of the same sign. Thus
d
n
d
n+1
=
2n + 1
2
log
_
1 +
1
2n+1
1
1
2n+1
_
1
=
1
3(2n + 1)
2
+
1
5(2n + 1)
4
+
1
7(2n + 1)
6
+ > 0 . (1)
So d
n
decreases as n decreases. Comparing the last series to a geometric one with ratio
(2n + 1)
2
gives
0 < d
n
d
n+1
<
(2n + 1)
2
3[1 (2n + 1)
2
]
=
1
3[(2n + 1)
2
1]
=
1
12n(n + 1)
=
1
12n

1
12(n + 1)
, so d
n

1
12n
< d
n+1

1
12(n + 1)
.
So we see that d
n
1/(12n) increases as n does. As n , d
n
decreases to some C with
C < + and d
n
1/(12n) increases up to some D with < D +. Since
1/(12n) converges to 0, we must have < C = D < +, and d
n
converges to a nite
limit C. By denition of d
n
we then have
n!/(n
n+1/2
)e
n
e
C
or n! e
C
n
n+1/2
e
n
.
The last step in the proof is to show that e
C
= (2)
1/2
. This will involve another
famous fact:
Theorem 2. Wallis product.

2
=
2
1

2
3

4
3

4
5

6
5

6
7

2m
2m1

2m
2m + 1
, or

2
= lim
m
2
4m
(m!)
4
(2m)!(2m+ 1)!
.
Remarks. To see the relationship between the two statements, rst note that 2 4 6 8
2m = (2 1)(2 2)(2 3) (2 m) = 2
m
m!, then that 1 3 5 7 (2m + 1) =
(2m + 1)!/(2 4 6 2m), etc. Note that the product converges to /2 rather slowly;
it would not give a good way to compute .
2
Proof. Integrating by parts gives, for n 2,
_
sin
n
xdx =
_
sin
n1
xd(cos x)
= cos x sin
n1
x + (n 1)
_
sin
n2
x cos
2
xdx
= cos x sin
n1
x + (n 1)
_
(sin
n2
x sin
n
x)dx , so
n
_
sin
n
xdx = cos x sin
n1
x + (n 1)
_
sin
n2
xdx , and
_
sin
n
xdx =
cos x sin
n1
x
n
+
n 1
n
_
sin
n2
xdx . Thus
(2)
_
/2
0
sin
n
xdx =
n 1
n
_
/2
0
sin
n2
xdx .
Then for m = 1, 2, . . . , iterating (2) gives
_
/2
0
sin
2m
xdx =
2m1
2m

2m3
2m2

1
2


2
since
_
/2
0
1dx =

2
.
_
/2
0
sin
2m+1
xdx =
2m
2m + 1

2m2
2m1

2
3
1 since
_
/2
0
sin xdx = 1 .
Let A
m
:=
_
/2
0
sin
2m
xdx/
_
/2
0
sin
2m+1
xdx. Then

2
= A
m

2
1

2
3

4
3

4
5

6
5

6
7

2m
2m1

2m
2m + 1
for all m = 1, 2, . . . .
Now we will prove lim
m
A
m
= 1. For 0 x /2,
0 sin
2m+1
x sin
2m
x sin
2m1
x, so
0 <
_
/2
0
sin
2m+1
xdx <
_
/2
0
sin
2m
xdx <
_
/2
0
sin
2m1
xdx .
Now by (2) above,
_
/2
0
sin
2m+1
xdx/
_
/2
0
sin
2m1
xdx =
2m
2m + 1
1 as m
and
_
/2
0
sin
2m
xdx, being between numerator and denominator, also has the ratio A
m
converging to 1, proving Wallis product.
3
Now to nish proving Stirlings formula, let B := e
C
. As n , n!e
n
/n
n+1/2

B, (2n)!e
2n
/(2n)
2n+1/2
B, and (n!)
2
e
2n
/n
2n+1
B
2
. Dividing gives
(n!)
2
2
2n+1/2
/[(2n)!n
1/2
] B. Now, Wallis product gives (n!)
2
2
2n
/[(2n)!(2n + 1)
1/2
]
(/2)
1/2
. Since (2n + 1)
1/2
(2n)
1/2
, we get b/2
1/2
= 2
1/2
(/2)
1/2
, B = (2)
1/2
,
proving Stirlings formula.
The proof provides further information on how good an approximation Stirlings for-
mula gives to n!. Since d
n
> C > d
n
1/(12n), where C = [log (2)]/2, so C < d
n
<
C + 1/(12n), we have the bounds
(3) (2)
1/2
n
n+1/2
e
n
< n! < (2)
1/2
n
n+1/2
e
n+[1/(12n)]
.
Even closer bounds are available. From (1),
d
n
d
n+1

j=1
1
3
j
(2n + 1)
2j
>
_
1
5

1
9
_
1
(2n + 1)
4
, so
d
n
d
n+1
>
3
1
(2n + 1)
2
1 3
1
(2n + 1)
2
+
4
45
(2n + 1)
4
=
1
3(2n + 1)
2
1
+
16
180(2n + 1)
4
=
1
12n
2
+ 12n + 2
+
16
180(4n
2
+ 4n + 1)
2
=
1
12n(n + 1)
_
1 +
1
6n(n + 1)
_
1
+
1
180n
2
(n + 1)
2
_
1 +
1
4n(n + 1)
_
2
>
1
12n(n + 1)
_
1
1
6n(n + 1)
_
+
1
180n
2
(n + 1)
2
_
1
1
2n(n + 1)
_
>
1
12n(n + 1)

3n(n + 1) + 1
360n
3
(n + 1)
3
(since
1
180

1
72
=
3
360
)
=
1
12
_
1
n

1
n + 1
_

1
360
_
1
n
3

1
(n + 1)
3
_
. So,
d
n

1
12n
+
1
360n
3
> d
n+1

1
12(n + 1)
+
1
360(n + 1)
3
,
and the sequence d
n
1/(12n) + 1/(360n
3
) decreases as n down to its limit, which
is also C, so d
n
1/(12n) + 1/(360n
3
) > C. Writing exp(x) := e
x
, we have the following
improvement on the left side of (3): for all n = 1, 2, . . . ,
(4)

2n
n+1/2
exp
_
n +
1
12n

1
360n
3
_
< n! <

2n
n+1/2
exp
_
n +
1
12n
_
.
As n , the ratio of the upper to lower bound converges to 1 rather fast since
1/(360n
3
) 0 rather fast.
4
There are further improvements, although they wont be proved here: Whittaker and
Watson, Modern Analysis, p. 252, gives an asymptotic expansion
d
n
C
1
12n

1
360n
3
+
1
1260n
5

1
1680n
7
+
1
1188n
9
.
The series does not converge for any n, but if the sum of the rst k terms is used as an
approximation to the left side d
n
C, the error in the approximation has the same sign
as, and smaller absolute value than, the next ((k + 1)st) term. This was proved above for
k = 0 by (3) and for k = 1 by (4).
Now Stirlings formula with error bounds can be used to give upper and lower bounds
for (n)
k
:= n(n 1) (n k + 1) = n!/(n k)! for integers 0 k n. Specically, (4)
implies
(n)
k
<
n
n+1/2
exp
_
n +
1
12n
+
1
360(nk)
3
_
(n k)
nk+1/2
exp
_
(n k) +
1
12(nk)
_ and
(n)
k
>
n
n+1/2
exp
_
n +
1
12n

1
360(nk)
3
_
(n k)
nk+1/2
exp
_
(n k) +
1
12(nk)
_ .
Let j(n, k) :=
n
n+1/2
(n k)
nk+1/2
exp
_
k
1
12n(n k)
_
. The above inequalities on (n)
k
show that it is approached by j(n, k) within a factor of exp[1/(360(n k)
3
)], which is
very close to 1 if n k is large. For n k large, exp[k/(12n(n k))] also approaches 1,
although not as fast.
Let p(n, k) := (n)
k
/n
k
, the probability that k numbers, chosen at random from
1, . . . , n with replacement, are all dierent. Then, to the accuracy of the above approxi-
mation for (n)
k
, p(n, k) is approximated by
e
k
_
n
n k
_
nk+1/2
exp
_

k
12n(n k)
_
.
For a simpler and rougher approximation, omit the exp. . . factor.
Now, suppose that for a given n and , with 0 < < 1, we want to nd the smallest
k such that p(n, k) < . For example, if n = 365 and = 1/2, the question is how
many people are needed to give an even chance that at least two of them have the same
birthday (neglecting leap years and assuming that births are evenly distributed throughout
the year).
To nd the desired k, one can compute p(n, k) and use trial and error. To speed up
the process one can use a simpler approximation where we can solve for k to get a good
rst approximation to k. Then most likely only a few values of k near the rst one need to
be tried. Here is how one can get such a simple approximation. For 0 k < n, the Taylor
series of log (1 x) gives
log
_
1
k
n
_
=
k
n

k
2
2n
2

k
3
3n
3
.
5
If k/n is small, later terms in the series can be neglected, and log p(n, k) is approximated
by
(5) log p(n, k) k
_
n k +
1
2
_

k
2
2n

k
3
6n
2
+
k
2n
+ ,
where the next largest terms would be of the order of k
4
/n
3
and k
2
/n
2
(and k/[12n(nk)]
is still smaller). Note that if we approximated log (1 k/n) bu just the rst term k/n,
we would not even get the rst term in (5) correct (the 2 in the denominator would be
missing). Using the rst term k
2
/(2n) in (5) as our rst approximation, solving for k
gives k
2
/(2n) log , or
(6) k [2nlog (1 )]
1/2
.
For such a k, the next two terms in the approximation are smaller by factors of the order
of 1/n
1/2
, so they can be reasonably be neglected if n is large. This gives a
Method. To nd the least k such that p(n, k) < , for given n and , rst try k as the
next larger integer than the number from (6). Compute p(n, k). If p(n, k) < , check that
p(n, k 1) . If not, consider k 2, etc. until a solution is found. If p(n, k) > , nd
whether p(n, k + 1) < . If so, the solution is k + 1. If not, try k + 2, . . . , until a solution
is found.
Example 1. The birthday problem. Here n = 365 and = 1/2. First try k as the
next integer larger than (2nlog 2)
1/2
, that is k = 23. Then we nd p(365, 23) < 1/2, so
we next compute p(365, 22) and nd it is larger than 1/2, so k = 23 is the solution: in a
group of 23 or more people, there is a better than even chance that at least two have the
same birthday.
Example 2. A computer pseudo-random number generator starts with a number s
called a seed and uses a function f to generate numbers s
1
= s, s
2
= f(s
1
), s
3
=
f(s
2
), . . . , s
j+1
= f(s
j
), l . . . . Suppose that the numbers s
j
will be integers from a to
n for some n, and f is a randomly chosen function from the set {1, 2, . . . , n} into itself,
where each of the n
n
such functions is equally likely. For how large r will there be an
even chance that s
r
= s
m
for some m < r? Once this happens, then s
r+1
= s
m+1
, etc.
and the s
i
will go round and round a closed cycle. So the event that s
r
= s
m
for some
m < r is the event that the s
j
for j r are not all dierent. The above method applies
with = 1/2. If n = 10
6
, for example, (6) gives r = 1178 and it can be checked that
p(n, 1178) < 1/2 < p(n, 1177). So in this case there is an even chance that the generator
will fall into a closed cycle after only 1178 of the 1,000,000 available numbers. By the way,
the average length of the closed cycle is just half of the rst number r such that s
r
= s
m
for some m < r.
So there is a paradox: a truly random function f makes a bad pseudo-random number
generator. Better generators are made by using number-theoretic methods to assure that
there are no short closed cycles.
Bibliographic Notes. James Stirling published his formula in Methodus Dierentialis
(1730). Abraham De Moivre, another mathematician and friend of Stirlings, discovered
6
the formula except for nding the value of the constant factor (2)
1/2
. The proof of the
formula and up through (3) above is due to Herbert Robbins, Amer. Math. Monthly 62
(1955) pp. 2629. The renement of the proof to give (4) is due to T. S. Nanjundiah, ibid.
66 (1959) pp. 701703. As mentioned, further terms in the asymptotic expansion (next
display after (4)) can be found from E. T. Whittaker and G. N. Watson, Modern Analysis
(Cambridge Univ. Press, 4th ed., 1927, repr. 1962) pp. 252-253. John Wallis published
his product (without a real proof) around 1650 (see his Opera Omnis, re-published in
1972). The above proof came from R. Courant, Dierential and Integral Calculus I, 2d.
ed., translated by E. J. McShane (Interscience, N. Y. , 1937).
Stirlings formula examples
Let S(n) = (n/e)
n
(2n)
1/2
. Then n! S(n) as n , meaning n!/S(n) 1, and
n!/[S(n)e
1/(12n)
] 1 faster. But n! S(n) does not converge to 0; in fact it increases very
fast, but not as fast as n! or S(n).
n n! S(n) n! S(n) n!/S(n) n!/[S(n)e
1/(12n)
]
5 1.2000 10
2
1.1802 10
2
1.9808 1.0168 .999978024
10 3.6288 10
6
3.5987 10
7
3.0104 10
4
1.0084 .999997299
20 2.4329 10
18
2.4228 10
18
1.0115 10
16
1.0042 .999999649
40 8.1592 10
47
8.1422 10
47
1.6980 10
45
1.0021 .999999948
60 8.3210 10
81
8.3094 10
81
1.1549 10
79
1.0014 .999999988
Note that the ratios in the next to last column decrease toward 1. They are approx-
imately 1 + 1/(12n). The ratios in the last column increase toward 1, faster. They are
approximately 1 1/(360n
3
). So as n becomes large, in terms of ratio (not dierence), n!
is fairly well approximated by S(n), much better approximated by S(n)e
1/(12n)
, and still
much better approximated by S(n) exp[1/(12n) 1/(360n
3
)].
For large n, one needs to take account of rounding error. In log(n
n+0.5
) = (n +
0.5) log n, a rounding error in log n is multiplied by n. If n is 10
k
, for example, this means
a loss of k decimal places of accuracy.
7