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=\sum_E g(E)e^{-\beta E}=\sum_{p,q} g(p,q)e^{-\beta H(p,q)}$$
Now the interesting thing and also a crucial step is to identify the the degeneracy $g(p,q)$ with 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 $$g(p,q)=\frac{1}{N!}\frac{dpdq}{w_0}=\frac{1}{N!}\frac{dpdq}{h^{3N}}$$ putting this back to the summation and converting into integration, we readily obtain the desired integral form for the partition function: $$Z=\frac{1}{N!h^{3N}}\int e^{-\beta H(p,q)}dqdp$$