= Solution
The Bose--Einstein and Fermi--Dirac occupation factors are
$$
\frac1{\exp((E-\mu)/(k_BT))\mp1}.
$$
In the stated nonrelativistic dilute regime, $E-\mu$ is dominated by $mc^2$, so the exponential is large and either denominator is asymptotic to $\exp((E-\mu)/(k_BT))$. Integrating the resulting Maxwell--Boltzmann occupation over momentum states gives
$$
n=\frac{4\pi g_s}{h^3}\int_0^\infty
p^2e^{-[E(p)-\mu]/(k_BT)}\,dp.
$$
For $E=mc^2+p^2/(2m)$, integration by parts or differentiation of the supplied Gaussian integral gives
$$
\int_0^\infty p^2e^{-p^2/(2mk_BT)}\,dp
=\frac{\sqrt\pi}{4}(2mk_BT)^{3/2}.
$$
Consequently
$$
\boxed{
n=g_s\left(\frac{2\pi mk_BT}{h^2}\right)^{3/2}
e^{(\mu-mc^2)/(k_BT)}}.
$$
Back to article page