How is Gosper's approximation to factorial derived?
Discussing Stirling's approximation, Wolfram Mathworld article mentions a modification of it due to Gosper: instead of the usual
$$n!\approx n^ne^{-n}\sqrt{2n\pi},$$
we have a tiny addition of $\frac13$ to $2n$ under the radical sign:
$$n!\approx n^ne^{-n}\sqrt{\left(2n+\frac13\right)\pi}.$$
This radically improves the precision of approximation for all $n\ge0$. But all the places I've seen discussing it don't explain how it was obtained. Wolfram Mathworld just says about Gosper's modification
...a better approximation to n! (i.e., one which approximates the terms in Stirling's series instead of truncating them)...
but it doesn't make me understand how this change was derived.
So, how to derive this modification? And is there a further improvement to the whole Stirling series, or is it just "whole series stuffed into one expression"?
As it is said it the Wolfram article, Gosper's formula approximates the Stirling series instead of truncating it.
To see that, let's take a look at the $\sqrt{2n+\frac{1}{3}}$ term which is itself a series:
$$ \sqrt{2n+\frac{1}{3}} = \sqrt{2n} \cdot\sqrt{1+\frac{1}{6n}} = \sqrt{2n} \cdot \left( 1 + \frac{1}{12n} + O \left( \frac{1}{n^2}\right) \right) $$
So in the in the end you have:
$$ \frac{n!}{n^n e^{-n}\sqrt{2 n \pi}} = 1 + O \left( \frac{1}{n}\right) $$
While:
$$ \begin{align}\frac{n!}{n^n e^{-n}\sqrt{ \left( 2 n + \frac{1}{3} \right) \pi}} & = \frac{n!}{n^n e^{-n}\sqrt{2 n \pi}} \cdot\frac{1}{\sqrt{1+\frac{1}{6n}}} \\ & = \left( 1 + \frac{1}{12n} + O \left( \frac{1}{n^2}\right) \right) \cdot \left( 1 - \frac{1}{12n} + O \left( \frac{1}{n^2}\right) \right) \\ & = 1 + O \left( \frac{1}{n^2}\right) \end{align} $$
That explains why this approximation is much better than the simple truncation.
Stirling series is
$$n!\sim n^ne^{-n}\sqrt{2n\pi}\left(1+\frac1{12n}+\cdots\right).$$
Taking only terms up to and including $\frac1{12n}$, we have:
$$\sqrt{2n}\left(1+\frac1{12n}\right)=\sqrt{2\left(n+\frac16+\frac1{12^2n}\right)}=\sqrt{2n+\frac13+\frac1{72n}}.$$
Neglecting the $\frac1{72n}$ term in the radical (since it vanishes as $n\to\infty$), we recover the result for the radical by Gosper:
$$\sqrt{2n+\frac13}.$$