How to extract coefficients of a generating function using a computer?

Question

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?

Answer 1

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.

Answer 2

The Cauchy product in R is handled with the discrete convolution function:

> (l= lapply(seq_len(4)-1, function(i) rep(c(1, integer(i)), length.out=6)))
[[1]]
[1] 1 1 1 1 1 1

[[2]]
[1] 1 0 1 0 1 0

[[3]]
[1] 1 0 0 1 0 0

[[4]]
[1] 1 0 0 0 1 0

> round(Reduce(f = function(x, y) convolve(x, y, type = "open"), x = l))

$[1] \text{0 0 0 0 1 1 2 3 5 6 6 8 8 8 6 6}\bullet \color{blue}{5}\bullet \text{3 2 1 1}$

This is consistent with the results in Wolfram Alpha:

Input:

$$\small (1+x+x^2+x^3+x^4+x^5) \times (1 + x^2 + x^4) \times (1 + x^3) \times (1 + x^4)$$

Output:

$$\small x^{16}+x^{15}+2 x^{14}+3 x^{13}+5 x^{12}+6 x^{11}+6 x^{10}+8 x^9+8 x^8+8 x^7+6 x^6+6 x^5+\large{\color{blue}{5}}\small x^4+3 x^3+2 x^2+x+1$$

---

EDIT CONTAINING THE ANSWER SOUGHT IN THE OP:

The answer I was looking for is as follows using Wolfram Alpha (not only free, but always readily accessible):

Example 1:

... Well, simply the problem in the OP - i.e. the number of partitions of an integer (in the example, $4$):

SeriesCoefficient[\prod_{i=1}^4\frac{1}{1-x^i},{x, 0, 4}]

---

Example 2:

From this presentation online:

If a fair coin is rolled 12 times, how many times will it sum to 30?

We want the 30th coefficient of

$$\begin{align}&[x^{30}]\;(x+x^2+x^3+x^4+x^5+x^6)^{12}\\[2ex] =&[x^{30}]\;x^{12}(1+x+x^2+x^3+x^4+x^5)^{12}\\[2ex] =&[x^{18}]\;(1+x+x^2+x^3+x^4+x^5)^{12}\\[2ex] =&[x^{18}]\;\left(\frac{1-x^6}{1-x}\right)^{12} \end{align}$$

SeriesCoefficient[(1 - x^6)^(12) (1 - x)^(-12),{x, 0, 18}]

Here is the alternative presented at the end of the analytical resolution of the problem in the presentation:

---

Example 3:

The GF for the Fibonacci sequence can be derived from the recursion $f_n=f_{n-1} + f_{n-2}$ with $f_0=0$ and $f_1=1:$

$$\begin{align} F(x) &= f_0 x^0 + f_1 x^1 + f_2 x^2 + \cdots \\[2 ex] &= x + f_1 x^2 + f_2 x^3 + \cdots\\ &\quad \quad+ f_0 x^2 + f_1 x^3+\cdots\\[2ex] &=x+ xF(x) +x^2 F(x) \end{align} $$

as

$$F(x)=\frac{x}{1-x-x^2}.$$

The 10th Fibonacci number is...

SeriesCoefficient[x / (1 - x - x^2), {x, 0, 10}]