Define $\delta (0)=1$ and $\delta (x) = 0$ when $x \neq 0$.

Can we find a "good" method or algorithm to work out $$ S(n) = \sum_{k_1=0}^{n-1} \sum_{k_2=0}^{n-1} \cdots \sum_{k_n=0}^{n-1} \delta (\mathrm{e}^{2 \pi k_1 \mathrm{i} /n} + \mathrm{e}^{2 \pi k_2 \mathrm{i} /n} + \cdots + \mathrm{e}^{2 \pi k_n \mathrm{i} /n})$$

with any positive integer $n$?

(It's easy to count this when $n=p^l$ or $n=pq$. In particular, $S(p)=p!$ is obvious.)