$p(\|x\|_2 < R)$ when $x$ follows a normal distribution

A friend of mine recently asked

What is the value of $p(|x|_{2}<R)$ where $x\sim\mathcal{N}(0,\Sigma)$ and $x\in\mathbb{R}^{D}$?

The question is not as easy as it looks. Although I did not get the answer, I arrived at some expression which is close to the goal. Before I show what I did, here is one result we will need.

Here is what I did. Note that $p(|x|_{2}<R)=p(x^{\top}x<R^2)$. The goal here is to find the CDF of $x^{\top}x$. Start by eigen-decomposing $\Sigma=UVU^{\top}$. Let $y:=U^{\top}x$. It follows that $x^{\top}x=y^{\top}y=\sum_{i=1}^{D}y_{i}^{2}$ because $U$ is an orthogonal matrix. So, $p(x^{\top}x<R^2)=p(y^{\top}y<R^2)$ . From the property of the normal distribution, we know that $y=U^{\top}x\sim\mathcal{N}(0,V)$. Equivalently, $y_{i}\sim\mathcal{N}(0,v_{i})$ where $v_{i}:=V_{ii}$. Let $z\sim\mathcal{N}(0,1)$ be the standard normal random variable. We can rewrite $y_{i}$ as $y_{i}=\sqrt{v_{i}}z$, and so $y_{i}^{2}=v_{i}z^{2}$, where $z^{2}\sim\chi^{2}(1)$, the Chi-squared distribution with one degree of freedom. From the previous result, we have $y_{i}^{2}\sim\mathrm{Gamma}(\frac{1}{2},2v_{i})$. This implies that the CDF of $y^{\top}y$ is the same as the CDF of a sum of independent Gamma random variables. According to Covo and Elalouf, 2014 and Moschopoulos, 1985, the CDF is not trivial to compute.