Let's briefly recall that, for a sequence $$$a_0, a_1, \dots$$$, it's ordinary generating function (OGF) is defined as a formal power series

and its exponential generating function (EGF) is defined as

Generating functions are used because it's very simple to represent convolutions with them.

For two sequences $$$a_0, a_1, \dots$$$ and $$$b_0, b_1, \dots$$$, their convolution is defined as a sequence $$$c_0, c_1, \dots$$$ such that

Conversely, their binomial convolution is defined as a sequence $$$c_0, c_1, \dots$$$, such that

If $$$A(x)$$$ and $$$B(x)$$$ are the OGFs of $$$a_i$$$ and $$$b_j$$$ then $$$A(x) B(x)$$$ is the OGF of their convolution. Conversely, if $$$A(x)$$$ and $$$B(x)$$$ are the EGFs of $$$a_i$$$ and $$$b_j$$$, then $$$A(x) B(x)$$$ is the EGF of their binomial convolution. You can read more about it here.

We will often need to extract the coefficient near $$$x^k$$$ in the expansion of an OGF $$$A(x)$$$. The conventional notation for this is

Note that if $$$A(x)$$$ is an EGF, it will return $$$\frac{a_k}{k!}$$$ instead. To avoid this, for an EGF $$$A(x)$$$ we write

meaning that $$$a_k$$$ is the coefficient near $$$\frac{x^k}{k!}$$$ in the expansion of $$$A(x)$$$.

The concepts above naturally generalize for multivariate generating functions.

For example, we can consider a two-dimensional sequence $$$a_{ij}$$$ and its bivariate OGF $$$A(x, y)$$$ defined as

Such two-dimensional sequences are often useful, for example, when analyzing sequence of polynomials. When extracting coefficients, it is possible to extract a coefficient near specific pair of powers $$$x^i y^j$$$, as well as near specific power $$$x^i$$$. In these cases,

In other words, $$$[x^i y^j]$$$ would extract the specific coefficient, while $$$[x^i]$$$ would treat $$$A(x, y)$$$ as a power series over $$$x$$$, whose coefficients are power series over $$$y$$$, and extract the coefficient near $$$x^i$$$, which in itself is a power series over $$$y$$$.

Let's substitute two of the binomial coefficients:

Taking $$$[x^n y^{2n}]$$$ out from the summation, we need to compute

Now that we regrouped and simplified the formula a bit, let's expand $$$(1+x+y)^n$$$ as $$$((1+x)+y)^n$$$ instead of $$$((1+y)+x)^n$$$:

Two last binomial coefficients can be regrouped as a multinomial coefficient:

which indeed can be computed in $$$O(1)$$$ with $$$O(n)$$$ pre-computation.

Expanding it as a binomial, we get

As we already seen above, it rewrites as

Let's now find the generating function for $$$[x^n]{}(1+x+x^2)^n$$$:

The inner sum reduces to

From this, we get the generating function

Conveniently, we can compute first $$$n$$$ terms of it in $$$O(n)$$$ using the same recurrence as for computing first $$$n$$$ terms of $$$(1+x+x^2)^n$$$ itself. Note that the derivation of the recurrence did not rely on the fact that $$$n$$$ is a positive integer, so we can just use $$$n=-\frac{1}{2}$$$ here.

Its genfunc has a somewhat nice closed form:

See my previous blog post for details on how compute $$$f(n,m)$$$. Note that it also requires series multisection.

The answer is

See my answer on math.stackexchange for details.

The answer is

See my answer on math.stackexchange for details.

Due to the Taylor expansion formula, the left-hand side essentially says that $$$F^{(i)}(G_0) = H^{(i)}(G_0)$$$ for $$$i=0,\dots, k-1$$$.

On the other hand, we can represent $$$G(x) = G_0 + x G_1(x)$$$, and consider the Taylor expansions

The terms for $$$F$$$ and $$$H$$$ match up to $$$i=k-1$$$, and all the consequent terms are divisible by $$$x^k$$$, meaning that they're equal modulo $$$x^k$$$.

Let $$$F(x)$$$ be the generating functions of Fibonacci numbers. It is generally known that

But the expression above assumes that the sequence is infinite. In our particular example, however, the sequence is finite, so we should account for it. Let's say that $$$f_0 = a$$$ and $$$f_1 = b$$$ instead of $$$f_0=f_1=1$$$. How would it change the genfunc? The denominator won't change, as it corresponds to the characteristic polynomial of the linear recurrence. But the numerator will change to $$$c+dx$$$ such that

which is equivalent to $$$(a+bx)(1-x) \equiv c+dx \pmod{x^2}$$$, thus $$$c=a$$$ and $$$d=b-a$$$. Hence, the new genfunc would be

Now we can use $$$a = f_m$$$ and $$$b=f_{m+1}$$$ and subtact $$$F'(x)$$$ multiplied by $$$x^m$$$, so that only terms from $$$0$$$ to $$$m-1$$$ remain:

Here we used the identity $$$f_{m+1} - f_m = f_{m-1}$$$ to simplify the expression. Now, we can compute $$$H(x)$$$, the remainder of $$$F_m(x)$$$ modulo $$$(x-1)^{k+1}$$$. To do this, we can substitute $$$t=x-1$$$, thus computing

after which we can do the Taylor shift back to $$$x = 1+t$$$ to get $$$H(x)=h_0 + h_1 x + \dots + h_k x^k$$$. After that, the answer is simply

Similar to $$$F_m(x)$$$, it is possible to find the exact expression for

truncated modulo $$$t^{k+1}$$$. It would be of the form $$$\frac{1-\alpha t^{k+1} - \beta t^{k+2}}{1-(1+t)-(1+t)^2}$$$. Then substitute $$$t=x-1$$$ back and compute it modulo $$$x^{k+1}$$$ in $$$O(k)$$$.

Solution is similar to example 8, except for now we also deal with limits and convergence. However, the answer is still

because LHS and RHS converge to the same limit value for $$$x \in \mathbb C$$$. So, using similar approach to the previous problem, we need to find

This is also done with substitution $$$t = x-1$$$, thus we get

Similar to the previous example, we can find this expression exactly instead of truncating modulo $$$t^{d+1}$$$:

From this we get

From the definition above, $$$H(x)$$$ can be computed exactly in $$$O(d)$$$, as the numerator is obtained from the binomial expansion, and consequent division by the denominator can also be done in $$$O(d)$$$ because it is the exact division, and the denominator is small, so one can simply use long division algorithm. After that, when you get $$$H(x) = h_0 + h_1 x + h_2 x^2 + \dots$$$, we compute

In the later sum, to compute it in $$$O(d)$$$, we need to find all the values $$$1^d, 2^d, \dots, d^d$$$ in $$$O(d)$$$. This is doable in $$$O(d)$$$ with linear sieving by only computing the value with binary exponentiation in primes and using factorization to compute the value in all the other numbers. I have implemented the solution to test it out: #144632.