Recovering a linear recurrence with the extended Euclidean algorithm

Revision en33, by adamant, 2022-04-08 15:24:39

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.

Tl'dr.

The procedure below is essentially a formalization of the extended Euclidean algorithm done on F(x) and xm+1.

If you need to find the minimum linear recurrence for a given sequence F0,F1,,Fm, do the following:

Let F(x)=Fm+Fm1x++F0xm be the generating function of the reversed F.

Compute the sequence of remainders r2,r1,r0,,rk such that r2=F(x), r1=xm+1 and

rk=rk2modrk1.

Let ak(x) be a polynomial such that rk=rk2akrk1.

Compute the auxiliary sequence q2,q1,q0,,qk such that q2=1, q1=0 and

qk=qk2+akqk1.

Pick k to be the first index such that degrk<degqk. Let qk(x)=a0xddk=1akxdk, then it also holds that

Fn=dk=1aka0Fnk

for any nd and d is the minimum possible. Thus, qk(x) divided by a0 is the characteristic polynomial of the minimum linear for F.

More generally, one can say for such k that

F(x)(1)krk(x)qk(x)(modxm+1).


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)=P(x)Q(x),

where P(x) and Q(x) are polynomials and degP<degQ. In this representation, Q(x)=1dk=1akxk corresponds to

Fn=dk=1akFnk.

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


Find Linear Recurrence. You're given F0,,Fm. Find a1,,ad with minimum d such that Fn=dk=1akFnk.
In formal power series terms it means that we're given F(x)=F0+F1x++Fmxm and we need to find P(x)Q(x) such that
F(x)P(x)Q(x)(modxm+1)

and d is the minimum possible. Note that in this particular problem it is not required for ad to be 0, hence degQd.

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)=k=0fkxk, its Padé approximant of order [p/q] is the rational function

R(x)=a0+a1x++apxpb0+b1x++bqxq=P(x)Q(x)

such that f(x)Q(x)P(x)(modxp+q+1). The Padé approximant is denoted [p/q]f(x)=R(x).

For normalization purposes, we will demand bq=1.


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 b0=1. However, to guarantee q=degQ and to easier deal with possible trailing zeros in the recurrence, we will use the normalization bq=1.

These definitions are different in terms of whether the Padé approximant for the given normalization exists. However, when degP<degQ, they're equivalent in a sense that when Q(x)=1b1xbqxq defines the linear recurrence

Fn=dk=1bkFnk,

the reversed polynomial xqQ(x1)=bqbq1xxq defines the recurrence

Fn=dk=1bkFn+k.
Formal proof

Thus to find the first kind recurrence for F0,F1,,Fm, we could find the second kind recurrence for Fm,Fm1,,F0 instead.


Online Judge — Rational Approximation. Given f(x), compute p(x) and q(x) of degrees at most m1 and n1 such that

f(x)q(x)p(x)0(modxm+n).

Extended Euclidean algorithm

Let's look again on the condition

F(x)Q(x)P(x)(modxm+1).

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

F(x)Q(x)=P(x)+xm+1K(x),

where K(x) is a formal power series. When P(x) is divisible by gcd(F(x),xm+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 r2=A and r1=B, we compute the sequence

rk=rk2modrk1=rk2akrk1.

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

A(x)B(x)=a0(x)+1a1(x)+1a2(x)+1=[a0(x);a1(x),a2(x),].

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

pk(x)qk(x)=[a0(x);a1(x),,ak(x)]=ak(x)pk1(x)+pk2(x)ak(x)qk1(x)+qk2(x),

starting with

p2(x)q2(x)=01,p1(x)q1(x)=10,p0(x)q0(x)=a0(x)1,

Now, one can prove the following identity:

(1)krk(x)=qk(x)A(x)pk(x)B(x).
Explanation

In other words, we get three sequences rk, pk and qk that define a family of solutions to the Bézout's identity.

Using A(x)=xm+1 and B(x)=F(x), we conclude that

F(x)qk(x)(1)krk(x)(modxm+1).

Enumerating approximants

Let's estimate the degrees of rk and qk to estimate how they relate with Padé approximants.

It follows from the recurrence that degqkdegqk1=degak=degrk2degrk1, therefore

degqk=dega0++degak

and

degrk=(m+1)dega0degakdegak+1=(mdegqk)(ak+11).

Multiplying both the numerator and the denominator of (1)krkqk by 1,x,x2,,xak+11, we get the Padé approximants [p/q]F for all q from degqk to degqk+ak+11=degqk+11, while also maintaining the inequality pmq.

Joining it together for all qk we see that all [(mq)/q]F for q from 0 to m are covered.

Thus, if you find k such that degqkq and degqk+1>q, assuming p+q=m, it will hold that

[p/q]F=(1)kxqdegqkrkxqdegqkqk.

Finding the linear recurrence

In this notion, the minimum recurrence is defined by the first k such that degrk<degqk and has a characteristic polynomial qk.

I have implemented the O(n2) 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(nlog2n) 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 f1(x) modulo h(x).


Let f(x)h(x)=[a0;a1,,ak] and pk1qk1=[a0;a1,,ak1], then fqk1hpk1=(1)k2rk=(1)k2gcd(f,p).

Therefore, if rk=gcd(f,h) is non-constant, the inverse doesn't exist. Otherwise, the inverse is (1)k2qk1(x).

In the problem, n5104, therefore you need to do something better than O(n2) Euclidean algorithm.

Tags berlekamp-massey, tutorial, linear recurrence, continued fraction

History

 
 
 
 
Revisions
 
 
  Rev. Lang. By When Δ Comment
en45 English adamant 2022-07-22 14:30:59 12
en44 English adamant 2022-07-22 14:30:30 31
en43 English adamant 2022-07-21 21:44:21 7000
en42 English adamant 2022-07-21 16:37:28 18
en41 English adamant 2022-04-08 17:18:50 140
en40 English adamant 2022-04-08 17:17:39 491
en39 English adamant 2022-04-08 17:07:37 37
en38 English adamant 2022-04-08 17:06:22 21
en37 English adamant 2022-04-08 17:04:53 773
en36 English adamant 2022-04-08 15:39:20 10
en35 English adamant 2022-04-08 15:38:39 5
en34 English adamant 2022-04-08 15:27:29 74
en33 English adamant 2022-04-08 15:24:39 22
en32 English adamant 2022-04-08 15:22:01 5090 explaining on pade approximants, b0=1 and bq=1
en31 English adamant 2022-04-07 22:03:19 253 rephrase tldr
en30 English adamant 2022-04-06 23:37:11 13
en29 English adamant 2022-04-06 22:19:33 334 problem example
en28 English adamant 2022-04-06 21:22:31 13
en27 English adamant 2022-04-06 21:21:31 105
en26 English adamant 2022-04-06 21:14:41 1082 Tldr
en25 English adamant 2022-04-06 18:50:52 34
en24 English adamant 2022-04-06 18:50:22 114
en23 English adamant 2022-04-06 17:23:32 1953
en22 English adamant 2022-04-06 17:05:46 1249
en21 English adamant 2022-04-06 16:43:26 89
en20 English adamant 2022-04-06 16:32:56 838 explanation
en19 English adamant 2022-04-06 14:57:04 216 it seems it is not necessarily unique?
en18 English adamant 2022-04-06 14:02:54 509 less continued fractions, they probably scare away people
en17 English adamant 2022-04-06 13:59:17 85
en16 English adamant 2022-04-06 13:57:34 47
en15 English adamant 2022-04-06 13:56:19 9
en14 English adamant 2022-04-06 13:53:57 1338 finally got the hang of it
en13 English adamant 2022-04-06 06:46:22 113 some edge case still need to be considered (published)
en12 English adamant 2022-04-06 06:00:12 0 (saved to drafts)
en11 English adamant 2022-04-06 03:49:58 15
en10 English adamant 2022-04-06 03:41:10 345
en9 English adamant 2022-04-06 03:23:06 77
en8 English adamant 2022-04-06 03:14:13 5
en7 English adamant 2022-04-06 02:13:21 86
en6 English adamant 2022-04-06 01:37:01 13
en5 English adamant 2022-04-06 01:28:10 57
en4 English adamant 2022-04-06 01:23:00 29
en3 English adamant 2022-04-06 00:10:26 9 I'm not sure
en2 English adamant 2022-04-05 23:29:13 22
en1 English adamant 2022-04-05 23:27:31 6370 Initial revision (published)