**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$$