arith_sum_of_squares

For $k = 1$ the result is 2 or 0. For $k = 2$ and $k = 4$ the function factors $n$ and uses Jacobi's formulas $r_2(n) = 4 \sum_{d \mid n,\, d \text{ odd}} (-1)^{(d-1)/2}$ and $r_4(n) = 8 \sum_{d \mid n,\, 4 \nmid d} d$ [Grosswald1985]. For $k = 3$ and $k = 5$ it uses the recursion $r_k(n) = \sum_{j \in \mathbb{Z}} r_{k-1}(n - j^2)$, with $O(\sqrt n)$ calls for $k - 1 \in \{2, 4\}$; $n$ is then assumed to fit in a ulong. For other $k$ it computes the coefficients of $\vartheta_3(q)^k$ up to $q^n$ with arith_sum_of_squares_vec, which requires $n$ to fit in a slong and throws otherwise.