The factorial factor $1/N!$ is exact but does not apply to all statistics.
Consider a two level system and let us call $\xi_i$ the grand-canonical partition function for the energy level $i$ ($i=0$ or $1$) and $z=\mathrm e^{\beta\mu}$ the fugacity. Keep in mind that $\xi_i$ is the partition function for a given energy level.
Classical particles are independent, indistinguishable, $\xi_i$ is computed as a usual partition function for undistinguishable particles. Here we divide by the symmetry $p!$ because the particles are undistinguishable and therefore the order in which we could label them should not matter.
$$\xi_i^{\text{classical}}=\sum_{p=0}^\infty \frac1{p!}z^p \mathrm e^{-p\beta E_i}=1+z\mathrm e^{-\beta E_i}+\frac{z^2}{2}\mathrm e^{-2\beta E_i}+\cdots=\exp\left(z\mathrm e^{-\beta E_i}\right).$$
For quantum particles, the sum runs over the number of particles which can occupy the energy level. There is no division by $p!$ because the quantum states already are (anti)-symmetric; for instance the quantum state where there is one fermion in each level is $\frac1{\sqrt2}\left(\left|01\right\rangle-\left|10\right\rangle\right)$.
$$ \xi_i^{\text{fermions}}=\sum_{p=0}^1z^p\mathrm e^{-p\beta E_i}=1+z\mathrm e^{-\beta E_i}$$
$$ \xi_i^{\text{bosons}}=\sum_{p=0}^\infty z^p\mathrm e^{-p\beta E_i}=\frac1{1-z\mathrm e^{-\beta E_i}}$$
The grand-canonical partition functions are
$$ \mathcal Z^{\text{classical}}=\xi_0^{\text{classical}}\xi_1^{\text{classical}}=\mathrm e^{z\left(1+\mathrm e^{-\beta E}\right)}
=\mathrm e^{zZ_1}=1+zZ_1+\frac{z^2}2Z_1^2+\cdots$$
$$\mathcal Z^{\text{fermions}}=\xi_0^{\text{fermions}}\xi_1^{\text{fermions}}=(1+z)\left(1+z\mathrm e^{-\beta E}\right)=1+zZ_1+z^2\mathrm e^{-\beta E}$$
$$\mathcal Z^{\text{bosons}}=\xi_0^{\text{bosons}}\xi_1^{\text{bosons}}=
\frac1{1-z}\frac1{1-z\mathrm e^{-\beta E}}=1+zZ_1-z^2\mathrm e^{-\beta E}+z^2Z_1^2+\cdots$$
with $Z_1=1+\mathrm e^{-\beta E}$ is the one-particle partition function.
Now looking at the coefficient $z^2$ in $\mathcal Z$ gives the two-particle canonical partition function $Z_2$. We have
$$ Z_2^{\text{classical}}=\frac{1}{2}Z_1^2,\quad Z_2^{\text{fermions}}=\mathrm e^{-\beta E},\quad Z_2^{\text{bosons}}=1+\mathrm e^{-\beta E}+\mathrm e^{-2\beta E}.$$
In none of these case was it necessary to make an approximation. Undistinguishability in the quantum cases is for quantum states, not for particles.
Remark
The difference between $Z_2^{\text{classical}}$ and $Z_2^{\text{bosons}}$ is evidence for the fact that the bosons and the classical particles have a difference. Indeed, classical particles have no correlations, which is expressed by the fact $Z_2\propto Z_1^2$, whereas bosons have correlations: compared to the uncorrelated classical particles, the relative weight of the states where the particles are both at the same energy level is larger: bosons prefer to be in the same state. (Fermions avoid being at the same energ level and classical particles don't care.)
You ask
how can you formally go from the summation over $s$ to the double sum over $s_1$ and $s_2$?
As we'll see in a moment, passing to the double sum relies on the mathematical fact (about tensor products of Hilbert spaces) that if $s_1$ labels a basis of states for system $1$, and if $s_2$ labels a basis of states for system $2$, then the set of all pairs $(s_1, s_2)$ labels a basis of states for the composite system.
Let's assume that the system at hand is a quantum system described by a state space (Hilbert space) $\mathcal H$. Let the (state) labels $s$ index an orthonormal basis of $\mathcal H$ consisting of eigenvectors $|s\rangle$ of the total system's Hamiltonian $H$;
\begin{align}
H|s\rangle = E(s)|s\rangle,
\end{align}
then the partition function is obtained by summing over these states;
\begin{align}
Z = \sum_s e^{-\beta E(s)}.
\end{align}
Notice that the state labels $s$ don't have to be numbers. They could, for example, be pairs of numbers or triples of numbers, whatever is most convenient to label the possible states for the system at hand.
Now, suppose that the system consists of a pair of subsystems $1$ and $2$, then the Hilbert space of the combined system can be written as a tensor product of the Hilbert spaces for the individual subsystems;
\begin{align}
\mathcal H = \mathcal H_1\otimes \mathcal H_2.
\end{align}
Let $H_1$ denote the hamiltonian for subsystem $1$, and let $H_2$ denote the hamiltonian for subsystem $2$. Let $\{|1, s_1\rangle\}$ be an orthonormal basis for $\mathcal H_1$ consisting of eigenvectors of $H_1$, and let $\{|2, s_2\rangle\}$ be an orthonormal basis for $\mathcal H_2$ consisting of eigenvectors of $H_2$, then we have the following basic fact;
The tensor product states
\begin{align}
|s_1, s_2\rangle := |1,s_1\rangle\otimes|2,s_2\rangle
\end{align}
yield an orthonormal basis for the full state space $\mathcal H = \mathcal H_1\otimes \mathcal H_2$.
In particular, the set of states one needs to sum over in the partition function can be enumerated by pairs $s=(s_1, s_2)$. Moreover, if the systems are non-interacting, then the Hamiltonian of the full system is essentially the sum of the Hamiltonians of the individual subsystems (with appropriate identity operator factors);
\begin{align}
H = H_1\otimes I_2 + I_1\otimes H_2,
\end{align}
so that if $E_1(s_1)$ and $E_2(s_2)$ denote the energy eigenvalues of $H_1$ and $H_2$ corresponding to the states $|1,s_1\rangle$ and $|2, s_2\rangle$, then the energy $E(s_1, s_2)$ of a tensor product basis state is their sum;
\begin{align}
H|s_1,s_2\rangle = E(s_1, s_2)|s_1, s_2\rangle = \big(E_1(s_1) + E_2(s_2)\big)|s_1,s_2\rangle,
\end{align}
and the partition function can be written as a sum over the tensor product basis states;
\begin{align}
Z = \sum_{(s_1,s_2)} e^{-\beta E(s_1, s_2)} = \sum_{s_1}\sum_{s_2} e^{-\beta\big(E_1(s_2) + E_2(s_2)\big)} = \sum_{s_1} e^{-\beta E_1(s_1)} \sum_{s_2} e^{-\beta E(s_2)} = Z_1Z_2
\end{align}
where in the second equality we have used the fact that a single sum over all possible pairs $(s_1, s_2)$ is equivalent to iterated sums over the possible values of $s_1$ and $s_2$.
Best Answer
In the example you gave above and in calculating certain thermodynamic variables such as pressure: $$P=-\frac{dF}{dV}=\frac{d}{dV}k_bT\ln(Z)$$ the $N!$ cancels: $$\frac{d}{dx}\ln(Z'/N!)=\frac{d}{dx}\left(\ln(Z')-\ln(N!)\right)=\frac{d}{dx}\ln(Z')$$
you still need to included the $N!$ when you do calculations for chemical potential or entropy because the $N!$ doesn't cancel. What you read above is just short hand for the calculation they were doing.