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.

Great thanks to nor for the proofreading and all the useful comments to make the article more accessible and rigorous.

Tl'dr.

The procedure below is essentially a formalization of the extended Euclidean algorithm done on $$$F(x)$$$ and

Unable to parse markup [type=CF_MATHJAX]

.

If you need to find the minimum linear recurrence for a given sequence

Unable to parse markup [type=CF_MATHJAX]

, do the following:

Let

Unable to parse markup [type=CF_MATHJAX]

be the generating function of the reversed $$$F$$$.

Compute the sequence of remainders

Unable to parse markup [type=CF_MATHJAX]

such that

Unable to parse markup [type=CF_MATHJAX]

,

Unable to parse markup [type=CF_MATHJAX]

and

Unable to parse markup [type=CF_MATHJAX]

Let

Unable to parse markup [type=CF_MATHJAX]

be a polynomial such that

Unable to parse markup [type=CF_MATHJAX]

.

Compute the auxiliary sequence

Unable to parse markup [type=CF_MATHJAX]

such that

Unable to parse markup [type=CF_MATHJAX]

,

Unable to parse markup [type=CF_MATHJAX]

and

Unable to parse markup [type=CF_MATHJAX]

Pick $$$k$$$ to be the first index such that $$$\deg r_k < \deg q_k$$$. Let

Unable to parse markup [type=CF_MATHJAX]

, then it also holds that

Unable to parse markup [type=CF_MATHJAX]

for any

Unable to parse markup [type=CF_MATHJAX]

and $$$d$$$ is the minimum possible. Thus,

Unable to parse markup [type=CF_MATHJAX]

divided by $$$a_0$$$ is the characteristic polynomial of the minimum linear for $$$F$$$.

More generally, one can say for such $$$k$$$ that

Unable to parse markup [type=CF_MATHJAX]

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:

---

Library Checker — 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$$$ is the minimum possible. Note that in this particular problem it is not required for

Unable to parse markup [type=CF_MATHJAX]

to be $$$0$$$, hence

Unable to parse markup [type=CF_MATHJAX]

.

In this terms, as we will see later, what the problem asks us to find is essentially 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

Unable to parse markup [type=CF_MATHJAX]

such that

Unable to parse markup [type=CF_MATHJAX]

. The Padé approximant is denoted $$$[p/q]_{f}(x) = R(x)$$$.

For normalization purposes, we will demand

Unable to parse markup [type=CF_MATHJAX]

, that is

Unable to parse markup [type=CF_MATHJAX]

, but it is acceptable to have

Unable to parse markup [type=CF_MATHJAX]

.

---

Explanation: The definition of Padé approximants requires some normalization, as you could obtain different representations of the same rational function by multiplying its numerator and denominator with the same polynomial. Most commonly used normalization is

Unable to parse markup [type=CF_MATHJAX]

. However, to guarantee

Unable to parse markup [type=CF_MATHJAX]

and to easier deal with possible trailing zeros in the recurrence, we will use the normalization

Unable to parse markup [type=CF_MATHJAX]

.

These definitions are different in terms of whether the Padé approximant for the given normalization exists. However, when $$$\deg P < \deg Q$$$, they're equivalent in a sense that when

Unable to parse markup [type=CF_MATHJAX]

defines the linear recurrence

Unable to parse markup [type=CF_MATHJAX]

the reversed polynomial

Unable to parse markup [type=CF_MATHJAX]

defines the recurrence

Unable to parse markup [type=CF_MATHJAX]

Thus to find the first kind recurrence for

Unable to parse markup [type=CF_MATHJAX]

, we could find the second kind recurrence for

Unable to parse markup [type=CF_MATHJAX]

instead.

---

Online Judge — Rational Approximation. Given $$$f(x)$$$, compute $$$p(x)$$$ and $$$q(x)$$$ of degrees at most $$$m-1$$$ and $$$n-1$$$ such that

Unable to parse markup [type=CF_MATHJAX]

---

Extended Euclidean algorithm

Let's look again on the condition

Unable to parse markup [type=CF_MATHJAX]

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$$$.

Formalizing the algorithm

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} = r_{k-2} \mod r_{k-1} = r_{k-2} - a_{k} r_{k-1}. $$$

If you're familiar with 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

Unable to parse markup [type=CF_MATHJAX]

and

Unable to parse markup [type=CF_MATHJAX]

, we conclude that $$$ F(x) q_k(x) \equiv (-1)^k r_k(x) \pmod{x^{m+1}}. $$$

Enumerating approximants

Let's estimate the degrees of $$$r_k$$$ and $$$q_k$$$ to estimate how they relate with Padé approximants.

It follows from the recurrence that

Unable to parse markup [type=CF_MATHJAX]

, therefore

Unable to parse markup [type=CF_MATHJAX]

and

Unable to parse markup [type=CF_MATHJAX]

Multiplying both the numerator and the denominator of

Unable to parse markup [type=CF_MATHJAX]

by

Unable to parse markup [type=CF_MATHJAX]

, we get the Padé approximants $$$[p/q]_F$$$ for all $$$q$$$ from $$$\deg q_k$$$ to

Unable to parse markup [type=CF_MATHJAX]

, while also maintaining the inequality

Unable to parse markup [type=CF_MATHJAX]

.

Joining it together for all $$$q_k$$$ we see that all

Unable to parse markup [type=CF_MATHJAX]

for $$$q$$$ from $$$0$$$ to $$$m$$$ are covered.

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

Unable to parse markup [type=CF_MATHJAX]

Note that the Padé approximant is not necessarily unique (for example,

Unable to parse markup [type=CF_MATHJAX]

). However, if

Unable to parse markup [type=CF_MATHJAX]

then it is also true that

Unable to parse markup [type=CF_MATHJAX]

Assuming that both approximants satisfy the definition of

Unable to parse markup [type=CF_MATHJAX]

, we can conclude that

Unable to parse markup [type=CF_MATHJAX]

and

Unable to parse markup [type=CF_MATHJAX]

otherwise the identity

Unable to parse markup [type=CF_MATHJAX]

won't hold. That being said,

Unable to parse markup [type=CF_MATHJAX]

found by the Euclidean algorithm provides the minimum possible

Unable to parse markup [type=CF_MATHJAX]

for a fixed

Unable to parse markup [type=CF_MATHJAX]

.

Finding the linear recurrence

Let

Unable to parse markup [type=CF_MATHJAX]

for some $$$P(x)$$$ and $$$Q(x)$$$ such that

Unable to parse markup [type=CF_MATHJAX]

.

If

Unable to parse markup [type=CF_MATHJAX]

then $$$\frac{P}{Q}$$$ satisfies the definition of

Unable to parse markup [type=CF_MATHJAX]

. Otherwise,

Unable to parse markup [type=CF_MATHJAX]

would provide another pair of $$$P(x)$$$ and $$$Q(x)$$$ with the same degree of $$$Q$$$ but even smaller degree of $$$P$$$.

Thus, a representation defining a minimum linear recurrence is present among the approximants generated by the Euclidean algorithm.

In this notion, the minimum recurrence is defined by the first $$$k$$$ such that $$$\deg r_k < \deg q_k$$$ and has a characteristic polynomial $$$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 library-free version, if you prefer it.

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

Unable to parse markup [type=CF_MATHJAX]

modulo $$$h(x)$$$.

---

Let

Unable to parse markup [type=CF_MATHJAX]

and $$$\frac{p_{k-1}}{q_{k-1}} = [a_0; a_1, \dots, a_{k-1}]$$$, then

Unable to parse markup [type=CF_MATHJAX]

.

Therefore, if

Unable to parse markup [type=CF_MATHJAX]

is non-constant, the inverse doesn't exist. Otherwise, the inverse is

Unable to parse markup [type=CF_MATHJAX]

.

In the problem,

Unable to parse markup [type=CF_MATHJAX]

, therefore you need to do something better than $$$O(n^2)$$$ Euclidean algorithm.