On Ramanujan's approximation, $n!\sim \sqrt{\pi}\big(\frac ne\big)^n\sqrt [6]{(2n)^3+(2n)^2+n+\frac 1{30}}$

This is question I could have been asking for years. I wonder if the basis of this magnificent approximation has been published anywhere.

What I have seen (the problem is that I do not remember where) are approximations built as $$n!\sim \sqrt{\pi}\left(\frac ne\right)^n\sqrt [2k]{P_k(n)}$$ where remains the mystery of the $\sqrt [2k]{.}$. The first case I saw is Gosper approximation.

The coefficients of the polynomials were obtained taking the logarithms of both sides and identified with the Stirling series. So, what was obtained are $$P_1(n)=2n+12$$ $$P_2(n)=4n^2+12n+18$$ $$\color{red}{P_3(n)=8 n^3+4 n^2+n+\frac{1}{30}}$$ and for sure, with the tools we have today, we could continue for ever $$P_4(n)=16 n^4+\frac{32 n^3}{3}+\frac{32 n^2}{9}+\frac{176 n}{405}-\frac{128}{1215}$$ $$P_5(n)=32 n^5+\frac{80 n^4}{3}+\frac{100 n^3}{9}+\frac{178 n^2}{81}-\frac{95 n}{972}+\frac{2143}{40824}$$ $$P_6(n)=64 n^6+64 n^5+32 n^4+\frac{128 n^3}{15}+\frac{8 n^2}{15}+\frac{8 n}{105}+\frac{596}{1575}$$ for more and more accuracy.

For example, for $n=5$, the magic formula given by the great Ramanujan gives $120.000147066$ while the last given here leads to $120.000000406$.

There is one thing interesting to notice : up to $k=3$ the coefficients of powers of $n >0$ are all integer numbers.

Edit

Staying with the power $\frac 16$ from Ramanujan, we could extend it as $$P_3(n)=8 n^3+4 n^2+n+\frac{1}{30}\left(1+\sum_{i=1}^\infty \frac {a_i}{n^i}\right)$$ and the sequence of the first $a_i$'s is $$\left\{-\frac{11}{8},\frac{79}{112},\frac{3539}{6720},-\frac{9511}{13440},-\frac{30 153}{71680},\frac{233934691}{212889600},\frac{3595113569}{5960908800},\cdots\right\}$$

The asymptotic expansion of the factorial may be written $$ n! \sim \left( {\frac{n}{e}} \right)^n \sqrt {2\pi n} \sum\limits_{k = 0}^\infty {( - 1)^k \frac{{\gamma _k }}{{n^k }}} , $$ where $\gamma_k$ denotes the Stirling coefficients (see, e.g., http://dx.doi.org/10.1017/S0308210513001558). The first three of them are $\gamma_0=1$, $\gamma_1=-\frac{1}{12}$ and $\gamma_2=\frac{1}{288}$. Now, if a function has an asymptotic expansion, then its 6th power also has one and it is the original expansion raised into the power of $6$. Accordingly, $$ \left( {n!/\left( {\frac{n}{e}} \right)^n \sqrt {2\pi n} } \right)^6 \sim \left( {\sum\limits_{k = 0}^\infty {( - 1)^k \frac{{\gamma _k }}{{n^k }}} } \right)^6 \sim 1 + \frac{1}{{2n}} + \frac{1}{{8n^2 }} + \frac{1}{{240n^3 }} - \frac{{11}}{{1920n^4 }} + \cdots , $$ and thus \begin{align*} n! & \sim \left( {\frac{n}{e}} \right)^n \sqrt {2\pi n} \sqrt[6]{{1 + \frac{1}{{2n}} + \frac{1}{{8n^2 }} + \frac{1}{{240n^3 }} - \frac{{11}}{{1920n^4 }} + \cdots }} \\ & = \left( {\frac{n}{e}} \right)^n \sqrt \pi \sqrt[6]{{8n^3 + 4n^2 + n + \frac{1}{{30}} - \frac{{11}}{{240n}} + \cdots }}\; . \end{align*} Therefore, Ramanujan's approximation is just a manipulation of the standard asymptotic expansion of the factorial. It is of course more accurate than the leading order asymptotics (known as Stirling's formula) since it uses additional terms from the asymptotic expansion. Let me add here that although the standard asymptotic expansion (and Ramanujan's as well) is divergent, for large $n$ the initial terms decrease in magnitude. There is a minimum occuring around $\left\lfloor {2\pi n} \right\rfloor$. Truncating the series at this point yields exponentially accurate approximations (see again the paper above for more details).