How can we justify, in deriving quantum statistics, the use of Stirling approximation in the form $\ln(x!)\approx x \ln x - x$? At first sight one can say "why not to use only one term, or maybe three or more terms"? Why use two terms? I see that books (see for example good books like Griffiths quantum mechanics or Atkins physical chemistry) use the $\ln(x!)\approx x \ln x - x$ formula, but no one justify this choice. I did calculus using the smaller formula $\ln(x!)\approx x \ln x$ finding, for identical fermions and bosons, the same results that books find using the longer formula (i.e. $N_n = \frac{d_n}{e^{\alpha+\beta e_n} \pm 1}$ where $N_n$ is the number of particles in energy level $E_n$, with degeneracy $d_n$, while $\alpha$ and $\beta$ are the Lagrange multipliers). 
The only reason I can speculate to choose $\ln(x!)\approx x \ln x - x$ instead of $\ln(x!)\approx x \ln x$ is that the approximation of the simpler formula is so bad that, despite it works for $x\to\infty$, it doesn't work in the physical context we are playing with. But if this is the reason, the choice should be accompanied by some considerations. A look at this plot suggest that this line of thinking is correct: even for a huge number like Avogadro's one, the difference with the two formulas is about 2\% (the ratio became bigger of 0.99 only for $2.688117141816367\cdot10^{43}$)

But these considerations are not sufficient to justify the use $\ln(x!)\approx x \ln x - x$. I see that $\ln(x!)\approx x \ln x$ is not appropriate for the proof of statistical quantum mechanics, but how can I see that $\ln(x!)\approx x \ln x - x$ is appropriate? The only way seems to be a comparison with the function $\ln(x!)$, but here problems arise.


*

*First. Surely we are interested to a number big and much smaller than Avogadro's number, but this sounds vague... what number should we elect as representative of the typical population of a typical energy level, to test the proper functioning of the approximated formula?

*Second. When the argument of $\ln(x!)$ starts becoming big, it put a strain even the powerful computing capabilities of computers


I did this plot

that suggests that this line of thinking too could be correct (for $x=170$, the maximum value hold by the junk of my computer before crashing, the error is about $0,5\%$), but the plot is not decisive and I'd like something stronger, solving satisfactorily the two problems I listed above. A way to proceed could be finding a way to maximize the error done by using $x\ln x - x$ instead of $\ln(x!)$ and then show that this error is not significative for every reasonable value of $x$, but I can't do it.
Adjunct
I'm asked to write down the proof, I'll report the final: it is possible to show under reasonable assumption that the quantum statistics can be found by finding the maximum (with constraints of conservation of particles and energy) of $Q$, where $Q$ for fermions and bosons is respectively
$$
Q_{ferm} = \prod_{n} \frac{d_n!}{N_n!(d_n-N_n)!}
$$
$$
Q_{bos} = \prod_{n} \frac{(N_n+d_n-1)!}{N_n!(d_n-1)!}
$$
and exploiting the hypothesis $1 \ll N_n \ll d_n$. Taking logarithm does not alter the positions of maxima (this step transforms products in sums and will allow us to exploit Stirling). Exploiting Lagrange method we have
$$
\frac{\partial}{\partial N_n}
\left[
\ln Q + \alpha \left( N - \sum_n N_n \right) + \beta  \left( E_n - \sum_n N_n E_n \right) 
\right]=0
$$
Going on we get to the crucial step, connected with this question: for fermions and bosons we have
$$
\frac{\partial }{\partial N_n} \sum_n 
[ \ln ( d_n! ) - \ln ( N_n ! ) - \ln(d_N - N_n)! ] = \alpha + \beta E_n
$$
$$
\frac{\partial }{\partial N_n} \sum_n 
[ \ln (N_n + d_n - 1)! - \ln ( N_n ! ) - \ln (d_n - 1)! ] = \alpha + \beta E_n
$$
from wich you can see that using indifferently $\ln(x!)\approx x \ln x$ or $\ln(x!)\approx x \ln x - x$ you will find $N_n = \frac{d_n}{e^{\alpha+\beta e_n} \pm 1}$. Anyway this is not the real problem that trouble me, the problem is "how can I see that taking two terms in $\ln(x!) \approx \dots$ is reasonably safe?".
 A: The essential point is: it depends on the necessary accuracy. However generally in the thermodynamic limit we can ignore terms subleading $N$ (see edit).
An example: the Sackur-Tetrode formula
Say we're trying to obtain the Sackur-Tetrode formula for entropy. Starting from the parition function for a $3D$ ideal gas,
$$Z=\frac{V^N}{\lambda^{3N}N!}$$
and using the formula for entropy $S = \frac{\partial}{\partial T}(kT \log{Z})$ one gets exactly
$$S =N k\log{\frac{V}{\lambda^{3}}}+\frac{3}{2}Nk-k\log{N!}$$
Now we can make the approximation
$$\log{N!} = N\log{N}+\mathcal{O}(N)$$
and we get:
$$S =N k\log{\frac{V}{\lambda^{3}N}}+\mathcal{O}(N)\tag{1}$$
which is the correct formula to $\mathcal{O}(N)$ because $V$ scales like $N$. If instead we took one further term:
$$\log{N!} = N\log{N}-N+\mathcal{O}(1)$$
The approximation would be more precise:
$$S =N k\left[\log{\frac{V}{\lambda^{3}N}}+\frac{5}{2}\right]+\mathcal{O}(1)\tag{2}$$
Equation $(2)$ is correct to one further order, but depending on what you're looking at $(1)$ may be just as good. Note that it's pointless to take further orders (say $\mathcal{O}(1/N)$) because we're in the large $N$ regime.
Another example: deriving the quantum statistics
I assume the derivation goes like in this question. You're trying to maximise an entropy of the form:
$$S = \log{\frac{N!}{n_i ! (N-n_i)!}}$$
subject to certain constraints. In order to take derivatives of $S$ we must apply Stirling's approximation:
$$S = \log{N!} - n_i \log{n_i}-(N-n_i) \log{(N-n_i)}+\mathcal{O}(n_i)$$
or to one further order after cancellation:
$$S = \log{N!} - n_i \log{n_i}-(N-n_i) \log{(N-n_i)}+N+\mathcal{O}(1)$$
In this case the correction is actually independent of $n_i$, so when taking derivatives both approximations yield the same result. However you do not know beforehand whether this is case, so you should use the $\mathcal{O}(1)$ approximation.
EDIT:
When doing statistical mechanics we almost always work in the so-called thermodynamic limit: $N,V \to \infty$ while keeping $N/V$ constant. This means that extensive quantities (for instance the entropy), which scale like $N$, will go to infinity. In order to extract useful results, it is useful to define stuff like "entropy per particle" $s = S/N$. This remains finite in the thermodynamic limit. It is clear at this point that the only terms which contribute in the thermodynamic limit are at least $\mathcal{O}(N)$, and subleading terms are irrelevant.
A: Let's start by writing two more terms in Stirling's approximation: $$\log N! = N \log N - N +\frac{1}{2}\log N + \frac{1}{2}\log 2\pi + \mathcal{O}(\frac{1}{N}).$$ 
On an intuitive level, the main point is that $\log$ increases very slowly, and thus for the typical numbers of statistical mechanics, i.e., Avogadro's number, $\log 10^{24} = 24 \log 10$ which is about 24, i.e., much much less than $10^{24}$ and not much more than 1. Thus, $\log N$ can be safely neglected compared to $N$, but if $\mathcal{O}(\log N)$ terms are kept it is often prudent to also keep $\mathcal{O}(1)$ terms. 
In fact, statistical mechanics is almost always accurate only "to leading order in the exponent" (before taking logs, i.e., when calculating partition functions or the number of states). So in the thermodynamic variables, terms of order $\log N$ are always neglected (think, e.g., of the difference between the Boltzmann and Gibbs entropies, or the difference between the log of the partition function and the Legendre transform of the entropy). The more mathematically sophisticated (or at least jargon-heavy) way of saying this uses the language of the theory of large deviations: thermodynamic potentials such as entropies and free energies are rate functions, defined (at least implicitly) by limits of the form $S(E) = \lim_{N\to\infty} \frac{1}{N} \log \Omega(E)$ (where $\Omega(E)$ is the number of states with energy $E$). In these limits, the higher order terms in Sterling's approximation do not contribute.
A further comment on why both $N \log N$ and $N$ are kept and not only the first of these: usually, thermodynamic quantities are either extensive ($\mathcal{O}(N)$) or intensive ($\mathcal{O}(1)$). So the $N \log N$ terms usually cancels out in the calculation and the leading term is $\mathcal{O}(N)$. This is what happens, e.g., in the Sackur-Tetrode example given by John Donne, or when approximating bynomial coefficients: $\log {{N}\choose{xN}}\simeq N[x \log x + (1-x)\log(1-x)]$ (where $0<x<1$). (In fact, when the $N\log N$ term does not cancel, people often subtract it by hand, e.g. by choosing to work with relative entropies rather than entropies). 
