Hi everyone!

Recently I've posted an article about continued fractions on CP-Algorithms. One very interesting detail about it is that it is possible to define a continued fraction analogue for formal power series and, due to the connection of continued fractions with the Euclidean algorithm, it would also allow to quickly compute greatest common divisor and best rational function approximation for power series.

The task of finding the minimum linear recurrence for the given starting sequence is typically solved with the Berlekamp-Massey algorithm. In this article I would like to highlight another possible approach, with the use of the extended Euclidean algorithm.

Linear recurrence interpolation

In the previous article on linear recurrences we derived that the generating function of the linear recurrence always looks like

$$$ G(x) = \frac{P(x)}{Q(x)}, $$$

where $$$P(x)$$$ and $$$Q(x)$$$ are polynomials and $$$\deg P < \deg Q$$$. In this representation, $$$Q(x) = 1 - \sum\limits_{k=1}^d a_k x^k$$$ corresponds to

$$$ F_n = \sum\limits_{k=1}^d a_k F_{n-k}. $$$

Typical competitive programming task of recovering linear recurrence is formulated as follows:

---

Find Linear Recurrence. You're given $$$F_0, \dots, F_{m}$$$. Find $$$a_1, \dots, a_d$$$ with minimum $$$d$$$ such that $$$F_n = \sum\limits_{k=1}^d a_k F_{n-k}$$$.

---

In formal power series terms it means that we're given $$$F(x) = F_0 + F_1 x + \dots + F_{m}x^{m}$$$ and we need to find $$$\frac{P(x)}{Q(x)}$$$ such that $$$ F(x) \equiv \frac{P(x)}{Q(x)} \pmod{x^{m+1}} $$$

and $$$d = \deg Q(x)$$$ is the minimum possible. In other words, what we're asked to do is to find the Padé approximant of $$$F(x)$$$.

Padé approximants

---

Given a formal power series $$$f(x) = \sum\limits_{k=0}^\infty f_k x^k$$$, its Padé approximant of order $$$[p/q]$$$ is the rational function

$$$ R(x) = \frac{a_0 + a_1 x + \dots + a_p x^p}{1+b_1 x + \dots + b_q x^q}=\frac{P(x)}{Q(x)} $$$

such that $$$f(x) \equiv R(x) \pmod{x^{p+q+1}}$$$. The Padé approximant is unique for any given $$$p$$$ and $$$q$$$ and is denoted $$$[p/q]_{f}(x) = R(x)$$$.

---

By definition, for Padé approximants generally holds

$$$ [p/q]_f \equiv [(p+k)/(q+l)]_f \pmod {x^{p+q+1}}. $$$

From this follows that the generating function of the linear recurrence is also the Padé approximant $$$[\deg P / \deg Q]_F$$$. Moreover, since

$$$ F(x) \equiv \frac{P(x)}{Q(x)} \pmod{x^{m+1}}, $$$

it also holds that $$$\frac{P(x)}{Q(x)} = [(\deg P + \alpha)/(\deg Q + \beta)]_F$$$ for any $$$\alpha + \beta \leq m - \deg P - \deg Q$$$.

Therefore, the needed fraction $$$\frac{P(x)}{Q(x)}$$$ can be found among $$$[1/m]_F, [2/(m-1)]_F,\dots,[m/1]_F$$$.

Extended Euclidean algorithm

Let's look again on the condition

$$$ F(x) \equiv \frac{P(x)}{Q(x)} \pmod{x^{m+1}}. $$$

It translates into the following Bézout's identity:

$$$ F(x) Q(x) = P(x) + x^{m+1} K(x), $$$

where $$$K(x)$$$ is a formal power series. When $$$P(x)$$$ is divisible by $$$\gcd(F(x), x^{m+1})$$$ the solution to this equation can be found with the extended Euclidean algorithm. Turns out, the extended Euclidean algorithm can also be used to enumerate all $$$[p/q]_F$$$ with $$$p+q=m$$$.

Let's formalize the extended Euclidean algorithm of $$$A(x)$$$ and $$$B(x)$$$. Starting with $$$r_{-2} = A$$$ and $$$r_{-1} = B$$$, we compute the sequence

$$$ r_{k+1} = r_{k-1} \mod r_{k} = r_{k-1} - a_k r_k. $$$

If you're familiar with the continued fractions, you could recognize, that this sequence corresponds to the representation

$$$ \frac{A(x)}{B(x)} = a_0(x) + \frac{1}{a_1(x) + \frac{1}{a_2(x) + \frac{1}{\dots}}} = [a_0(x); a_1(x), a_2(x), \dots]. $$$

Same as with rational numbers, for such sequence it is possible to define the sequence of convergents

$$$ \frac{p_k(x)}{q_k(x)} = [a_0(x); a_1(x), \dots, a_k(x)] = \frac{a_{k}(x) p_{k-1}(x) + p_{k-2}(x)}{a_{k}(x) q_{k-1}(x) + q_{k-2}(x)}, $$$

starting with

$$$ \begin{matrix} \frac{p_{-2}(x)}{q_{-2}(x)} = \frac{0}{1}, & \frac{p_{-1}(x)}{q_{-1}(x)} = \frac{1}{0}, & \frac{p_{0}(x)}{q_{0}(x)} = \frac{a_0(x)}{1},& \dots \end{matrix} $$$

Now, one can prove the following identity:

$$$ (-1)^{k}r_k(x) = q_k(x) A(x) - p_k(x) B(x). $$$

In other words, we get three sequences $$$r_k$$$, $$$p_k$$$ and $$$q_k$$$ that define a family of solutions to the Bézout's identity.

Using $$$A(x) = x^{m+1}$$$ and $$$B(x) = F(x)$$$, we conclude that

$$$ F(x) \equiv \frac{(-1)^k r_k(x)}{q_k(x)} \pmod{x^{m+1}}. $$$

Note that $$$\deg r_k$$$ monotonously decreases until $$$r_k=0$$$ while $$$\deg q_k$$$ monotonously increases until $$$q_k = \frac{B(x)}{\gcd(B(x), x^{m+1})}$$$. Moreover,

$$$\deg q_{k} - \deg q_{k-1} = \deg r_{k-1} - \deg r_k = \deg a_{k}.$$$

From this follows that if you find $$$k$$$ such that $$$\deg q_k \leq q$$$ and $$$\deg q_{k+1} > q$$$, assuming $$$p+q=m$$$, it will hold that

$$$ [p/q]_F = \frac{(-1)^k r_k}{q_k}. $$$

And, finally, the optimal linear recurrence is defined by the first $$$k$$$ such that $$$\deg r_k < \deg q_k$$$ and $$$q_k(0) \neq 0$$$. However, note that this solution assumes that we're looking for the non-degenerate linear recurrence. In other words, we say that the recurrence

$$$F_n = \sum\limits_{k=1}^d a_k F_{n-k}$$$

is valid only when $$$a_d \neq 0$$$, that is all $$$d$$$ previous values have to be used. Otherwise you'd need to deal with some degenerate cases.

Half-GCD algorithm

The notion above provides the basis to construct the $$$O(n \log^2 n)$$$ divide and conquer algorithm of computing $$$\gcd(P, Q)$$$ in polynomials and finding the minimum linear recurrence. I have a truly marvelous demonstration of this proposition that this margin is, unfortunately, too narrow to contain. I hope, I will be able to formalize the process for good and write an article about it sometime...