In An Introduction to Quantum Field Theory by Peskin and Schroeder, section 9.2, they calculate the four-point correlation function for a free real scalar field $\phi(x)$ using the path integral formulation: $$ \langle \Omega |T\{\phi(x_1) \phi(x_2) \phi(x_3) \phi(x_4)\}|\Omega\rangle = \lim_{T\to\infty(1-i\epsilon)} \frac{\int \mathcal{D}\phi\ \phi(x_1) \phi(x_2) \phi(x_3) \phi(x_4) \exp \left[i\int d^4x \mathcal{L}\right]}{\int \mathcal{D}\phi\ \exp\left[i\int d^4x \mathcal{L}\right]} $$
Setting the fields to be defined on the points $x_i$ in a discrete lattice of volume V and performing a Fourier series $\phi(x_i)\to \phi_n(k_n)$, he arrives at the following:
$$ \text{Numerator} = \frac{1}{V^4}\sum_{m,l,p,q}e^{-i(k_m\cdot x_1 + k_l\cdot x_2+ k_p\cdot x_3+ k_q\cdot x_4)}\left( \prod_{k_n^0 >0} \int d\text{Re} \phi_n \, d\text{Im} \phi_n \right)\times \\ \times (\text{Re}\phi_m + i\text{Im}\phi_m) (\text{Re}\phi_l + i\text{Im}\phi_l) (\text{Re}\phi_p + i\text{Im}\phi_p) (\text{Re}\phi_q + i\text{Im}\phi_q)\times \\ \times \exp\left[ -\frac{i}{V}\sum_{k_n^0>0} (m^2-k_n^2)\left[(\text{Re}\phi_n)^2 + (\text{Im}\phi_n)^2\right] \right] $$
Now, most of the terms vanish, since the integrand would be odd, but there are some values of $m,l,p,q$ for which the integral is non-zero. For example, if $k_l = -k_m$ and $k_q = -k_p$, then the fact that $φ (x_i)$ is real, i.e. $\phi(k_n) = (\phi(-k_n))^*$, means that we get:
$$ \frac{1}{V^4}\sum_{m,p}e^{-ik_m\cdot (x_1 - x_2)}e^{i k_p\cdot (x_3+ x_4)}\left( \prod_{k_n^0 >0} \int d\text{Re} \phi_n \, d\text{Im} \phi_n \right)\times \\ \times \left[(\text{Re}\phi_m)^2 + (\text{Im}\phi_m)^2\right] \left[(\text{Re}\phi_p)^2 + (\text{Im}\phi_p)^2\right]\times \\ \times \exp\left[ -\frac{i}{V}\sum_{k_n^0>0} (m^2-k_n^2)\left[(\text{Re}\phi_n)^2 + (\text{Im}\phi_n)^2\right] \right] = \\ = \frac{1}{V^4}\sum_{m,p}e^{-ik_m\cdot (x_1 - x_2)}e^{i k_p\cdot (x_3+ x_4)}\left( \prod_{k_n^0>0}\frac{-i\pi V}{m^2-k^2_n} \right)\frac{-iV}{m^2-k_m^2 - i\epsilon} \frac{-iV}{m^2-k_p^2 - i\epsilon}, $$
which, when returning to the continuum, becomes:
$$ \left( \prod_{k_n^0>0}\frac{-i\pi V}{m^2-k^2_n} \right)D_F(x_1-x_2)D_F(x_3-x_4), $$
where $D_F(x-y)$ is the Feynman propagator. However, I have a problem with this argument. For the values corresponding to $k_m = k_p$, we would have a quartic term $\left[(\text{Re}\phi_m)^2 + (\text{Im}\phi_m)^2\right]^2 = (\text{Re}\phi_m)^4 + (\text{Im}\phi_m)^4 + 2 (\text{Re}\phi_m)^2 (\text{Im}\phi_m)^2 $ in the integral which, using the results from Gaussian integration:
$$ \int dx\ x^2 e^{-ax^2} = \frac{1}{2a}\sqrt{\frac{\pi}{a}}, \quad \int dx\ x^4 e^{-ax^2} = \frac{3}{4a^2}\sqrt{\frac{\pi}{a}}, $$
would mean that we get an additional factor of 2 only for the case when $k_m = k_p$. This would mess up the result with the Feynman propagators, right? Can anyone see where my mistake lies?