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
where $$$P(x)$$$ and $$$Q(x)$$$ are some polynomials. In particular, $$$Q(x) = 1 - \sum\limits_{k=1}^d a_k x^k$$$ corresponds to the linear recurrence
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
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
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
From this follows that the generating function of the linear recurrence is also the Padé approximant $$$[\deg P / \deg Q]_F$$$. Moreover, since
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
It translates into the following Bézout's identity:
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
If you're familiar with the continued fractions, you could recognize, that this sequence corresponds to the representation
Same as with rational numbers, for such sequence it is possible to define the sequence of convergents
starting with
Now, one can prove the following identity:
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
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,
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
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$$$. Proof by AC.
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...