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.)
Equilibrium statistical mechanics is not really about the partition function per se, it just aims at finding (according to one interpretation at least) what is the least biased probability distribution (for the microstates) that satisfies the known constraints on the system.
In your case, you assumed that the probability to get either head or tail for each coin is 1/2 which implicitely uses a statistical inference reasoning to get the probability a priori (for instance Pascal's principle of indifference).
In my view, stat. mech. is concerned with finding those probabilities a priori for each statistical ensemble (microcanonical, canonical, grand canonical and so on). Once this is done, you just apply the usual tools of probability theory to compute averages, variances etc..
In particular, the partition function is closely related to the probability generating function in mathematics.
Let us look at the canonical ensemble describing a statistical ensemble where a system has a conserved volume, number of particles and is in contact with a thermostat at inverse temperature $\beta$.
Equilibrium statistical mechanics tells us that the probability for the system to be in a given microstate $m$ is:
$p(m|\beta) = \frac{e^{-\beta E_m}}{Z(\beta)}$
where normalization of this probability yields:
$\sum_m p(m|\beta) = 1 \: \Rightarrow \: Z(\beta) = \sum_m e^{-\beta E_m}$
Let us now consider a macrovariable $X(m)$ which is a function of the microstate $m$. Its probability distribution can be seen as the marginal distribution of the microstate probability distribution. If the random variables are discrete for simplicity:
$p(X=x|\beta) = \sum_m p(m|\beta) \delta_{X(m),x}$ where $\delta_{x,y}$ is a Kronecker delta function which is one only when $x=y$.
From that point on we can compute the moments for $X$ by using either the probability generating function $G(z) \equiv \mathbb{E}_X[z^X]$ or the moment generating function $M(t) \equiv \mathbb{E}_X[e^{tX}]$ from which we can the moments bt taking derivatives with respect to $z$ or $t$ when the latter are evaluated at $1$ and $0$ respectively.
When $X$ is the energy itself, then the sum over microstates having an energy $E$ can be factorized by $e^{-\beta E}$ and we get:
$p(E=e|\beta) = \frac{g(e)e^{-\beta e}}{Z(\beta)}$
Because both $g(e)$ and $e$ are assumed independent of $\beta$, we can actually notice that all the moments:
$\mathbb{E}_E[E^n] \equiv \sum_e \:e^n \frac{g(e)e^{-\beta e}}{Z(\beta)}$ can be gotten by looking at the quantity $\frac{(-1)^n}{Z(\beta)}\frac{\partial^n Z(\beta)}{\partial \beta^n}$ which is very similar to what we would have done with either the PGF or the MGF or even the CGF.
The reason why, when it comes to moments of the energy variable, we prefer doing this way instead of using the other generating function tools is first because we can (thanks to the fact that both $g(e)$ and $e$ do not depend on $\beta$ or at least we assume so) and second because there is a link with derivatives of the free energy in thermodynamics and energy fluctuations.
Otherwise, in general, it is quite common to introduce an "external field" coupled linearly with the variable $X$ in the energy of the system (which is almost always possible if we imagine the field strength small enough) and that plays the same role as the variable $t$ in the MGF or the CGF approaches. In magnetic systems for instance, it is very common to introduce an external magnetic field to get the moments of the magnetization.
Best Answer
That definition of $\Omega$ only works if you have distinguishable particles and you can put as many as you want in any state. That gives you the Maxwell-Boltzmann distribution.
In quantum mechanics, particles of the same type are not distinguishable; e.g. there's no way to tell one electron from another. Moreover, for Fermions (like electrons), you can have at most one particle in a given state. Bosons (like photons) can have as many particles in one state as you please.
These give rise to two '$\Omega$'s, both different than the one you gave.
As for the partition function, it takes the same form, but you sum over different things. Define
$$Z=\sum_se^{-E_s\beta}$$
where the sum is over the possible microstates of the system. (A state is an energy level that a particle can occupy; e.g. the ground state of a hydrogen atome. A microstate is a configuration of all the particles in the system; some examples are given below.)
Say you have three particles and three states available to you. For Bosons, you have many different microstates available to you: all three particles in state one, one particle in each state, all three particles in state three, etc. So, your sum has many terms. For Fermions, you have only one option: one particle in each state, since there are three particles, three states, and only one particle can be in each state. So there is only one microstate available for Fermions in this system, and thus the sum only has one term.
Short version: the equation for the partition function stays the same, but the microstates you can sum over are very different.