How to extract coefficients of a generating function using a computer? Say that as a simple example I want to extract the fourth coefficient (partitions of $4$) of $$[x^4]\,G(x)= \prod_{i=1}^4 \frac{1}{1 - x^i}$$.
What is the best way to go about it?
Is it Wolfram Alpha? If so, how is it entered?
Is it usually calculated with ad hoc code? Suggestions?
 A: Here are some ways to compute this on a computer. First of all if you have access to powerful mathematical software like Mathematica this can be done very easily. For example like this:
n = 4;
gseries = Product[1/(1 - x^i), {i, 1, 4}];
a4 = Residue[gseries/x^(n + 1), {x, 0}]

where Residue is an instruction to extract the coefficient of the $\frac{1}{x}$ term in a power-series. We can also Series expand a function about $x=0$ and simply read off the coefficient using Series[ gseries , {x, 0, n + 1} ].
This solution relies on having all these complex functionality in the programming environment, but we can also quite easily solve this without having it.

We can represent a polynomial (or power-series) $A(x) = \sum_{n\geq 0} a_n x^n$ on a computer as a sequence $\{a_0,a_1,\ldots,a_{n_{\rm max}}\}$. The product of two power-series is given by the Cauchy-product
$$A(x)B(x) = \sum_{n\geq 0} c_n x^n~~~\text{where}~~~c_n = \sum_{k=0}^na_n b_{n-k}$$
which tells us how the sequences transform under multiplication. We will use this as a basis to compute the power-series. Below is a implementation in Mathematica (but it uses no Mathematica specific functions so it can be implemented in any programming language without needing special libraries). 
(* The maximum exponent we will cover *)
nmax = 20;

(* Returns the Cauchy product of two power-series *)
CauchyProduct[an_, bn_] := Module[{cn},
   cn[m_] := Sum[an[[k]] bn[[m + 2 - k]], {k, 1, m + 1}];
   Table[cn[m], {m, 0, nmax - 1}]
];

(* The function 1 + x^i + x^2i + ...  = (1,0,..,0,1,0,..,0,1,0,..) *)
gfunc[i_] := Table[If[Mod[j - 1, i] == 0 , 1, 0], {j, 1, nmax}];

(* Expand the product Prod gfunc_i = Prod ( 1/(1-x^i)) to order nmax *)
gseries = gfunc[1];
Do[
  gseries = CauchyProduct[gseries, gfunc[i]]
 , {i, 2, 4}];

(* The coefficient of x^n *)
n = 4;
Print["an = ", gseries[[n + 1]]];

Note that in this example we only need to take nmax = 5 to get the answer.
