So, maybe this is too programming-oriented, but the obvious way to do the computation, with access to a computer, is to count how many subsets there are of each residue mod 5, and update that progressively as you increase the value of n. You can say that there is one subset of the empty set that's 0 mod 5, and no subsets mod anything else. List these as a list of residues, 1 of 0, 0 of 1, 0 of 2, 0 of 3, and 0 of 4, which I will represent with the (space-separated) ordered tuple 0 0 0 0 1. Adding 1 to the set will give a subset with residue 1, or 0 0 0 1 1; adding 2 to the set will give subsets with residues 2 and 3, or 0 1 1 1 1; and in general, adding a new number with residue k takes us from tuple of residues T to tuple of residues T + rot_(k mod 5)(T). You could run this program with a computer no problem, you could write it in Python and have the answer in a fraction of a second after having gotten the program right.

These tuples are clearly secretly polynomials in the ring Z[x]/(x^5-1), and the T + rot_(k)(T) operation is actually multiplication of T by (x^k + 1). I guess this is the same leap of faith as the generating function in this problem, but it seems so very natural to write out the problem like this first and then look at this expression "shift by some amount, with wraparound, and add" and turn it into operations in the ring Z[x]/(x^5 - 1).

So the problem becomes "what is the units place value of (2(x+1)(x^2+1)(x^3+1)(x^4+1))^400 mod (x^5 - 1)?" And so we can evaluate the inside of that expression really fast, and we get 2 (4 + 3x + 3x^2 + 3x^3 + 3x^4) which is _extremely_ convenient for us, because this is actually 2 (1 + 3(1+x+x^(2)+x^(3)+x^(4))) = 2(1+3Z) where Z is the fifth cyclotomic polynomial. Well, of course Z is a zero divisor of our ring, and we find out that Z^2 = 5 Z (rotated copies of Z are equal to Z, and there's five of them). So we're reducing this polynomial mod Z^2 = 5 Z and then evaluating it at Z = 1, sure, great, but that's not a convenient form. But we can just evaluate a polynomial f(Z) at Z = 5 and at Z = 0, and (f(Z) % (Z^(2)-5Z)) % Z - 1) = (f(5) + 4 f(0))/5. If f_0 is the constant term, and f_1 is the Z coefficient, and f_2 the Z^2 coefficient, and so on, we are looking for f_0 + f_1 + 5 f_2 + 25 f_3 + ..., which is 1/5 f(5) + 4/5 f(0).

And since f = 2(1+3Z)^400, we can perform the evaluation before the exponentiation, to get the answer, 2^400(16^400+4)/5 = 2^2000/5 + 2^402/5, a lot of excess but not proportionally a lot more. And we know that the polynomial we get, working backwards, is 2^(400) (16^400-1)/5 Z + 2^(400), so we not only have a nice solution for the n = 5k for some k, but we can perform these multiplications for a relatively short representation of the very large solution for arbitrary n. I do have though that the sixth roots of unity don't multiply out to such a nice polynomial since six is not prime.

Does this zero divisor in the ring I'm talking about give a _natural_ connection to the roots of unity being considered here? Because the generating function for n = 4 is the fifth cyclotomic and has this convenient form?