A sum of partitions A friend of mine presented a problem I found interesting:

Compute the following: $$\sum_{n=0}^\infty\left(\prod_{k=1}^j(1+x^k)\right)[x^{mn}]$$ where $P(x)[x^n]$ denotes the $x^n$ coefficient of $P$.

This problem is meant to be solved with the assistance of a program for large $j$ (larger than 1000) and $m\ll j$. I'm interested in general, but if the case of $m~|~j$ is special, then that suffices for me.
My first thought was to brute force compute all of the coefficients, but this takes $\mathcal O(j^3)$ time, which is too slow.
My second thought was that this is essentially searching for distinct partitions. The issue is the finite upper bound, meaning that this becomes inaccurate after than $x^{j+1}$ coefficient. I'm not super familiar with partitions to know how to adjust for this without brute force like above.
My third thought was to tackle this more generally with:
$$\sum_{n=0}^\infty P(x)[x^{mn}]=\frac1m\sum_{n=0}^{m-1}P(e^{2\pi in/m})$$
But this almost completely ignores the form of $P$, so I don't feel like this is how the problem should be tackled.
How can I tackle this problem?

Update:
I just realized that when $m~|~j$ we have
$$\sum_{n=0}^\infty\left(\prod_{k=1}^j(1+x^k)\right)[x^{mn}]=\frac1m\sum_{n=0}^{m-1}\left[\prod_{k=0}^{m-1}(1+e^{2\pi ink/m})\right]^{j/m}$$
which is far fewer computations. Ruby code.
Still interested in alternative solutions.
 A: I haven't checked this, but it looks correct.
When $m\mid j$ and $\gcd(m,n)=d$ then:
$$\prod_{k=0}^{m-1}\left(1-e^{2\pi ink/m}z\right)=\left(1-z^{m/d}\right)^{d}$$
Setting $z=-1,$ you get:
$$\prod_{k=0}^{m-1}\left(1+e^{2\pi ink/m}\right)=(1-(-1)^{m/d})^{d}$$
So $$\begin{align}\frac{1}{m}\sum_{n=0}^{m-1} P\left(e^{2\pi i n/m}\right)&=\frac1m\sum_{d\mid m}\phi(m/d)\left(1-(-1)^{m/d}\right)^{jd/m}\\
&=\frac{1}m\sum_{d\mid m\text{ odd}}\phi(d)2^{j/d}
\end{align}$$
In particular, if $m$ is an odd prime, and $m\mid j$, then the value you get is:
$$\frac{1}{m}\left(2^j+(m-1)2^{j/m}\right)$$
In the general case when$m\mid j$ and $m$ might not be prime, you get an approximation:
$$\frac{1}{m}2^j + O\left(2^{j/3}\right)$$
This is not surprising, since we'd expect, for large $j,$ that approximately $\frac{1}{m}$ of the sums is a multiple of $m.$ And you might expect this or something like it to be true even when $m\not\mid j.$
A: Since the only criterion for a term to contribute to the sum is for the power to be a multiple of $m$, we can consider all exponents mod $m$ i.e. $\mathbb R[x]/(x^m-1)$. Thus if $m~|~j$ then we have
$$\sum_{n=0}^\infty\left(\prod_{k=1}^j(1+x^k)\right)[x^{mn}]=\sum_{n=0}^\infty\left(\prod_{k=0}^{m-1}(1+x^k)\right)^{j/m}[x^0]$$
This can then be interpreted as matrix exponentiation:
$$\sum_{n=0}^\infty\left(\prod_{k=1}^j(1+x^k)\right)[x^{mn}]=(A^{j/m})_{0,0}$$
where
$$A_{a,b}=\left(\prod_{k=0}^{m-1}(1+x^k)\right)[x^{(a-b)\%m}]$$
which is probably the intended solution.
