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*}
Best Answer
Hint: read Section 1.2 of this text.