Find coefficient of $x^{20}$

Dividing by $x^{78}$ and substituting $x\mapsto\frac1x$ shows that $$ \left[x^{70}\right]\prod_{k=1}^{12}\left(x^k-k\right)=\left[x^8\right]\prod_{k=1}^{12}\left(1-kx^k\right) $$ and $$ \left[x^{20}\right]\prod_{k=1}^{12}\left(x^k-k\right)=\left[x^{58}\right]\prod_{k=1}^{12}\left(1-kx^k\right) $$ It is not too difficult to derive the recursion $$ \overbrace{\left[x^m\right]\prod_{k=1}^{n\vphantom{-1}}\left(1-kx^k\right)}^{a(n,m)} =\overbrace{\left[x^m\right]\prod_{k=1}^{n-1}\left(1-kx^k\right)}^{a(n-1,m)} -n\,\overbrace{\left[x^{m-n}\right]\prod_{k=1}^{n-1}\left(1-kx^k\right)}^{a(n-1,m-n)} $$ So if we define $$ a(n,m)=\left\{\begin{array}{} 0&\text{if $m\lt0$ or $m\gt\frac{n(n+1)}2$}\\ 1&\text{if $m=0$}\\ a(n-1,m)-n\,a(n-1,m-n)&\text{otherwise} \end{array}\right. $$ then $$ a(12,8)=4 $$ and $$ a(12,58)=1152600 $$ This requires less work than computing $\prod\limits_{k=1}^{12}\left(x^k-k\right)$, but it still is not simple.

Although, when doing this by hand, the table of $a(n,m)$ gets big, it is no more difficult than computing Pascal's Triangle. For example, computing $a(12,58)$ required only $75$ values to be computed; that is less than computing $12$ rows of Pascal's Triangle.

Mathematica Code

a[n_,m_]:=If[m==0,1,If[m<0||m>n(n+1)/2,0,a[n,m]=a[n-1,m]-n a[n-1,m-n]]]

Let's formalize you process, which is valid. Your polynomial is $$ \prod_{k=1}^{12} (x^k -k) $$

Now define $\mathbb B[12]$ the space of boolean vectors of length $12$. (that is the set sequences of $0$ or $1$ of length $12$).

The expansion is

$$ \sum_{b\in\mathbb B[12]} \prod_{k=1}^{12} (b_k x^k) + (1-b_k)(-k) $$

You now want to select the terms such that $$ \prod_{k=1}^{12} (b_k x^k) + (1-b_k)(-k) = c x^{20} $$

for some $c \in \mathbb R$. Note that if $b_k = 1$, then $(1-b_k)(-k)$ is $0$. So let's ignore that term for now. You're looking for the $b \in \mathbb B[12]$ such that $$ \prod_{\substack{k\in\{1..12\}\\b_k \ne 0}} b_k x^k = x^{20}. $$

Which is $$ \prod_{k=1}^{12} x^{b_k \cdot k} = x^{20}. $$

Since $x^{a} \cdot x^{b} = x^{a+b}$, you can translate that equation on the coefficients.

$$ \sum_{k=1}^{12} b_k k = 20 $$

This is a rather recursive problem. The number of solutions $b \in \mathbb B[12]$ is the number of solutions $b \in \mathbb B[11]$ that sums to $20$ (if $b_{12} = 0$) plus the number of solutions $b \in \mathbb B[11]$ that sums to $20-12 = 8$ (if $b_{12} = 1$).

\begin{align} \left\{b \in \mathbb B[12]: \sum_{k=1}^{12} b_k k = 20 \right\} &= \left\{[b,0] : b \in \mathbb B[11] \text{ and } \sum_{k=1}^{11} b_k k = 20 \right\} \\ &\cup \left\{[b,1] : b \in \mathbb B[11] \text{ and } \sum_{k=1}^{11} b_k k = 8 \right\} \end{align}

I'll let you search a solution for the recursion if it exists (I don't have the time sorry, I feel the solution cannot be simply written anyway).

Note that for $b \in \mathbb B[12]$ that is a solution, its associated coefficient is $$\prod_{\substack{k\in\{1..12\}\\b_k \ne 1}} (-k).$$

Hence sum for each $b \in \mathbb B[12]$ that is a solution to the above recursion, their associated coefficients to your final answer.

And Wolfy says the answer is

$$ x^{78} - x^{77} - 2 x^{76} - x^{75} - x^{74} + 5 x^{73} + x^{72} + 13 x^{71} + 4 x^{70} + 2 x^{68} - 8 x^{67} - 61 x^{66} - 18 x^{65} + 14 x^{64} - 181 x^{63} - 19 x^{62} - 45 x^{61} + 228 x^{60} + 111 x^{59} + 281 x^{58} + 683 x^{57} + 1248 x^{56} + 486 x^{55} + 379 x^{54} - 1668 x^{53} - 384 x^{52} + 404 x^{51} - 7810 x^{50} - 6207 x^{49} - 9633 x^{48} - 3231 x^{47} - 6517 x^{46} + 14633 x^{45} + 416 x^{44} + 7701 x^{43} + 51725 x^{42} + 61060 x^{41} + 31028 x^{40} + 99864 x^{39} + 27934 x^{38} - 20684 x^{37} - 121316 x^{36} - 158122 x^{35} - 33316 x^{34} - 434638 x^{33} - 232730 x^{32} - 472780 x^{31} - 363980 x^{30} - 546300 x^{29} + 595776 x^{28} + 371604 x^{27} + 2874060 x^{26} + 1092936 x^{25} - 85176 x^{24} + 1736472 x^{23} + 999408 x^{22} + 1959960 x^{21} + 1152600 x^{20} + 611616 x^{19} - 2386848 x^{18} - 7381296 x^{17} - 14072832 x^{16} - 22188528 x^{15} - 4414320 x^{14} + 40194720 x^{13} - 9756000 x^{12} + 1978560 x^{11} - 2756160 x^{10} - 380160 x^{9} - 7413120 x^{8} + 39346560 x^{7} - 3991680 x^{6} + 103783680 x^{5} + 39916800 x^{4} + 79833600 x^{3} - 239500800 x^{2} - 479001600 x + 479001600 $$

so robjohn is correct.

We have that $$ P(x,n) = \prod\limits_{1\, \le \,k\, \le \,n} {\left( {x^{\,k} - k} \right)} = \left. {{\partial \over {\partial y}}\prod\limits_{1\, \le \,k\, \le \,n} {\left( {x^{\,k} - y^{\,k} } \right)\;} } \right|_{\,y\, = \,1} $$

and $$ \eqalign{ & P(x,y,n) = \prod\limits_{1\, \le \,k\, \le \,n} {\left( {x^{\,k} - y^{\,k} } \right)\;} = x^{\,n\left( {n + 1} \right)/2} \prod\limits_{1\, \le \,k\, \le \,n} {\left( {1 - \left( {{y \over x}} \right)^{\,k} } \right)\;} = \cr & = x^{\,n\left( {n + 1} \right)/2} \prod\limits_{0\, \le \,k\, \le \,n - 1} {\left( {1 - {y \over x}\left( {{y \over x}} \right)^{\,k} } \right)\;} = x^{\,n\left( {n + 1} \right)/2} \left( {{y \over x};\;{y \over x}} \right)_{\,n} \cr} $$ where $\left( {a;\;z} \right)_{\,q} $ denotes the q-Pochhammer .

It is known that, in terms of q-Binomial coefficents $\left( \matrix{ m \cr k \cr} \right)_{\,q} $ the series expansion of the $P(x,y,n)$ can be written as $$ \eqalign{ & P(x,y,n) = x^{\,n\left( {n + 1} \right)/2} \left( {{y \over x};\;{y \over x}} \right)_{\,n} = \cr & = x^{\,\left( {\scriptstyle n + 1 \atop \scriptstyle 2} \right)} \sum\limits_{0\, \le \,j\, \le \,n} {\left( { - 1} \right)^{\,j} \left( \matrix{ n \cr j \cr} \right)_{\,y/x} \left( {{y \over x}} \right)^{\,j} \left( {{y \over x}} \right)^{\,\left( {\scriptstyle n \atop \scriptstyle 2} \right)} } \cr & = \sum\limits_{0\, \le \,j\, \le \,n} {\left( { - 1} \right)^{\,j} \left( \matrix{ n \cr j \cr} \right)_{\,y/x} y^{\,\left( {\scriptstyle n \atop \scriptstyle 2} \right) + j} x^{\,n - j} } \cr} $$

But the derivative of the q-Binomial is not easily expressible.

An alternative way is to rewrite $P(x,y,n)$ as $$ \eqalign{ & P(x,y,n) = \prod\limits_{1\, \le \,k\, \le \,n} {\left( {x^{\,k} - y^{\,k} } \right)\;} = y^{\,n\left( {n + 1} \right)/2} \prod\limits_{1\, \le \,k\, \le \,n} {\left( {\left( {{x \over y}} \right)^{\,k} - 1} \right)\;} = \cr & = y^{\,n\left( {n + 1} \right)/2} Q(x/y,n) \cr} $$ and in turn, express $Q(z,n)$ as product of roots of unity $$ \eqalign{ & Q(z,n) = \prod\limits_{1\, \le \,k\, \le \,n} {\left( {z^{\,k} - 1} \right)\;} = \prod\limits_{0\, \le \,k\, \le \,n - 1} {\prod\limits_{0\, \le \,j\, \le \,k - 1} {\left( {z - e^{\,i\;\,\left( {j/\left( {k + 1} \right)} \right)\,2\pi } } \right)} \;} = \cr & = \prod\limits_{0\, \le \,j\, \le \,k\, \le \,n - 1} {\left( {z - \omega _{\,n,\,j,\,k} } \right)\;} = \prod\limits_{0\, \le \,l\, \le \,h - 1} {\left( {\,z - s_{\,l} } \right)} = \sum\limits_{0\, \le \,q\, \le \,h} {c_{\,h,\,h - q} \,z^{\,q} } \quad \left| {\;h = \left( \matrix{ {n+1} \cr 2 \cr} \right)} \right. \cr} $$ The coefficients $c_{\,h,\,h - q}$ are given by the OEIS sequence A231599. There is a recurrence indicated therein, that is (putting the indices in bracket for better readability) $$c(n,k) = c(n-1, k) - c(n-1, n*(n-1)/2-k)*(-1)^n \quad |\; n > 1$$

Another interesting recurrence is provided in this post

The $c_{\,h,\,h - q}$ coefficients could also be computed by resorting to Vieta's formulas.
But in doing that we shall be careful to note that the fraction $j/(k+1)$ when reduced will provide repeated values, and consequently repeated are the $s_l$. $$ \eqalign{ & \prod\limits_{0\, \le \,l\, \le \,h - 1} {\left( {\,z - s_{\,l} } \right)} = \sum\limits_{0\, \le \,j\, \le \,h} {c_{\,h,\,h - q} \,z^{\,q} } \quad \Rightarrow \cr & \Rightarrow \quad c_{\,h,\,m} = \left( { - 1} \right)^{\,m} \sum\limits_{\,\left\{ {k_{\,0} \,,\,k_{\,1} \,,\, \cdots \,,\,k_{\,m - 1} } \right\}\,\, \subset \;\left\{ {0\,,\,1\,,\, \cdots \,,\,h - 1} \right\}\;\;} {\prod\limits_{0\, \le \,j\, \le \,m - 1} {s_{\,k_{\,j} } } } \cr} $$ After which , going back to $P(x,y,n)$ and derive wrt $y$ is quite straight.

Otherwise, it is possible to start directly from $P(x,n)$, rewritten as $$ \eqalign{ & P(x,n) = \prod\limits_{1\, \le \,k\, \le \,n} {\left( {x^{\,k} - k} \right)} = \cr & = \prod\limits_{0\, \le \,k\, \le \,n - 1} {\;\prod\limits_{0\, \le \,j\, \le \,k - 1} {\left( {x - \left( {k + 1} \right)^{\,1/\left( {k + 1} \right)} e^{\,i\;\,\left( {j/\left( {k + 1} \right)} \right)\,2\pi } } \right)} \;} \cr} $$ and applying Vieta's formulas as above.

Although algebraically interesting, I am however concious that this approach is not computationally viable for large $n$.