Higher Order Terms in Stirling's Approximation

The Euler-MacLaurin summation formula can be used to approximate $n!$ conveniently.

The following is taken from Concrete Mathematics by D.E. Knuth ad R.L. Graham. We find in section 9.5 among other representations this version of the Euler-MacLaurin summation formula for integer values $a\leq b$

\begin{align*} \sum_{a\leq k < b}f(k)&=\int_a^bf(x)dx-\frac{1}{2}\left.f(x)\right|_a^{b}+\sum_{k=1}^{m}\frac{B_{2k}}{(2k)!}\left.f^{(2k-1)}(x)\right|_a^b\tag{1}\\ &\qquad+ \mathcal{O}\left((2\pi)^{-2m}\right)\int_a^b\left|f^{(2m)}(x)\right|dx \end{align*}

with $B_k$ the Bernoulli numbers starting with \begin{align*} B_0=1,B_1=-\frac{1}{2},B_2=\frac{1}{6},B_4=-\frac{1}{30},B_6=\frac{1}{42},\ldots \end{align*} and $B_k=0$ for odd $k$ greater $1$.

It is also shown, that if the function $f$ has the nice property $$f^{2m}(x)\geq 0$$ for $a\leq x \leq b$ then it turns out that the remainder in (1) can be replaced with \begin{align*} \Theta_m\frac{B_{2m+2}}{(2m+2)!}\left.f^{(2m+1)}(x)\right|_a^b\qquad\qquad \text{for some }0<\Theta<1 \end{align*}

and we obtain \begin{align*} \sum_{a\leq k < b}f(k)&=\int_a^bf(x)dx-\frac{1}{2}\left.f(x)\right|_a^{b}+\sum_{k=1}^{m}\frac{B_{2k}}{(2k)!}\left. f^{(2k-1)}(x)\right|_a^b\\ &\qquad+\Theta_m\frac{B_{2m+2}}{(2m+2)!}\left.f^{(2m+1)}(x)\right|_a^b\tag{2} \end{align*}

---

Example: $f(x)=\ln x$

In order to find the Stirling formula for $n!$ we consider $f(x)=\ln x$ and derive following approximation for $(n-1)!$

\begin{align*} \sum_{1\leq k < n}\ln k&=n\ln n - n +\sigma-\frac{1}{2}\ln n+\sum_{k=1}^m\frac{B_{2k}}{2k(2k-1)n^{2k-1}}\\ &\qquad+\Theta\frac{B_{2m+1}}{(2m+2)(2m+1)n^{2m+1}}\qquad\qquad 0<\Theta<1 \end{align*}

Another instructive example in section 9.6, namely the asymptotic estimation for $\sum_k\binom{2n}{k}$ provides a formula for $\sigma$ \begin{align*} e^\sigma=\sqrt{2\pi} \end{align*}

Now we have everything we need to estimate $n!$ as precise as we like:

Let's assume we want to approximate $n!$ up to terms of order $n^{-3}$. We put $m=2$ and observe

\begin{align*} \sum_{1\leq k < n}\ln k=n\ln n-n+\sqrt{2\pi}-\frac{1}{2}\ln n+\frac{1}{12n}-\frac{1}{360n^3} +\mathcal{O}\left(\frac{1}{n^5}\right)\tag{3} \end{align*}

Using the generating function of $\exp(x)=\sum_{n\geq 0}\frac{x^n}{n!}$ we obtain from (3) \begin{align*} n!&=\exp\left(n\ln n-n+\sqrt{2\pi}+\frac{1}{2}\ln n+\frac{1}{12n}-\frac{1}{360n^3} +\mathcal{O}\left(\frac{1}{n^5}\right)\right)\tag{4}\\ &=\sqrt{2\pi n}\left(\frac{n}{e}\right)^n \exp\left(\frac{1}{12n}-\frac{1}{360n^3}+\mathcal{O}\left(\frac{1}{n^5}\right)\right)\\ &=\sqrt{2\pi n}\left(\frac{n}{e}\right)^n \left(1+\left(\frac{1}{12n}-\frac{1}{360n^3}\right)+\frac{1}{2}\left(\frac{1}{12n}-\frac{1}{360n^3}\right)^2\right.\\ &\qquad\qquad\qquad\qquad\left.+\frac{1}{6}\left(\frac{1}{12n}-\frac{1}{360n^3}\right)^3+\mathcal{O}\left(\frac{1}{n^4}\right)\right)\\ &=\sqrt{2\pi n}\left(\frac{n}{e}\right)^n \left(1+\left(\frac{1}{12n}-\frac{1}{360n^3}\right)+\frac{1}{2}\left(\frac{1}{144n^2}\right)\right.\\ &\qquad\qquad\qquad\qquad\left.+\frac{1}{6}\left(\frac{1}{1728n^3}\right)+\mathcal{O}\left(\frac{1}{n^4}\right)\right) \end{align*}

Observe the change from $-\frac{1}{2}\ln n$ to $+\frac{1}{2}\ln n$ from (3) to (4) due to the change from $(n-1)!$ to $n!$.

We finally conclude \begin{align*} n!=\sqrt{2\pi n}\left(\frac{n}{e}\right)^n \left(1+\frac{1}{12n}+\frac{1}{288n^2}-\frac{139}{51840n^3}+\mathcal{O}\left(\frac{1}{n^4}\right)\right) \end{align*}