There's a very pretty combinatorial proof of the general identity $$\sum_{k \geq 0} \binom{n}{rk} = \frac{1}{r} \sum_{j=0}^{r-1} (1+\omega^j)^n,$$
for $\omega$ a primitive $r$th root of unity, in Benjamin, Chen, and Kindred, "Sums of Evenly Spaced Binomial Coefficients," Mathematics Magazine 83 (5), pp. 370-373, December 2010.
They show that both sides count the number of closed $n$-walks on $C_r$ beginning at vertex 0, where $C_r$ is the directed cycle on $r$ elements with the addition of a loop at each vertex, and a walk is closed if it ends where it starts.
Left-hand side: In order for an $n$-walk to be closed, it has to take $kr$ forward moves and $n-kr$ stationary moves for some $k$.
Right-hand side: The number of closed walks starting at vertex $j$ is the same regardless of the choice of $j$, and so it suffices to prove that the total number of closed $n$-walks on $C_r$ is $\sum_{j=0}^{r-1} (1+\omega^j)^n.$ For each $n$-walk with initial vertex $j$, assign each forward step a weight of $\omega^j$ and each stationary step a weight of $1$. Define the weight of an $n$-walk itself to be the product of the weights of the steps in the walk. Thus the sum of the weights of all $n$-walks starting at $j$ is $(1+\omega^j)^n$, and $\sum_{j=0}^{r-1} (1+\omega^j)^n$ gives the total weight of all $n$-walks on $C_r$. The open $n$-walks can then be partitioned into orbits such that the sum of the weights of the walks in each orbit is $0$. Thus the open $n$-walks contribute a total of $0$ to the sum $\sum_{j=0}^{r-1} (1+\omega^j)^n$. Since a closed $n$-walk has weight $1$, $\sum_{j=0}^{r-1} (1+\omega^j)^n$ must therefore give the number of closed $n$-walks on $C_r$.
They then make a slight modification of the argument above to give a combinatorial proof of
$$\sum_{k \geq 0} \binom{n}{a+rk} = \frac{1}{r} \sum_{j=0}^{r-1} \omega^{-ja}(1+\omega^j)^n,$$
where $0 \leq a < r$.
Benjamin and Scott, in "Third and Fourth Binomial Coefficients" (Fibonacci Quarterly, 49 (2), pp. 99-101, May 2011) give different combinatorial arguments for the specific cases you're asking about, $\sum_{k \geq 0} \binom{n}{3k}$ and $\sum_{k \geq 0} \binom{n}{4k}$. I prefer the more general argument above, though, so I'll just leave this one as a link and not summarize it.
I compute the number as:
22962613905485090484656664023553639680446354041773904009552854736515325227847406277133189726330125398368919292779749255468942379217261106628518627123333063707825997829062456000137755829648008974285785398012697248956323092729277672789463405208093270794180999311632479761788925921124662329907232844394066536268833781796891701120475896961582811780186955300085800543341325166104401626447256258352253576663441319799079283625404355971680808431970636650308177886780418384110991556717934409897816293912852988275811422719154702569434391547265221166310540389294622648560061463880851178273858239474974548427800576
I found this as follows:
Let $S_j(n)$ be the number of subsets of $\{1,2,\ldots,n\}$ whose sum is congruent to $j$ modulo $5$.
When $n \geq 5$, since any subset of $\{1,2,\ldots,n\}$ is the union of a unique pair of subsets, one subset of $\{1,2,\ldots,n-5\}$ and one subset of $\{n-4,n-3,n-2,n-1,n\} \equiv \{0,1,2,3,4\} \pmod 5$, we have the recurrence $$S_j(n) = \sum_{i=0}^4 S_i(n-5)S_k(5)$$ where $k$ is the solution to the congruence $i+k \equiv j \pmod 5$.
I implemented this in GAP using the following code (note: GAP uses indices starting at 1):
S:=List([1..5],i->List([1],j->Number(Combinations([1,2,3,4,5]),S->Sum(S) mod 5=i mod 5)));
The above initialises S with the boundary conditions, i.e., when $n=5$.
M:=[ [ 5, 1, 2, 3, 4 ],
[ 4, 5, 1, 2, 3 ],
[ 3, 4, 5, 1, 2 ],
[ 2, 3, 4, 5, 1 ],
[ 1, 2, 3, 4, 5 ] ];
The matrix M gives the solutions to $i+k \equiv j \pmod 5$.
for r in [2..400] do
for j in [1..5] do
S[j][r]:=Sum([1..5],i->S[i][r-1]*S[M[i][j]][1]);
od;
od;
The above just implements the recurrence.
As some quick sanity checks: (a) the code below checks that the counts sum to $2^n$ for any given $n$ (i.e. the number of subsets of $\{1,2,\ldots,n\}$).
for r in [2..400] do
s:=Sum(List([1..5],i->S[i][r]));
t:=2^(5*r);
if(s<>t) then Print(r," failed!\n"); fi;
od;
(b) we check that the answer is divisible by $2^{400}$ (as pointed out by Ross Millikan).
gap> S[5][400] mod (2^400);
0
Best Answer
Here is a calculation using something similar to the polynomial you consider. Set $\epsilon = \cos(2\pi /n)+i\sin(2\pi /n)$. For every integer $k\geq 1$, there is a polynomial factorization $$ \prod_{j=1}^{an} \left(x-\epsilon^{jk}\right) = \left(x^{n/(n,k)}-1\right)^{a(n,k)}. $$ Also, we have $$ \sum_{j=1}^{n}\epsilon^{jb}=\begin{cases}n&:n|b,\\0&n\nmid b.\end{cases} $$ So, the number of subsets $B\subseteq \{1,\ldots,an\}$ with sum divisible by $n$ is $$ \begin{align*} \frac{1}{n}\sum_{B\subseteq\{1,\ldots,an\}}\sum_{j=0}^{n-1}\epsilon^{js(B)}&=\frac{1}{n}\sum_{j=1}^{n}\prod_{k=1}^{an}\left(1+\epsilon^{jk}\right)\\ &=\lim_{x\to 1}\frac{1}{n}\sum_{j=1}^{n}\prod_{k=1}^{an}\left(x+\epsilon^{jk}\right)\\ &=\lim_{x\to 1}\frac{1}{n}\sum_{j=1}^{n}\prod_{k=1}^{an}\left(\frac{x^2-\epsilon^{2jk}}{x-\epsilon^{jk}}\right)\\ &=\frac{1}{n}\sum_{j=1}^{n} \lim_{x\to 1}\frac{\left(x^{2n/(n,2j)}-1\right)^{a(n,2j)}}{\left(x^{n/(n,j)}-1\right)^{a(n,j)}} \end{align*}. $$ The $j$-th term in the sum is $0$ if $(n,2j)>(n,j)$ (equivalently, $n/(n,j)$ is even), and is $2^{a(n,j)}$ if $(n,2j)=(n,j)$ (equivalently, $n/(n,j)$ is odd). So, if we write $n=2^km$ with $m$ odd, the number of subsets in question is $$ \begin{align*} \frac{1}{n}\sum_{\substack{j=1\\n/(n,j)\text{ odd}}}^{n} 2^{a(n,j)}=\frac{1}{n}\sum_{j=1}^m 2^{a2^k(m,j)}=\frac{1}{n}\sum_{d|m}\varphi(m/d)2^{2^kad}. \end{align*} $$ I don't know if this sum can be further simplified.