An ultimate extremely accurate formula for approximation of the factorial function (Q1039939)

Beginning with the identity \[ n! =\left(\frac ne\right)^n \sqrt {2\pi n} \exp\left(\sum_{k=n}^\infty\sum_{j=2}^\infty \frac{j - 1}{2 j (j + 1)} (- \frac1k)^j\right), \] the author determines a series of positive terms, \(\sum_{n=1}^\infty a_n\), such that \[ n! = \left(\frac ne\right)^n \sqrt {2\pi n} \exp \left(\sum_{k=n}^\infty a_k\right), \] from which follows the identity \[ n! = \sqrt {\frac{2\pi} {e}} \left(\frac {n+1}{e}\right)^{n + \frac 12} \exp\left(\sum_{k = n + 1}^\infty a_k \right), \] the estimate \[ n! \approx \sqrt {\frac {2\pi}{e}}\left(\frac {n+1}{e}\right)^{ n + \frac 12} = : \alpha_n \] and the inequalities \[ \left(\frac ne\right)^n \sqrt {2\pi n} < \sqrt {\frac {2\pi} {e}} \left(\frac{n+1}{e}\right)^{n+ \frac12} < n! \] A result established by \textit{W. Burnside} [Messenger 46, 157--159 (1917; JFM 46.0340.02)] and subsequently rediscovered by \textit{Y. Weissman} [Am. J. Phys. 51, No. 1, 9 (1983), \url{http://dx.doi.org/10.1119/1.13412}], has \[ n! \approx \sqrt{2\pi} \left(\frac{ n+ {\frac 12}}{e}\right)^{n+\frac 12} =:\beta_n. \] Both \(\alpha_n\) and \(\beta_n\) are estimates of the form \(n! \approx \lambda(\frac{n + p}{e})^{n + q} = : \tau_n\) where \(\lambda\), \(p\) and \(q\) are constants, and this has led the author to determine values for these constants which, in the presence of the condition \(\lim_n\frac{\tau_n}{n!}= 1\), yield superior estimates. The author's analysis of this problem finally leads to the inequalities, valid for all nonnegative \(x\), \[ \sqrt{2\pi e}\;e^{-\omega}\left(\frac{x+\omega}{e}\right)^{x+\frac12}<\Gamma(x+1)\leq \alpha\sqrt{2\pi e}\;e^{-\omega}\left(\frac{x+\omega}{e}\right)^{x+\frac12} \] and \[ \beta\sqrt{2\pi e}\;e^{-\zeta}\left(\frac{x+\zeta}{e}\right)^{x+\frac12}\leq\Gamma(x+1)<\sqrt{2\pi e}\;e^{-\zeta}\left(\frac{x+\zeta}{e}\right)^{x+\frac12}, \] where \(\omega=\frac{3-\sqrt3}{6}\), \(\alpha=\frac{e^\omega}{\sqrt{2\pi e}(\frac{\omega}{e})^{\frac12}} = 1.072042464\dots\), \(\zeta=\frac{3+\sqrt3}{6}\), and \(\beta=\frac{e^\zeta}{\sqrt{2\pi e}(\frac{\zeta}{e})^{\frac12}} = 0.988503589\dots\).