How is Gosper's approximation to factorial derived?

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}.$$