Let $M$ be a symmetric integer valued positive definite matrix with $\det(M) \geq 1$. I would like write code to compute

$$S_M= \sum_{x \in \mathbb{Z}^n} e^{-x^TMx}.$$

One option is to simply iterate over the vectors $x$ starting with ones with small coefficients and stop if things seem to be converging. Apart from the obvious problem of determining when to stop, this method is tremendously slow if $M$ is even $10$ by $10$.

It is tempting to look at the integral $\int_{x \in \mathbb{R}^n} e^{-x^TMx}\;dx$ instead but this equals $\sqrt{\frac{\pi^n}{\det(M)}}$ which is potentially a terrible approximation (it can for example be much less than $1$ where $S_M \geq 1$).

How can one compute a good approximation for $S_M$? Even an algorithm that runs $2^n$ time would be a huge improvement over what I have currently.