Distribution of the output of AWGN channel
I see this question is open for a long time, I am not sure if I understand it correctly.
PDF of the symbols
The first probability is a $P_{\mathbf{Y}}: \mathbb{C}^n \to \mathbb{R}$, you are not considering phase noise, the probability  of the $n$ symbols are identically distributed, let's say with $P_{\mathbf{Y}_i}: \mathbb{C}^n \to \mathbb{R}$, since the joint probability of independent variables are the product of the individual variables this distribution can be written as
$$P_{\mathbf{Y}}(y) = \prod_{i=1}^{N}P_{\mathbf{Y}_i}(y_i)$$
Where $P_{\mathbf{Y}_i}(y_i) = \sum_{x\in \mathbf{X}_i}p_{\mathbf{X}}(x_i) e^{-||y_i - x_i||^2}/\pi$ which leads to
$$P_{\mathbf{Y}}(y) = \frac{1}{\pi^N}\prod_{i=1}^{N} \sum_{x\in \mathbf{X}_i} p_{\mathbf{X}}(x_i) e^{-||y_i - x_i||^2}$$
The term $p_{\mathbb{X}}(x_i)$ can be used to model probabilistic shaping, but for the original quesiton it will probably be a the inverse of the number of symbols in the constellation $1/c = 1/|\mathbf{X}_i|$, that by it's turn is the same for all symbols, thus
$$P_{\mathbf{Y}}(y) = \frac{1}{c^N\pi^N}\prod_{i=1}^{N} \sum_{x\in \mathbf{X}_i} p_{\mathbf{X}}(x_i) e^{-||y_i - x_i||^2}$$
The complex of evaluating this will be $c N$ exponentials. The number of additions and multiplications will be linear in $cN$ as well, differently from the initial formula that would compute $c^N$ terms.
For some applications you may want the logarithm of that expression, in such situations it is preferable to use the sum of logarithms
$$log\left(P_{\mathbf{Y}}(y)\right) = \sum_{i=1}^{N}\log\left(\frac{1}{c \pi} \sum_{x\in \mathbf{X}_i} p_{\mathbf{X}}(x_i) e^{-||y_i - x_i||^2}\right)$$
All of this is valid for QAM modulation as well, and don't have the exponential number of terms in the initial formula.
PDF of the energy
In the second part you mention the PDF of the energy of the symbols
That is the probability distribution of
$$ \sum_{i=0}^{N} \sum_{x_i \in \mathbf{X}_i} p_{\mathbf{X}_i}(x_i) ||x_i + w_i||^2$$
Assuming $w$ having a variance $\sigma$ for both real and imaginary parts, since the norm of a complex number $||a + 1j b||^2 = a^2 + b^2$, assuming uniform probability for the symbols in the constelation, we can divide the whole expression the symbols $\sigma$ (since the the formula is given for a covariance matrix being the identity). And the distribution of the energy will be of the form $\sum_{k=1}^{2N} X_{i}^2$ as assumed for the noncentral $\chi^2$ distribution with $2N$ degrees of freedom.
The noncentrality parameter given by
$$\lambda =\frac{N}{c \sigma^2} \sum_{x \in \mathbf{X}} ||x||^2$$
The number of degrees of freedom $k=2N$, because we separated real and imaginary parts.
Let $z$ be the energy of one sample, taking in to account that the $\chi^2$ distribution was normalized by dividing by the variance, the PDF of $z$ will be
$$p_{||\mathbf{Y}||^2}(z) = 
\frac{1}{2\sigma} e^{\left(\frac{z}{2\sigma} + \frac{\lambda}{2}\right)}
\left(\frac{z}{\lambda \sigma}\right)^{(N-1)/2}
I_{N-1}\left( \sqrt{\frac{\lambda z}{\sigma}}\right)
$$
Note:
In spite of being a closed formula it depends on the computation of $I_{\alpha}(x)$, a modified Bessel function, that for positive integer $\alpha$ and positive $x$ reduces to a series with rational coefficients in $x$, thus
$$ I_{N-1}\left( \sqrt{\frac{\lambda z}{\sigma}}\right) = \sum_{m=0}^{\infty} \frac{\left(\sqrt{\frac{\lambda z}{\sigma}}\right)^{2m+N-1}}{m!(m+N)!} $$