I recall my teacher mentioned that you can only use the Stirling formula with Maxwell-Boltzmann statistics, not in FD or BE statistics, because of large numbers. How can you then approximate $\ln(x!)$ with FD or BE?
Example approximation with Bose-Einstein statistics
The probability function for Bose-Einstein with $N_i$ particles with $M-1$ walls is
$$P=\Pi_i \frac{(N_i+M-1)!}{N_i!(M-1)!}$$
Now my lecture slides for the physics 2062 in Aalto University claims (p.146 here)
$$\ln P \approx \sum_i \left[\left(N_i+M-1\right)\ln\left(N_i+M-1\right)-(N_i+M-1)-N_i \ln(N_i)+N_i\right]-\ln\left(\left(M-1\right)!\right)$$
where the premise is $\ln(n!)\approx n \ln(n)-n$ but a small err in the constant term, more in the comment.
Now Stirling formula is not precise with small amount of particles but the teacher still uses it with both fermions and bosons. I am uneasy about this: more accurate approximation with smaller systems of particles?
You can derive better Stirling-like approximations of the form $$n! = \sqrt{2 \pi n} \left(\dfrac{n}{e} \right)^n \left(1 + \dfrac{a_1}n + \dfrac{a_2}{n^2} + \dfrac{a_3}{n^3} + \cdots \right)$$ using Abel summation technique (For instance, see here), where $$a_1 = \dfrac1{12}, a_2 = \dfrac1{288}, a_3 = -\dfrac{139}{51740}, a_4 = - \dfrac{571}{2488320}, \ldots$$ The hard part in Stirling's formula is deriving the $\sqrt{2 \pi}$ leading coefficient. The remaining coefficients $a_k$'s are not hard to obtain.
Posting this as an answer since it was a bit too long for a comment.