When computing partition functions for classical systems with $N$ particles with a given Hamiltonian $H$ I've seen some places writing it as
$$Z = \dfrac{1}{h^{3N} N!}\int e^{-\beta H(p,q)}dpdq$$
where the integral is over the available phase space. Now, in my book there is not this factor $1/h^{3N} N!$. Where this factor comes from? Why do we need to include it there?
The factor of $1/h^{3N}$ is a total hack. The integral over phase space has dimensions, whereas $Z$ only makes sense if it's dimensionless. The $h$ factors are there to make $Z$ dimensionless. Suppose you have a system with only one particle in one dimension. Then the integral in phase space goes over one position variable, $dq$ and one momentum variable $dp$. The dimensions of position times momentum, i.e. $dpdq$, is called action.$^{[a]}$ So, to make $Z$ dimensionless in the one-particle one-dimensional case, you have to divide by something with dimensions of action. Planck's constant is a perfectly acceptable choice because it happens to have dimensions of action.$^{[b]}$ However, as much as it's acceptable it's also totally arbitrary and should make you highly suspicious that someone is trying to fool you or not explain what's really going on.
Now, if you go to three dimensions the integral over phase space has three factors of length and three factors of momentum. Therefore, the dimensions of the integral are $\text{action}^3$. To get rid of those dimensions you divide by $h^3$.
When you have $N$ particles you have $N$ integrals and so you have $3N$ powers of $h$.
The thing is,when you really come down to it, $Z$ doesn't actually have to be dimensionless. Whenever you compute anything physical you wind up dividing $Z$ by itself. Suppose we define $Z$ without the $h$ factors (and without the $N$ factors for now): $$Z \equiv \int e^{-\beta H(q,p)}\,dqdp \, .$$ Consider the average energy of the system \begin{align} \langle H \rangle &= \frac{1}{Z} \int H(p,q) e^{-\beta H(q,p)} \, dqdp \\ &= \frac{1}{Z} \int - \frac{d}{d\beta} e^{-\beta H(q,p)} \, dqdp \\ &= - \frac{1}{Z} \frac{d}{d\beta} \int e^{-\beta H(q,p)} \, dqdp \\ &= - \frac{1}{Z} \frac{dZ}{d\beta} \, . \end{align} The dimensions of $Z$ cancel out! This always happens. So, in classical physics the $h$ factors used to remove the dimensions for $Z$ don't really need to be there.
The same thing happens when you deal with entropy. At first, it seems like it has the wrong dimensions because it's the logarithm of the available phase space volume. If phase space volume has dimensions you can't take the log! But when you realize that only ratios of available phase space volumes, and therefore differences in entropy that actually mean anything in a classical system, it's no longer an issue.$^{[c]}$
The $N$ factor is more interesting. First, you need to understand the Gibb's paradox. The basic idea is that if I put a divider in the middle of a box of gas molecules, I cut the available positions of each one in half, so the entropy would seem to be reduced. Think of it like this: there's less possible ways for the system to be arranged if I add the divider, because molecule $A$ is restricted to one half of the box. The problem is that I now violated the 2nd law of thermo by reducing the entropy.
Dividing $Z$ by $N!$ fixes this. The number $N!$ is the number of ways a set of $N$ things can be permuted among themselves. So, dividing $Z$ by $N!$ is basically saying to treat the states related by swapping molecules with one another as the same state. You are un-counting states which are related by particle interchange. Another way to say this is to say that you're treating the particles as "indistinguishable" because a situation with particle $A$ in state 1 and particle $B$ in state 2 is the same as particle $A$ in state 2 and particle $B$ in state 1. A better way to say this is that the notion of $A$ and $B$ doesn't have meaning, and instead you should just say "the number of particles in state 1 is 1 and the number of particles in state 2 is 1".
This indistinguishibility is really weird in the context of classical mechanics. Where did it come from? Why can't we talk about two molecules as having their own identity? The answers come from quantum mechanics where you find that "particles" are really best thought of as excitations of modes.
I highly recommend reading this other Physics.SE post.
$[a]$: The SI unit of action is the Joule-second. You can think of that as a momentum times a length, as we saw with $dqdp$, or an energy times a time.
$[b]$: Of course, this is not an accident at all, but since we're not really talking about quantum mechanics here I don't get into it.
$[c]$: If someone wants to point out that this is wrong, please do!
For the $h$ factor: (the $N!$ is the Gibbs correction as already explained by others )
Step 1: Semi-classical setting
The $h$ factor can actually be understood in a semiclassical way. We are talking about classical ideal gas using classical phase space formulation of classical statistical mechanics, but still from a semi-classical point of view we are not able to pinpoint exactly each point in the phase space simply due to the uncertainty relation: $\Delta p \Delta q \sim h$. (In fact other than this, we have also used quantized energy levels for classical ideal gas in a box as another part of the semi-classical approach.) This basically means each phase point effectively occupies some small volume $w_0$ in the phase space, which is called the fundamental volume by some authors. If we can some how figure out $w_0$, then the total number of relevant phase points would just be $\Omega =\Gamma /w_0$, where $\Gamma$ denotes the total phase space volume for some very small energy range $(E,E+\Delta)$.
Step 2: Fundamental volume $w_0$ defined in microcanonical ensemble
Note that the reason an energy range is used is that if the energy is at a particular fixed value, then the relevant phase space would just be a hypersurface with volume 0, and the reason that the range is chosen to be small is that we will still be able to work in the microcanonical ensemble, where each microstate is equally likely, which validates the introduction of the fundamental volume. (meaning we need a uniform phase space density to have the same $w_0$ for all the phase points)
Now for your system with $N$ particles, the phase space volume can be calculated as $$\Gamma = \frac{1}{N!}\int_{E\leq \sum_i {p_i}^2/2m \leq E+\Delta} dpdq=\frac{\Delta}{E}\frac{(2\pi mE)^{3N/2}V^N}{(3N/2-1)!N!}$$ On the other hand, using the conventional $N$-particle in a box with total energy varying from $E$ to $E+\Delta$, we can calculate the total number of microstates to be $$\Omega=\frac{\Delta}{E}\frac{(2\pi mE)^{3N/2}V^N}{(3N/2-1)!N!h^{3N}}$$ Consequently, the fundamental volume will be given by $$w_0=\frac{\Gamma}{\Omega}=h^{3N}$$ i.e. the fundamental volume will be different for phase space of different dimensions.
Step 3: From microcanonical multiplicity $\Omega$ to canonical partition function $Z$
Partition function is just summation over states. In the discrete case (including degeneracy): $$Z=\frac{1}{N!}\sum_E g(E)e^{-\beta E}=\frac{1}{N!}\sum_{p,q}e^{-\beta H(p,q)}$$
When going from the discrete sum to the continuous integration, we have to consider the number of microstates (or equivalently phase points) contained in a small phase space volume $dqdp$ around the phase point $(p,q)$. Then we immediately have $$\sum_{p,q} = \int \frac{dpdq}{w_0} = \int \frac{dpdq}{h^{3N}}$$ after which we readily obtain the desired integral form for the partition function: $$Z=\frac{1}{N!h^{3N}}\int e^{-\beta H(p,q)}dqdp$$
Now, in my book there is not this factor $1/h^{3N} N!$. Where this factor comes from? Why do we need to include it there?
$\mathbf{1/h^{3N}}~$:
Some people like the value of the partition integral to be independent of units, so they add the denominator $h^{3N}$ to the formula where $h$ has the same units as $pq$. This makes the expression dimensionless.
However, this has no consequence for the significant results of calculations based on this formula. For example, the probability distribution of velocities of ideal gas particles (the Maxwell-Boltzmann distribution) can be derived with help of the above expression for $Z$, whether the denominator $h^{3N}$ is used or not; the resulting distribution does not depend on $h$ (see also Daniel Sank's answer). Consequently, $h$ may be given any value.
Since Planck introduced quantity of the same dimensions, (Planck constant $h_{Planck}=$6.6E-34 J.s), people naturally tend to put Planck's $h_{Planck}$ for $h$.
Another reason for the use of the $h$-dependent denominator is the idea that the phase space is a wrong concept and should be used only as an approximate tool to evaluate the number of discrete states. This number is often estimated as $(phase~space~volume)/h^{3N}$ (the last expression has generally real-number value, so it is only an estimation of the count, which should be natural number). Value of $h$ has some significant consequences in quantum statistical physics, I think - Jaynes mentions chemical constants.
I never found this procedure of imaginary phase cells to be convincing and reliable - it may work in some cases, but in its core it seems to mix two things that are just too different to play well together (distributions on phase space and the idea of discrete states).
Good read with some related information (section 9.):
Now, in my book there is not this factor $1/h^{3N} N!$. Where this factor comes from? Why do we need to include it there?
$\mathbf{1/N!}~$:
There are two basic but different concepts of $Z$. The first is the partition function of classical statistical physics:
$$ Z_{int} = \int dq_1...dq_N \int dp_1 ... dp_N~e^{-\beta H(\mathbf q,\mathbf p)}. $$
There is no $h$ and no $N!$. The meaning of $Z_{int}$ is that it is the normalization factor in the formula for probability density in the phase space:
$$ \rho(\mathbf q,\mathbf p) = \frac{e^{-\beta H(\mathbf q,\mathbf p)}}{Z_{int}}. $$
The other meaning is that $Z$ is a sum over all distinct discrete states $k$ with energy $\epsilon_k$:
$$ Z_{sum} = \sum_k e^{-\beta \epsilon_k}. $$ There is no $h$ and no $N!$ here either. This second formula is often easier to work with.
There is a view that the first formula $Z_{int}$ is deficient in that when it is discretized so that one calculates a sum instead of an integral, it counts the states of the system related by permutation of the particle states in the sequence 1,...,N $N!$ times. Since the view regards these as the same state of the system, it should, according to this view, be counted only once in the formula $Z_{sum}$. Hence the modifying factor $1/N!$ is added to the integral, to suppress the contribution of the superfluous permutations.
This is connected to change of statistical formula for entropy too. In order to preserve the relation
$$ Z= e^{-\beta (U - TS)} $$
if $Z$ changes to $Z'=Z/N!$, also some quantity on the right needs to change. This quantity is entropy, which changes into $S'=S-k_B\ln N!$ so that the relation
$$ Z' = e^{-\beta (U - TS')} $$ holds. Whether one uses $Z,S$ or $Z',S'$, the results of the calculations are the same. However, if connection to conventions of thermodynamics is to be made, $S'$ is closer to thermodynamic entropy $S_{th}$, because they are both additive. The entropy based on the formula
$$ Z= e^{-\beta (U - TS)} $$ (with $Z$ summation including all permutations without reduction) is not additive.
