General Theta Functions in Mathematica?

Given a positive-definite or negative-definite matrix $A$ of size $r times r$, one can define the corresponding theta function

$$f_{A}(q) = sum_{v in mathbb{Z}^{r}} q^{v^{t}Av} = 1 + a_{1}q + a_{2}q^{2} + cdots$$

whose coefficients count lattice points of a given square with respect to $A$. For a particular $A$, in my case $16 times 16$, I’m trying to use Mathematica to compute $f_{A}(q)$ out to some reasonable order in $q$. The following seems like the natural first attempt, but it suffers from a few issues. First, here’s my matrix:

A={{-8, -4, -4, -4, -4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, {-4, -4, -2, -2, -2, -1, 0, -1, -1, -1, 0, 0, 0, -1, -1, -1}, {-4, -2, -4, -2, -2, -1, -1, 0, -1, -1, 0, -1, -1,0, 0, -1}, {-4, -2, -2, -4, -2, -1, -1, -1, 0, -1, -1, 0, -1, 0, -1, 0}, {-4, -2, -2, -2, -4, -1, -1, -1, -1, 0, -1, -1, 0, -1, 0, 0}, {-1, -1, -1, -1, -1, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {-1, 0, -1, -1, -1, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {-1, -1, 0, -1, -1, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0}, {-1, -1, -1, 0, -1, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0}, {-1, -1, -1, -1, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0}, {-1, 0, 0, -1, -1, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0}, {-1, 0, -1, 0, -1, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0}, {-1, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0}, {-1, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0}, {-1, -1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0}, {-1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2}}

And an attempt to approximate the theta function:

Theta0[k_, q_] := Sum[q^(-{n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n13, n14, n15, n16}.A.{n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n13, n14, n15, n16}), {n1, -k, k}, {n2, -k, k}, {n3, -k, k}, {n4, -k, k}, {n5, -k, k}, {n6, -k, k}, {n7, -k, k}, {n8, -k, k}, {n9, -k, k}, {n10, -k, k}, {n11, -k, k}, {n12, -k, k}, {n13, -k, k}, {n14, -k, k}, {n15, -k, k}, {n16, -k, k}];

This is summing over the lattice points in the box $[-k,k]^{16}$ for any $k$. Now, obviously a box does not at all respect the symmetry of the problem, but in principle, for large enough $k$ I will be able to correctly recover some number of coefficients of $q$.

The issue here is that this is extremely inefficient. Even just $k=2$ wouldn’t compute letting my program run. So I’m wondering if there is a way of streamlining what I’m trying to do here? Or maybe someone can see an alternative method? I realize in a rank 16 lattice there are a ton of lattice points to deal with, but I’m hoping that can be overcome somehow!