Berlekamp-Massey is a powerful tool that can knock out almost all linear recurrence problems, but it’s often explained in the context of BCH decoding in many online tutorials, making it difficult to understand in a more general sense. There’s TLE’s Codeforces blog, which contains all the core concepts of the algorithm but is a bit terse in my opinion. There’s also Massey’s paper, which is more on the mathematically rigorous side, and I have incorporated some of Massey’s proofs into the proof section of this article. The goal of this article is to explain Berlekamp-Massey using my intuition/understanding of it.

I promise this article isn’t hard, just dense. You won’t need an understanding of mathematics beyond high school, but it’s very easy to get lost if you don’t read line-by-line and just skim it. I tried my best to make it as clear as possible, but there’s only so much I can convey with words.

Table of Contents

1. Definition of a Linear Recurrence
2. Finding the $k$th Term of a Linear Recurrence
3. What is the Berlekamp-Massey algorithm?
4. An Example
5. How does our process for choosing $d$ work?
6. Which $b$ do we choose?
7. The First Step
8. Putting it all together in code…
9. Test for Understanding
10. Proofs
11. Applications in Problems
12. Problems
13. References

Definition of a Linear Recurrence

Let’s say we have some arbitrary sequence of numbers:

Now, let’s say we use the following rule to generate an infinite sequence from this initial sequence:

Using this rule, what will the next element of this sequence be? What will $s_5$ be?

And what will $s_6$ be?

Following this rule, here are the first 10 elements of our sequence:

We say the following sequence is defined by the linear recurrence $s_i = s_{i-1} + 2s_{i-2} + 5s_{i-3} - 3s_{i-4} - s_{i-5}$. The first 5 terms, $s_0, s_1, \dots, s_4$, are the base case or initial conditions of the recurrence. In general, a linear recurrence is a recurrence relation of the form:

where $c_j$ are constants, and $n$ is the length of the linear recurrence. Technically, what we defined above is a homogeneous linear recurrence. A linear recurrence could also be non-homogeneous, such as if we tack on a constant:

However, we will only be dealing with homogeneous linear recurrences in this article. After all, we can transform a $n$-length linear recurrence with an added constant into a $(n+1)$-th length homogeneous linear recurrence anyways.

Here are some more examples (and one non-example) of linear recurrences:

Once we obtain the linear recurrence, we can compute the $k$th term of a $n$ length linear recurrence in $\mathcal O(nk)$ naively. Or, if $k$ is larger (say $k \leq 10^{18}$), we can instead compute it using matrix exponentiation in $\mathcal O(n^3 \log k)$ time. However, what might be less known is an even faster algorithm that works in $\mathcal O(n \log n \log k)$. Let’s learn about it in the next section.

Finding the $k$th Term of a Linear Recurrence

This isn’t part of the Berlekamp-Massey algorithm, but it’s still useful to know. TLE’s blog covers this algorithm pretty well, but I’ll include it here for the sake of completeness.

First, for a polynomial $f(x) = \sum_{i=0}^n c_i x^i$, we define $G(f) = \sum_{i=0}^n c_i s_i$, where $s_i$ are the terms of our infinite sequence generated by our recurrence. $G$ satisfies the properties $G(f + g) = G(f) + G(g)$ for polynomials $f$ and $g$ and $G(kf) = k G(f)$ for scalar $k$. The $k$th term of our sequence is just $G(x^k)$:

This alone is kind of silly, because $G(x^k)$ requires knowledge of $s_k$ to evaluate, which is what we’re trying to find. However, let’s consider a special polynomial. For a linear recurrence $\{c_1, c_2, \dots, c_n\}$, we define the characteristic polynomial of that recurrence as $f(x) = x^n - c_1 x^{n-1} - c_2 x^{n-2} - \dots - c_{n-1} x - c_n$.

The characteristic polynomial has interesting properties when applying $G$. What is $G(f)$?

And similarly, what is $G(x f(x))$? Or $G(x^2 f(x))$? Or $G(f(x) g(x))$ for any arbitrary polynomial $g(x)$?

What is the implication of this? Consider $x^k / f(x)$. We can express any two polynomials $a(x)$ and $b(x)$ as $a(x) = b(x) q(x) + r(x)$, where $q(x)$ is the quotient after dividing $a(x)$ by $b(x)$ and $r(x)$ is the remainder. Applying this to $x^k / f(x)$, we get

We use the notation $r(x) = a(x) \bmod b(x)$. So it suffices to compute $G(x^k \bmod f(x))$! If we use binary exponentiation, we can compute this in $\mathcal O(n \log n \log k)$:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
template<typename T>
T solve(const vector<T> &c, const vector<T> &s, long long k) {
    int n = (int) c.size();
    assert(c.size() <= s.size());

    vector<T> a = n == 1 ? vector<T>{c[0]} : vector<T>{0, 1}, x{1};
    for (; k>0; k/=2) {
        if (k % 2)
            x = mul(x, a);  // mul(a, b) computes a(x) * b(x) mod f(x)
        a = mul(a, a);
    }
    x.resize(n);

    T ret = 0;
    for (int i=0; i<n; i++)
        ret += x[i] * s[i];
    return ret;
}

If you want to know how to compute $a(x) \bmod b(x)$ in $\mathcal O(n \log n)$, you can refer to this cp-algorithms article.

In practice, we usually perform mul naively in $\mathcal O(n^2)$, giving us an $\mathcal O(n^2 \log k)$ algorithm instead. This is for two reasons.

1. Berlekamp-Massey works in $\mathcal O(n^2)$ and is usually the bottleneck anyways, so the $\mathcal O(n \log n \log k)$ algorithm has limited usage.
2. The $\mathcal O(n \log n \log k)$ algorithm uses FFT, so we would have to implement FFT with arbitrary modulus if our modulus was not FFT-friendly (e.g. $10^9 + 7$).

1
2
3
4
5
6
7
8
9
10
11
12
13
vector<T> mul(const vector<T> &a, const vector<T> &b) {
    vector<T> ret(a.size() + b.size() - 1);
    // ret = a * b
    for (int i=0; i<(int)a.size(); i++)
        for (int j=0; j<(int)b.size(); j++)
            ret[i+j] += a[i] * b[j];
    // reducing ret mod f(x)
    for (int i=(int)ret.size()-1; i>=n; i--)
        for (int j=n-1; j>=0; j--)
            ret[i-j-1] += ret[i] * c[j];
    ret.resize(min((int) ret.size(), n));
    return ret;
}

Now it’s time for the meat of the article: obtaining the linear recurrence with Berlekamp-Massey.

What is the Berlekamp-Massey algorithm?

Given some finite sequence of numbers $s$, the Berlekamp-Massey algorithm can find the shortest linear recurrence that satisfies $s$. For example, if we plug in the sequence $\{0, 1, 1, 3, 5, 11, 21\}$, the algorithm will return the sequence $\{1, 2\}$, corresponding to the linear recurrence $s_i = s_{i-1} + 2s_{i-2}$.

Before I explain the algorithm, I will unfortunately have to do a quick definition/notation dump. This is mainly so that I can make the rest of the article less verbose. Don’t worry, most of the notation is very natural in my opinion.

• Sequences will be denoted using curly brackets (e.g. $\{1, 2, 3\}$ or $\{43, 123, 2, 10\}$).
• 0-based indexing will be used for the sequence, and 1-based indexing will be used for the sequence of recurrence relation coefficients. The sequence will be $s$, and the sequence of recurrence relation coefficients will be $c$. So we will have $\{s_0, s_1, s_2, \dots\}$ and $\{c_1, c_2, \dots, c_n\}$ satisfying $s_i = \sum_{j=1}^n c_j s_{i-j}$ for $i \geq n$.
• The length of a sequence is denoted as $|c|$ or $n$.
• Say I have some recurrence relation sequence $c = \{c_1, c_2, \dots, c_n\}$. I will denote “evaluating the recurrence relation at index $i$” as $c(i) = \sum_{j=1}^n c_j s_{i-j}$. If the recurrence relation correctly evaluates $s_i$, then $c(i) = s_i$.
  • A recurrence relation is correct if it successfully evaluates $c(i) = s_i$ for all $i \geq n$.
  • To avoid negative indices for $c(i)$ when $i < n$, we will simply say $c(i) = s_i$ is always true for $i < n$. This is an acceptable definition since those $s_i$ form the base case anyways.
  • If $c$ is empty, $c(i) = 0$.
• Say we have two sequences $c$ and $d$. I define adding two sequences $c$ and $d$ as $c + d = \{c_1 + d_1, c_2 + d_2, \dots, c_n + d_n\}$.
  • If $c$ and $d$ are of different lengths, we pad with zeros. So say $|c| = n, |d| = m > n$. Then $c + d = \{c_1 + d_1, c_2 + d_2, \dots, c_n + d_n, d_{n+1}, \dots, d_m\}$.
  • From this definition, we have the identity $c(i) + d(i) = (c + d)(i)$. You can prove this by explicitly writing the terms out.
• I define scalar multiplication with a sequence $c$ as $kc = \{kc_1, kc_2, \dots, kc_n\}$.
  • From this definition, we have the identity $k \cdot c(i) = (kc)(i)$. You can prove this by explicitly writing the terms out.
• The general gist of Berlekamp-Massey is as follows: we have some initial $c$. Each time we process the next element $s_i$, we will check if $c$ correctly evaluates $s_i$. If it’s correct, we keep $c$. If it’s wrong, we know $c$ isn’t the correct recurrence relation for our sequence. We will have $c(i) \neq s_i$. Then, we will adjust $c$ to some new $c’$ such that $c’(j) = s_j$ for all $j \leq i$.
  • We say a recurrence relation sequence $c$ fails at index $i$ if $i$ is the first index where $c(i) \neq 0$.

Phew! Ok, let’s start.

An Example

It is much easier to explain with an example. So let’s use the sequence $s = \{1, 2, 4, 8, 13, 20, 28, 215, 757, 2186\}$. We will begin with an empty recurrence relation sequence $c = \{\}$.

Let’s start by processing $i = 0$. We have $c(0) = 0 \neq 1 = s_0$ (recall that if $c$ is empty, $c(i) = 0$). Hmm, that’s not right. Since this is the first time we’re initializing $c$, we will merely set $c = \{1\}$. It’s very likely that we’ll have to change $c$ in the next step anyways, so the main purpose of initializing $c$ like this is simply to ensure $s_0$ is the base case. For all we care, we could set $c = \{0\}$ or $c = \{69\}$ and it wouldn’t matter.

Next, let’s go to $i = 1$. We have $c(1) = c_1 \cdot s_{1-1} = 1 \cdot 1 = 1 \neq 2 = s_1$. Ok, unfortunately the recurrence relation was wrong. We’ll have to adjust. Fortunately, there’s a simple fix to this: change to $c = \{2\}$. Now, $c(1) = 2 \cdot 1 = 2$ like we wanted.

Next, we process $i = 2$. We have $c(2) = c_1 \cdot s_{2-1} = 2 \cdot 2 = 4 = s_2$. It appears our $c$ works, so there’s no need to modify it.

Next, $i = 3$. $c(3) = c_1 \cdot s_{3-1} = 2 \cdot 4 = 8 = s_3$. This works, so we make no changes.

Next, $i = 4$. $c(4) = c_1 \cdot s_{4-1} = 2 \cdot 8 = 16 \neq 13 = s_3$. We need to adjust. What if we made $c = \{13 / 8\}$? The issue is, now $c(4) = 13$ like we wanted, but $c(3)$ is wrong! We can no longer assume the linear recurrence is just of length 1. We’ll need something more sophisticated.

Let’s be more specific about what we want. We want some sequence $c’ = c + d$ such that $c’(i)$ evaluates correctly for all $i \leq 4$. So we want some $d$ such that $d(i) = 0$ for $i < 4$ and $d(4) = s_4 - c(4) = -3$. We denote our desired value for $d(4)$ as $\Delta$.

Here’s the trick to generate such a $d$: let’s keep track of each previous version of $c$ and which index it failed on. So for example, we have $\{\}$ which failed at index 0 and $\{1\}$ which failed at index 1. Let’s consider the second version of $c$, the $\{1\}$ sequence that failed at index 1. Let’s denote the failure index as $f$. Here’s what we’ll do:

1. Set $d$ equal to that sequence.
2. Multiply the sequence by $-1$.
3. Insert a 1 on the left.
4. Multiply the sequence by $\frac{\Delta}{d(f + 1)} = \frac{-3}{1} = -3$.
5. Insert $i - f - 1 = 4 - 1 - 1 = 2$ zeros on the left.

So we have $d = \{0, 0, -3, 3\}$. I know this seems super arbitrary, but let’s evaluate $d$ and see what happens:

And $d$ isn’t defined for $i < 4$, since it’s a length 4 sequence. So I guess this $d$ works. Our new recurrence is thus $c + d = \{2\} + \{0, 0, -3, 3\} = \{2, 0, -3, 3\}$.

Let’s keep going. Our new $c$ will work for $i = 5$ and $i = 6$, but fail for $i = 7$.

We need some $d$ to add to $c$ such that $d(i) = 0$ for $i < 7$ and $d(7) = s_7 - c(7) = 174$.

Once again, let’s look at the past versions of $c$:

• $\{\}$ which failed at index 0.
• $\{1\}$ which failed at index 1.
• $\{2\}$ which failed at index 4.

This time, we’ll consider the third version of $c$, the $\{2\}$ sequence that failed at index 4 (I’ll explain which version you pick in a moment). We will apply the following seemingly arbitrary steps:

1. Set $d$ equal to our chosen sequence.
2. Multiply the sequence by $-1$.
3. Insert a 1 on the left.
4. Multiply the sequence by $\frac{\Delta}{d(f + 1)} = \frac{174}{-3} = -58$.
5. Insert $i - f - 1 = 7 - 4 - 1 = 2$ zeros on the left.

Does our $d$ work? Let’s find out.

Holy cow this actually works! So we add $d$ to our old $c$ to get our new $c$: $\{2, 0, -3, 3\} + \{0, 0, -58, 116\} = \{2, 0, -61, 119\}$.

Finally, we process $i = 8$ and $i = 9$, and find the recurrence is correct for both of them. So $c = \{2, 0, -61, 119\}$ is our final answer.

How does our process for choosing $d$ work?

Let’s examine those “arbitrary” steps closely, shall we?

1. Set $d$ equal to our chosen sequence.
2. Multiply the sequence by $-1$.
3. Insert a 1 on the left.
4. Multiply the sequence by $\frac{\Delta}{d(f + 1)}$.
5. Insert $i - f - 1$ zeros on the left.

First, we choose some older version of $c$, which we will denote as $b$ and having failed at index $f$. Notice that $b(j) = s_j$ for $j < f$ and $b(f) \neq s_f$. That’s true by the definition of a failed sequence. Let’s instead consider it as $s_j - b(j) = 0$ for $j < f$ and $s_f - b(f) \neq 0$. That’s essentially what steps 1-3 are doing: they set $d$ such that $d(j + 1) = s_j - b(j)$.

We know that $d(j) = 0$ for $j \leq f$ and $d(f + 1) \neq 0$. After multiplying the sequence by $\frac{\Delta}{d(f + 1)}$ per step 4, we get that $d(j) = 0$ for $j \leq f$ still holds true, but now $d(f + 1) = \Delta$.

Finally, step 5 makes us insert $i - f - 1$ zeros on the left. To understand the effect of that, let’s start by inserting one zero. Let $d’$ be $d$ with an extra zero inserted on the left. We have that $d’(j) = d(j - 1)$. To prove this, just write out the terms explicitly. So inserting zeros on the left is just a shift! If we insert $k$ zeros, we get $d’(j) = d(j - k)$. If we insert $i - f - 1$ zeros, we get $d’(j) = d(j - (i - f - 1))$. Specifically, $d’(i) = d(i - (i - f - 1)) = d(f + 1) = \Delta$, and $d’(j) = 0$ for $j < i$. This is exactly the condition we wanted to impose on our $d$ sequence!

Which $b$ do we choose?

While choosing any $b$ will work, let’s not forget that in addition to finding a linear recurrence that works, we want a linear recurrence of minimal length.

Take a look at the diagram below.

We are currently processing $i = 8$. The current linear recurrence sequence is the red one, is of length 4, and has just failed at index 8. We have three previous versions of $c$: a yellow sequence of length 1 that failed at index 1, a green sequence of length 2 that failed at index 4, and a blue sequence of length 3 that failed at index 6. You can also imagine there’s one more sequence: an empty sequence that failed at index 0. I didn’t draw it since it’s of length 0.

Let’s say we pick the blue sequence as our $b$. I have now drawn a purple segment on our diagram:

This purple segment represents the length of $d$. You can verify that one unit of increased length comes from the inserted 1, and one unit of increased length comes from the $i - f - 1 = 8 - 6 - 1 = 1$ zeros inserted. And the final updated $c’$ comes from overlaying the red and purple segments, representing us adding $c$ and $d$.

This tells us what the optimal $b$ is: the one whose segment has the rightmost left endpoint! So we only need to keep track of that best $c$ throughout the algorithm, not every previous version.

The First Step

Selecting a previous version of $c$ to use only works if you’re not updating $c$ for the first time. So what do you do if you do update $c$ for the first time? Let’s find out with the following example: $s = \{0, 0, 0, 0, 1, 0, 0, 2\}$. Once again, $c = \{\}$ initially.

Process $i = 0$. $c(0) = 0 = s_0$. So we keep going.

Process $i = 1$. $c(0) = 0 = s_1$. So we keep going.

This repeats until we reach $i = 4$. $c(0) = 0 \neq 1 = s_4$. Since we are updating $c$ for the first time, there are no previous versions of $c$ to rely on. So we will initialize $c$ to an arbitrary length 5 sequence. The reason why it must be at least length 5 is because that is the minimal length required to include $s_4$ in the base case. If $s_4$ was not part of the base case, then it must satisfy some recurrence $s_4 = \sum_{j=1}^n c_j s_{4-j}$. But since $s_0, s_1, s_2, s_3$ are all zero, such a recurrence would always evaluate to 0. So $s_4$ must be part of the base case. For simplicity, we will assign $c = \{0, 0, 0, 0, 0\}$, although like I stated in the previous example, the values could be arbitrary.

Now, when you update $c$ a second time when failing at index 7, you can rely on $b = \{\}$ with failure index 4. You can try by hand to confirm that this works.

Notice that if the sequence is all 0s, this algorithm will simply return $c = \{\}$, which is fine since we defined $c(i) = 0$ for an empty sequence.

Putting it all together in code…

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
template<typename T>
vector<T> berlekampMassey(const vector<T> &s) {
    vector<T> c;    // the linear recurrence sequence we are building
    vector<T> oldC; // the best previous version of c to use (the one with the rightmost left endpoint)
    int f = -1;     // the index at which the best previous version of c failed on
    for (int i=0; i<(int)s.size(); i++) {
        // evaluate c(i)
        // delta = s_i - \sum_{j=1}^n c_j s_{i-j}
        // if delta == 0, c(i) is correct
        T delta = s[i];
        for (int j=1; j<=(int)c.size(); j++)
            delta -= c[j-1] * s[i-j];   // c_j is one-indexed, so we actually need index j - 1 in the code
        if (delta == 0)
            continue;   // c(i) is correct, keep going
        // now at this point, delta != 0, so we need to adjust it
        if (f == -1) {
            // this is the first time we're updating c
            // s_i was the first non-zero element we encountered
            // we make c of length i + 1 so that s_i is part of the base case
            c.resize(i + 1);
            mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
            for (T &x : c)
                x = rng();  // just to prove that the initial values don't matter in the first step, I will set to random values
            f = i;
        } else {
            // we need to use a previous version of c to improve on this one
            // apply the 5 steps to build d
            // 1. set d equal to our chosen sequence
            vector<T> d = oldC;
            // 2. multiply the sequence by -1
            for (T &x : d)
                x = -x;
            // 3. insert a 1 on the left
            d.insert(d.begin(), 1);
            // 4. multiply the sequence by delta / d(f + 1)
            T df1 = 0;  // d(f + 1)
            for (int j=1; j<=(int)d.size(); j++)
                df1 += d[j-1] * s[f+1-j];
            assert(df1 != 0);
            T coef = delta / df1;   // storing this in outer variable so it's O(n^2) instead of O(n^2 log MOD)
            for (T &x : d)
                x *= coef;
            // 5. insert i - f - 1 zeros on the left
            vector<T> zeros(i - f - 1);
            zeros.insert(zeros.end(), d.begin(), d.end());
            d = zeros;
            // now we have our new recurrence: c + d
            vector<T> temp = c; // save the last version of c because it might have a better left endpoint
            c.resize(max(c.size(), d.size()));
            for (int j=0; j<(int)d.size(); j++)
                c[j] += d[j];
            // finally, let's consider updating oldC
            if (i - (int) temp.size() > f - (int) oldC.size()) {
                // better left endpoint, let's update!
                oldC = temp;
                f = i;
            }
        }
    }
    return c;
}

I tried to annotate the code as best as I can, using the same terminology as in the article. The method uses generic types, so you can plug in modint if you’re computing something under mod, or double, if you’re working with real numbers (warning that numerical stability is not guaranteed). If you need a working modint template, you can find mine here. The complexity of this algorithm is clearly $\mathcal O(n^2)$.

Test for Understanding

Do you really understand everything I just explained? If so, try your hand at some of these exercises! For the purpose of this exercise, assume the first step fills with $i + 1$ zeros instead of random values.

1. Find the intermediate sequences of $c$ after each step of the algorithm on $s = \{0, 2, 3, 4, 5, 6, 7, 8\}$.

   Answer

   $c$ after each $i$ is processed:

   \[\begin{align*} i = 0&: c = \{\} \\ i = 1&: c = \{0, 0\} \\ i = 2&: c = \{3/2, 0\} \\ i = 3&: c = \{3/2, -1/4\} \\ i = 4&: c = \{3/2, -1/4, -1/8\} \\ i = 5&: c = \{2, -1, 0\} \\ i = 6&: c = \{2, -1, 0\} \\ i = 7&: c = \{2, -1, 0\} \end{align*}\]

   ---

2. Find the intermediate sequences of $c$ after each step of the algorithm on $s = \{1, 8, 10, 26, 46\}$.

   Answer

   $c$ after each $i$ is processed:

   \[\begin{align*} i = 0&: c = \{0\} \\ i = 1&: c = \{8\} \\ i = 2&: c = \{8, -54\} \\ i = 3&: c = \{1, 2\} \\ i = 4&: c = \{1, 2\} \end{align*}\]

   ---

3. Find the intermediate sequences of $c$ after each step of the algorithm on $s = \{1, 3, 5, 11, 25, 59, 141, 339\}$.

   Answer

   $c$ after each $i$ is processed:

   \[\begin{align*} i = 0&: c = \{0\} \\ i = 1&: c = \{3\} \\ i = 2&: c = \{3, -4\} \\ i = 3&: c = \{1, 2\} \\ i = 4&: c = \{1, 1, 3\} \\ i = 5&: c = \{3, -1, -1\} \\ i = 6&: c = \{3, -1, -1\} \\ i = 7&: c = \{3, -1, -1\} \end{align*}\]

   ---

Proofs

I know not everyone will be interested in this section, so I’ll put them under a spoiler.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
template<typename T>
vector<T> berlekampMassey(const vector<T> &s) {
    int n = (int) s.size(), l = 0, m = 1;
    vector<T> b(n), c(n);
    T ld = b[0] = c[0] = 1;
    for (int i=0; i<n; i++, m++) {
        T d = s[i];
        for (int j=1; j<=l; j++)
            d += c[j] * s[i-j];
        if (d == 0)
            continue;
        vector<T> temp = c;
        T coef = d / ld;
        for (int j=m; j<n; j++)
            c[j] -= coef * b[j-m];
        if (2 * l <= i) {
            l = i + 1 - l;
            b = temp;
            ld = d;
            m = 0;
        }
    }
    c.resize(l + 1);
    c.erase(c.begin());
    for (T &x : c)
        x = -x;
    return c;
}

Applications in Problems

Ok, this is what you probably came for. Aside from contrived problems asking us to explicitly find the minimum length linear recurrence satisfying a sequence, what can we do with Berlekamp-Massey? The answer is that we pretty much never have to figure out an explicit recurrence ever again!

Take this Codeforces problem for example. First, let’s come up with a naive DP solution. Feel free to think for a while, then click the spoiler.

This is obviously way too slow. But we can make the following observation: the transitions are the same for all $i$. When this is the case, we can turn the $m$ into a $\log m$ factor with matrix exponentation or the polynomial method. The number of states we have per $i$ is $\mathcal O(n^2 L^2) \leq 1600$. Matrix exponentiation will be too slow, but our new method using polynomials can handle this just fine!

Ok, so we know a linear recurrence of size $\mathcal O(n^2 L^2)$ exists, but what exactly is it? Let’s ask Berlekamp-Massey! Essentially, we will run the naive DP to generate the first few terms, then plug them into Berlekamp-Massey to get the full linear recurrence. It’s magic! Here’s a submission.

By the way, how many terms do we need to plug into Berlekamp-Massey before we can be sure we got the right recurrence? For example, say we’re working with the Fibonacci sequence $\{1, 1, 2, 3, 5, \dots\}$. If we just feed $\{1, 1, 2\}$ into Berlekamp-Massey, the algorithm could return $s_i = s_{i-1} + s_{i-2}$. But it could also return $s_i = 3s_{i-1} - s_{i-2}$. Or $s_i = 4s_{i-1} - 2s_{i-2}$. How do we ensure Berlekamp-Massey gives us the right recurrence?

It turns out, for a linear recurrence of length $n$, you need to feed at least $2n$ terms into Berlekamp-Massey to guarantee getting the same or equivalent recurrence.

Now as an aside, it’s worth noting that in this problem, you need way less than $\mathcal O(n^2 L^2)$ terms. According to the editorial, there are only at most 241 states that matter, because we never reach certain combinations of states. Furthermore, in practice we don’t even need to calculate an exact bound. Since extra terms can’t hurt, it’s much easier to just generate as many terms as the problem’s time limit will allow, then plug all of them into Berlekamp-Massey and hope it gets AC. How easy is that?

Problems

While none of these problems require Berlekamp-Massey to solve, Berlekamp-Massey trivializes them significantly in my opinion.

Codeforces Round 286, Div 1 E: Mr. Kitayuta’s Gift

Atcoder Beginner Contest 198F: Cube

Codechef July Challenge 2021, PARTN01: Even Odd Partition

References

https://codeforces.com/blog/entry/61306

https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=1054260