Recently I am learning generating function and some application about it, and found out this algorithm is very elegant but seems to only be well-known in Japan, also I have nothing to do during lunar new year, so I decide to write a blog to share this algorithm with CF community!

edit: Actually this algorithm have been discussed before, you can also see the other explanation in the link of this comment.

If there have some mistake or question, please comment below to let me know :)

Prerequistes

• FFT/NTT
• basic understanding of generating function
• basic understanding of linear recurrence
• knowing how to compute inv of sparse formal power series in $$$O(nk)$$$

the connection between $$$\frac{P(x)}{Q(x)}$$$ and linear recurrence

In this section, I will show the relation between $$$\frac{P(x)}{Q(x)}$$$ and linear recurrence, by proving the following lemma.

lemma: A sequence $$$(a_n)_{n\ge 0}$$$ is linearly recurrent of order $$$d$$$ $$$\Leftrightarrow$$$ its generating function is $$$\frac{P(x)}{Q(x)}$$$ where $$$deg(P(x)) \le d - 1, deg(Q(x)) = d$$$.

$$$\Rightarrow$$$: By definition, we have $$$a_i = \sum\limits_{j = 1}^{d}c_j a_{i - j}$$$ for $$$i \ge d$$$, and let $$$A(x)$$$ be the power series of sequence $$$(a_n)_{n\ge 0}$$$, i.e., $$$A(x) = \sum\limits_{i = 0}^{\infty} a_i x^i$$$

Consider construct a polynomial $$$Q(x) = 1 - \sum\limits_{i = 1}^{d}c_i x^i$$$, let $$$P(x) = Q(x)A(x)$$$, then for the $$$i$$$-th term of $$$P(x)$$$ satisfy $$$i \ge d$$$, we have

Since $$$p_i = 0$$$ for all $$$i \ge d$$$, we have $$$deg(P(x)) \le d - 1$$$, and obviously $$$deg(Q(x)) = d$$$. And since $$$q_0 = 1$$$, $$$Q(x)^{-1}$$$ exist, we can multiply $$$Q(x)^{-1}$$$ on both side of $$$P(x) = Q(x)A(x)$$$, given us $$$\frac{P(x)}{Q(x)} = A(x)$$$.

$$$\Leftarrow$$$: Similarly, consider multiply $$$Q(x)$$$ on both side of $$$\frac{P(x)}{Q(x)} = A(x)$$$, which give us $$$A(x)Q(x) = P(x)$$$, then again, consider $$$i$$$-th term of $$$P(x)$$$ for all $$$i \ge d$$$, we have

which imply sequence $$$(a_n)_{n\ge 0}$$$ is linearly recurrent of order $$$d$$$.

the algorithm

Given $$$P(x)$$$ and $$$Q(x)$$$ of the linear recurrence, Bostan-Mori try to find the $$$k$$$-th term of linear recurrence with the following inductive structure.

base case, $$$k = 0$$$

From the equation $$$P(x) = A(x)Q(x)$$$, we know $$$p_0 = a_0q_0 = a_0 = P(0)$$$, so it is just $$$P(0)$$$.

inductive step, $$$k \ge 1$$$

Consider multiply $$$\frac{P(x)}{Q(x)}$$$ with $$$Q(-x)$$$ on both denominator and numerator, which give us $$$\frac{P(x)Q(-x)}{Q(x)Q(-x)}$$$, consider the $$$t$$$-th term of $$$Q(x)Q(-x)$$$ for some odd $$$t$$$, which is $$$\sum\limits_{i + j = t \newline j \text{ is even}} q_i q_j + \sum\limits_{i + j = t \newline j \text{ is odd}} q_i (-q_j)$$$, we can see that they cancel each other, thus it is $$$0$$$, which means $$$Q(x)Q(-x)$$$ is a polynomial with only even terms, and there exist some $$$V(x^2) = Q(x)Q(-x)$$$.

Since $$$V(x^2)$$$ is a even polynomial, if $$$k$$$ is even(odd), then odd(even) terms of $$$P(x)Q(-x)$$$ is not important. Why? Let $$$U(x) = P(x)Q(-x)$$$, then the new linear recurrence will be $$$\frac{U(x)}{V(x^2)}$$$, where $$$deg(U(x)) \le 2d - 1, deg(V(x^2)) = 2d$$$, by the above lemma, this is a linear recurrence with order $$$2d$$$, and let try to derive the recurrence relation from the definition of $$$Q(x)$$$.

We can see that $$$c'$$$ is the coefficient of the new linear recurrence, and only have non-zero coefficient in even terms, so if we want to compute $$$k$$$-th term of linear recurrence, we only need to know $$$k - 2, k - 4, k - 6...$$$ terms, and the terms with another parity will be useless, so we only need even(odd) term of $$$P(x)$$$.

So let try to split $$$U(x)$$$ into two polynomial with odd/even terms, $$$U_e(x) = \sum\limits_{i = 0}^{d - 1}u_{2i} x^i, U_o = \sum\limits_{i = 0}^{d - 1}u_{2i + 1} x^i$$$, then we have $$$\frac{P(x)Q(-x)}{Q(x)Q(-x)} = \frac{U(x)}{V(x^2)} = \frac{U_e(x^2) + xU_o(x^2)}{V(x^2)}$$$.

Finally, consider case working for parity of $$$k$$$, we have

the implementation (psudo-code) and complexity analysis

Actually, it is no need to force $$$q_0 = 1$$$, we can let $$$q_0$$$ be any non-zero value, the above inductive step will still hold, which just need an little modification at base case, instead of returning $$$P(0)$$$, return $$$\frac{P(0)}{Q(0)}$$$.

Bostan-Mori(P, Q, k)
  while(k >= 1)
    U(x) = P(x)Q(-x)
    if k is even
      P(x) = U_e(x)
    else
      P(x) = U_o(x)
    V(x^2) = Q(x)Q(-x)
    Q(x) = V(x)
    k = floor(k/2)
  return P(0) / Q(0)

For each phrase, all we need to do is just FFT/NTT $$$O(1)$$$ times, and since after each phrase, $$$k$$$ will be divide by $$$2$$$, so there are $$$O(\lg k)$$$ phrases, and the time complexity is $$$O(d\lg d \lg k)$$$, where $$$d$$$ is the order the linear recurrence.

Finally, here is my implementation, though it is not very fast, it works.

compare with other algorithm

• using matrix exponentiation can compute $$$k$$$-th term of linear recurrence in $$$O(d^3 \lg k)$$$, if multiply matrix naively.

• using kitamasa method can compute $$$k$$$-th term of linear recurrence in $$$O(d^2 \lg k)$$$, or $$$O(d \lg d \lg k)$$$ using FFT/NTT.

reference

A Simple and Fast Algorithm for Computing the N-th Term of a Linearly Recurrent Sequence (original paper)

線形漸化的数列のN項目の計算 (blog written by Ryuhei Mori, one of the creator of this algorithm, written in Japanese)