Is the $N$ factorial in the Partition function for $N$ indistinguishable particle an approximation? I suspect that the $N$ factorial in the partition function for N indistinguishable particles 
$$ Z = \frac{ Z_0^N } {N!} $$
is an approximation. Please someone correct me if I am wrong and why or why not. Thanks.

A simple case: 
each particle has two states with energy $0$ and $E$. The partition function for a single particle is 
$$ Z_0 = 1 + e^{- \beta E} . $$
If there are only two particles, there is the total partition function
$$
Z = \frac{ Z_0^2 } {2}.
$$
But regarding the whole system consisting of these two particles, we can also write
$$
Z = 1 + e^{- \beta E}  + e^{-2 \beta E} .
$$
And it is certain that 
$$
  \frac{ Z_0^2 } {2} \neq  1 + e^{- \beta E}  + e^{-2 \beta E}
$$
 A: 
In the figure above, consider the different configurations that are possible with 3 particles and 5 energy levels. Dividing by 3! gets the symmetry factor correct only for configurations of type 1 but is wrong for configurations of type 2 and 3. You can see this by explicitly writing out $Z$ and comparing with $z^3/3!$. That is why the OP's statement that $z^3/3!$ is an approximation is correct. I can add details to this, if necessary. (Notation: The lower-case $z$ is the single particle partition function)
To add to Josh's remark above, even for fermions where terms of type I are only allowed, the expansion of $z^N/N!$ contains terms of type 2 and 3 which are not present in $Z$. Nevertheless, the dominant contribution at large $N$ (as well as number of energy levels) is from terms of type 1. Hence my statement that it is a fairly good approximation holds.
A: 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.)
A: It is indeed an approximation, in particular an approximation that is usually good at high enough temperatures.
For $N$ distinguishable, non-interacting particles the partition function is $Z(\mathrm{dist.}) = {Z_0}^N$, where $Z_0$ is the single-particle partition function. If the $N$ particles are all in different quantum states then there are $N!$ microstates of distinguishable particles which all correspond to the same microstate of indistinguishable particles. However, if the particles are not all in different states then the number of microstates of distinguishable particles corresponding to the same microstate of indistinguishable particles is less than $N!$.
So in general the relationship between the sum over microstates for distinguishable particles and for indistinguishable particles is quite complex. However, if the temperature is high enough then (for most systems) the number of accessible quantum states will be much larger than the number of particles, so it is unlikely more than one particle will occupy the same quantum state. In this case we can make the approximation that the partition function for distinguishable particles is just $N!$ times that for indistinguishable particles. I.e. $Z(\mathrm{indist.}) = {Z_0}^N/N!$ as you wrote.
The claim that at high enough temperatures the number of accessible quantum states will be much larger than the number of particles is not true for all systems, for example it is false for those with bounded energy spectrum (such as the ideal paramagnet or the example the OP gives). Also, how high the temperature needs to be depends upon the system. For an ideal gas at atmospheric densities it's a good approximation for temperatures much above 1K.
The above is true for quantum particles (Fermions or Bosons). In the classical limit the quantum energy-level separation is taken to zero and the spectrum becomes continuous. For a continuous spectrum the probability that two particles are in exactly the same state is zero, and so the approximation becomes exact in this limit. Note however, that there is no such thing as a classical two-level system, such as the example given by the OP. A classical system with two local minima can in some cases be treated as a two level system, but this is in itself an approximation.
A: It is indeed an approximation. As suresh says, the approximation holds only when each particle is in a different single-particle state, or more precisely, when the number of states in which each particle is in a different single-particle state vastly outweighs the number of states in which more than one particle is in the same SP state.
The single particle partition function is
$$ Z_0 = 1 + e^{-\beta E}$$
The two distinguishable particle partition function is
$$ Z_\mathrm{dist} = 1 + 2e^{-\beta E} + e^{-2\beta E} = Z_0^2$$
The two indistinguishable particle partition function is
$$ Z_\mathrm{in} = 1 + e^{-\beta E} + e^{-2\beta E}$$
You can hopefully see that $Z_\mathrm{in}$ is not equal to $Z_0^2/2!$ from the above equations. The logic for the $N!$ division is as follows: for each overall state of the two distinguishable particle system there are $N!$ ways of permuting the particles, but for identical particles these are all indistinguishable. However, for any overall state of the two distinguishable particle system with two particles in the same single-particle state, we are already considering the possible permutations as indistinct. Even for distinguishable particles, "particle A in state X and particle B in state X" is considered to be the same state as "particle B in state X and particle A in state X", as should be fairly obvious if you think about it! So only for those overall states for which all particles are in different single-particle states is the full division by $N!$ appropriate. 
For systems with enormous numbers of states --- such as systems with translation degrees of freedom --- this is a very good approximation, since the typical number of available single-particle states vastly outnumbers the typical number of particles. On the other hand, high temperature is not generally a limit in which this approximation holds. You should be able to see this from the above expressions for partition functions. Only for systems with an infinite number of single-particle states does the approximation hold in this limit.
A: Actually, it's exact. The flaw is "regarding the whole system consisting of these two particles, we can also write" $Z = 1 + e^{- \beta E}  + e^{-2 \beta E}.$ Assuming the two particles are distinguishable, we have 
$$Z=\sum_ig_ie^{-\beta E_i}=1+2e^{-\beta E}+e^{-2\beta E}=Z_0^2,$$
with the $2e^{-\beta E}$ since the state of energy $E$ is doubly-degenerate.
The additional factor of 1/2 in $Z=Z_0/2!$ accounts for indistinguishability.
EDIT: This seems to be wrong. See other answers.
