Imagine a ball on the number line, each \(A_i\) represents choosing a non-negative integer \(x_i\) and moving that ball \(A_ix_i\) spaces to the right. This can be represented with the series \(\displaystyle \prod _{i = 1} ^ N \sum _{k = 0} ^\infty x^{A_ik}\). If we define \(\displaystyle S(x) := \sum _{k = 0} ^\infty x^k\), we can write our solution as:

\[\mathrm {ans} = \sum _{m = 0} ^M [x^m] \prod _{i = 1} ^ N S(x^{A_i})\]

The “prefix sum” of a series can be obtained by multiplying it by \(S(x)\), therefore:

\[\mathrm {ans} = [x^M] S(x) \prod _{i = 1} ^ N S(x^{A_i})\]

The generating function for \(S(x) = \dfrac 1 {1 - x}\). This is a well-known identity (geometric series), but can also be derived by solving the recurrence \(S(x) = 1 + xS(x)\). Therefore:

\[\mathrm {ans} = [x^M] \dfrac 1 {(1 - x) \prod_i (1 - x^{A_i})}\]

We now introduce the Bostan-Mori algorithm. This allows us to compute the form \([x^k] \dfrac {P(x)} {Q(x)}\), i.e. obtain the coefficient of \(x^k\) of a rational power series (where \(P(x)\) and \(Q(x)\) are polynomials), and runs in \(O(d \log d \log k)\) time, where \(d := \max (\deg P, \deg Q)\).

The idea for Bostan-Mori is this: we consider a transform \(\Xi _i: i \in \{ 0, 1 \}\) that extracts and reindexes the even-order or odd-order terms of a generating function:

\[\Xi_i \left(\sum_{n \geq 0} a_n x^n \right) = \sum_{m \geq 0} a_{2m + i} x^m\]

If we could obtain \(\Xi_i \left( \dfrac {P(x)} {Q(x)} \right)\) from \(\dfrac {P(x)} {Q(x)} \) in time \(t\), we could use the following relation:

\[\lbrack x^k \rbrack \dfrac {P(x)} {Q(x)} = \lbrack x^{\lfloor k/2 \rfloor} \rbrack \Xi_{k \text{ mod } 2} \left( \dfrac {P(x)} {Q(x)} \right)\]

to repeatedly halve the order of the desired term, and obtain it in time \(O(t \log k)\).

\(\Xi_i \left( \dfrac {P(x)} {Q(x)} \right)\) can be computed as follows. First multiply \(Q(-x)\) to the numerator and denominator:

\[\frac{P(x)}{Q(x)} = \frac{P(x)Q(-x)}{Q(x)Q(-x)}\]

Since \(Q(x)Q(-x)\) is an even polynomial, we may define a polynomial \(V(x)\) such that \(V(x^2) = Q(x)Q(-x)\).

Using this \(V\), we may conclude:

\[\begin{aligned} \Xi_i\left( \frac{P(x)}{Q(x)} \right) &= \Xi_i\left( \frac{P(x)Q(-x)}{V(x^2)} \right) \\ &= \frac{\Xi_i (P(x) Q(-x))}{V(x)} \end{aligned}\]

We’ve thus obtained \(\Xi_i\left( \dfrac{P(x)}{Q(x)} \right)\) in \(O(d \log d)\) time.

Therefore, \([x^k] \dfrac {P(x)} {Q(x)}\) can be obtained by the following algorithm in \(O(d \log d \log k)\) time:

Input: \(P(x), Q(x), k\)

Output: \([x^k] \dfrac {P(x)} {Q(x)}\)

\(\texttt {while } k>0:\)

\(\qquad P(x) \larr \Xi_{k \text{ mod } 2 } (P(x)Q(-x))\)

\(\qquad Q(x) \larr \Xi_{0 } (Q(x)Q(-x))\)

\(\qquad k \larr \lfloor \frac k 2 \rfloor\)

\(\texttt {end while}\)

\(\texttt {return } \dfrac {P(0)} {Q(0)}\)

See this tutorial or this tutorial for further explanation of Bostan-Mori as well as its applications to linear recurrences.

The degree of the denominator polynomial in our case is \(D := 1 + \sum _i A_i\). Obtaining the polynomial takes \(O(D \log D \log N)\) time, so our solution should run in time \(O(D \log D (\log N + \log M))\).