Hi everyone!
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 polynomials and $$$\deg P < \deg Q$$$. In this representation, $$$Q(x) = 1 - \sum\limits_{k=1}^d a_k x^k$$$ corresponds to
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 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 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
Therefore,
Thus, $$$\frac{(-1)^kr_k}{q_k}$$$ corresponds to $$$[p/q]_F$$$ for $$$q$$$ from $$$\deg q_k$$$ to $$$m - \deg r_k$$$, starting with $$$\deg q_0 = 0$$$ and $$$m - \deg B$$$.
Note that $$$\deg q_{k} - \deg q_{k-1} = a_k = \deg r_{k-2} - \deg r_{k-1}$$$, therefore
and
Therefore, the set of fractions $$$\frac{p_k}{q_k}$$$ covers all possible approximants for $$$q$$$ from
to
which altogether constitutes the full set of approximants.
Thus, 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
Note that it is not guaranteed that $$$q_k(0) = 1$$$. Moreover, it is possible that $$$q_k(0)=0$$$. This is due to the fact that with the Euclidean algorithm we can only guarantee the normal form of the denominator
rather than
To deal with it, one should execute the algorithm for $$$x^m F(x^{-1}) = F_{m} + F_{m-1} x + \dots + F_0 x^m$$$ and reverse the resulting $$$Q(x)$$$.
Finding the linear recurrence
In this notion, the minimum recurrence is defined by the first $$$k$$$ such that $$$\deg r_k < \deg q_k$$$.
I have implemented the $$$O(n^2)$$$ algorithm as a part of my polynomial algorithms library:
// Returns the characteristic polynomial
// of the minimum linear recurrence for
// the first d+1 elements of the sequence
poly min_rec(int d = deg()) const {
// R1 = reversed (Q(x) mod x^{d+1}), R2 = x^{d+1}
auto R1 = mod_xk(d + 1).reverse(d + 1), R2 = xk(d + 1);
auto Q1 = poly(T(1)), Q2 = poly(T(0));
while(!R2.is_zero()) {
auto [a, nR] = R1.divmod(R2); // R1 = a*R2 + nR, deg nR < deg R2
tie(R1, R2) = make_tuple(R2, nR);
tie(Q1, Q2) = make_tuple(Q2, Q1 + a * Q2);
if(R2.deg() < Q2.deg()) {
return Q2 / Q2.lead(); // guarantee that the highest coefficient is 1
}
}
assert(0);
}
You can see this Library Judge submission for further details.
Also here is a shorter version, if you prefer a library-free code.
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...
As a teaser, here's an example of the problem, that (probably) requires Half-GCD:
Library Checker — Inv of Polynomials. You're given $$$f(x)$$$ and $$$h(x)$$$. Compute $$$f^{-1}(x)$$$ modulo $$$h(x)$$$.
Let $$$\frac{f(x)}{h(x)} = [a_0; a_1, \dots, a_k]$$$ and $$$\frac{p_{k-1}}{q_{k-1}} = [a_0; a_1, \dots, a_{k-1}]$$$, then $$$f q_{k-1} - h p_{k-1} = (-1)^{k-2} r_k = (-1)^{k-2} \gcd(f, p)$$$.
Therefore, if $$$r_{k} = \gcd(f, h)$$$ is non-constant, the inverse doesn't exist. Otherwise, the inverse is $$$(-1)^{k-2} q_{k-1}(x)$$$.
In the problem, $$$n \leq 5 \cdot 10^4$$$, therefore you need to do something better than $$$O(n^2)$$$ Euclidean algorithm.
