Let $z\in \mathbb{R}^D$ and $\Sigma\in\mathbb{R}^{D\times D}$ be a positive semi-definite matrix. Then I would like to compute the following integral $$ \int_{\mathbb{R}^D}\exp\left\lbrace-\dfrac{1}{2}z^T\Sigma^{-1}z\right\rbrace\prod_{i=1}^D(z_i+a_i)^{n_i}dz$$ where $z_i$ is the $i^{th}$ component of $z$, $n_i$ is an integer power, and $a_i\in \mathbb{R}$. Is there a nice way of solving this?
For simple cases, when $D$ is low, I can do this in a nasty, complicated way by introducing the transformation $\hat{z}=\Sigma^{-\frac{1}{2}}z$ to get
$$ \int_{\mathbb{R}^D}\exp\left\lbrace-\dfrac{1}{2}\hat{z}^T\hat{z}\right\rbrace\prod_{i=1}^D(\Sigma_i^{\frac{1}{2}}\hat{z}+a_i)^{n_i}\text{det}\left(\Sigma^{-\frac{1}{2}}\right)d\hat{z}$$
where $\Sigma_i^{\frac{1}{2}}$ is the $i^{th}$ row of $\Sigma^{\frac{1}{2}}$. Then I can expanding out each polynomial using the binomial theorem, collecting like terms and getting integral into the sum of integrals of the form of
$$ \int_{\mathbb{R}^D}\exp\left\lbrace-\dfrac{1}{2}\hat{z}^T\hat{z}\right\rbrace c_k^T\hat{z}^kd\hat{z}$$
where $c_k\in \mathbb{R}^D$ and then solving these integrals one dimension at a time. This is not really a nice approach as it doesn't easily generalize for $D$ due to not having a nice form of expanding out the polynomial or knowing how to generalize the final integral above without considering each dimension individually. Any tips, tricks, or resources for being able to solve this integral?