AIOJ

The Berlekamp–Massey algorithm finds the shortest recurrence of a sequence. Given a sequence of length $n$ whose shortest recurrence has order $m$ , Berlekamp–Massey finds the shortest recurrence for every prefix of the sequence in $O(nm)$ time. In the worst case, $m = O(n)$ , so the worst-case complexity is $O(n^2)$ .

Definition

For a sequence $\{a_0 \dots a_{n - 1} \}$ , define a recurrence as a sequence $\{r_0\dots r_m\}$ satisfying

$\sum_{j = 0} ^ m r_j a_{i - j} = 0, \forall i \ge m$

where $r_0 = 1$ . The value $m$ is called the order of the recurrence.

The shortest recurrence of the sequence $\{a_i\}$ is the recurrence with the smallest order.

Method

For convenience, we use a slightly different set of recurrence coefficients, $\{f_0 \dots f_{m - 1}\}$ , satisfying

$a_i = \sum_{j = 0} ^ {m - 1} f_j a_{i - j - 1}, \forall i \ge m$

It is easy to see that $f_i = -r_{i + 1}$ , and that the order $m$ is the same under the two definitions.

We can construct the recurrence incrementally. Consider the terms of $\{a_i\}$ in order, and adjust the coefficients $\{f_i\}$ whenever the current recurrence gives an incorrect result. For convenience, let the shortest recurrence for the first $i$ terms be $F_i = \{f_{i, j}\}$ .

Initially, $F_0 = \{\}$ . Suppose that $F_{i - 1}$ is valid for the first $i - 1$ terms of $\{a_i\}$ . For the $i$ -th term, there are two cases:

The recurrence also gives the correct value for $a_i$ . No adjustment is needed; set $F_i = F_{i - 1}$ .
The recurrence does not give the correct value for $a_i$ . Adjust $F_{i - 1}$ to obtain a new $F_i$ .

Let $\Delta_i = a_i - \sum_{j = 0} ^ m f_{i - 1, j} a_{i - j - 1}$ be the difference between $a_i$ and the value predicted by $F_{i - 1}$ .

If this is the first adjustment of the recurrence coefficients, then $a_i$ is the first nonzero term of the sequence. In this case, set $F_i$ to $i$ zeroes, that is, to $0$ coefficients. This is clearly a valid shortest recurrence.

Otherwise, let $k$ be the number of terms of $\{a_i\}$ considered when the recurrence coefficients were last adjusted. Suppose there is a sequence $G = \{g_0 \dots g_{m' - 1}\}$ satisfying

$\sum_{j = 0} ^ {m' - 1} g_j a_{i' - j - 1} = 0, \forall i' \in [m', i)$

and $\sum_{j = 0} ^ {m' - 1} g_j a_{i - j - 1} = \Delta_i$ . Then adding $F_k$ and $G$ term by term produces a valid recurrence coefficient sequence $F_i$ .

To construct $G$ , one possible choice is

$G = \{0, 0, \dots, 0, \frac{\Delta_i}{\Delta_k}, -\frac{\Delta_i}{\Delta_k}F_{k-1}\}$

There are $i - k - 1$ leading $0$ coefficients. The final term, $-\frac{\Delta_i}{\Delta_k} F_{k-1}$ , means that every term of $F_{k-1}$ is multiplied by $-\frac{\Delta_i}{\Delta_k}$ and appended to the sequence.

Indeed, $\sum_{j = 0} ^ {m' - 1} g_j a_{i - j - 1} = \Delta_k \frac{\Delta_i}{\Delta_k} = \Delta_i$ , so this construction gives a valid $G$ . Set $F_i$ to the termwise sum of $F_k$ and $G$ .

If the recurrence is required in the original form $\{r_i\}$ , negate every $\{f_j\}$ and insert $r_0 = 1$ at the beginning.

The algorithm above has complexity $O(nm)$ when the shortest recurrence has order $m$ . Since $m = O(n)$ in the worst case, the worst-case complexity is $O(n^2)$ .

When implementing the algorithm, each adjustment only needs the recurrence coefficients $F_k$ from the previous adjustment. Therefore, if only the shortest recurrence of the whole sequence is required, it is enough to store the current coefficients and the coefficients from the previous adjustment, using $O(n)$ space.

Reference implementation

vector<int> berlekamp_massey(const vector<int> &a) {
  vector<int> v, last;  // v is the answer, 0-based, p is the module
  int k = -1, delta = 0;

  for (int i = 0; i < (int)a.size(); i++) {
    int tmp = 0;
    for (int j = 0; j < (int)v.size(); j++)
      tmp = (tmp + (long long)a[i - j - 1] * v[j]) % p;

    if (a[i] == tmp) continue;

    if (k < 0) {
      k = i;
      delta = (a[i] - tmp + p) % p;
      v = vector<int>(i + 1);

      continue;
    }

    vector<int> u = v;
    int val = (long long)(a[i] - tmp + p) * power(delta, p - 2) % p;

    if (v.size() < last.size() + i - k) v.resize(last.size() + i - k);

    (v[i - k - 1] += val) %= p;

    for (int j = 0; j < (int)last.size(); j++) {
      v[i - k + j] = (v[i - k + j] - (long long)val * last[j]) % p;
      if (v[i - k + j] < 0) v[i - k + j] += p;
    }

    if ((int)u.size() - i < (int)last.size() - k) {
      last = u;
      k = i;
      delta = a[i] - tmp;
      if (delta < 0) delta += p;
    }
  }

  for (auto &x : v) x = (p - x) % p;
  v.insert(v.begin(), 1);

  return v;  // $\forall i, \sum_{j = 0} ^ m a_{i - j} v_j = 0$
}

The naive Berlekamp–Massey algorithm finds the shortest recurrence of a finite sequence. If the sequence is infinite and an upper bound on the order of its shortest recurrence is known, taking the first $2m$ terms is enough to determine the shortest recurrence of the entire sequence. (Proof omitted.)

Applications

Because the Berlekamp–Massey algorithm has relatively poor numerical stability, it is rarely used for problems over the real numbers. For convenience, all computations below are assumed to take place in the residue field of a prime $p$ .

Finding the shortest recurrence of a vector or matrix sequence

To find the shortest recurrence of a vector sequence $\boldsymbol{v}_i$ , let the vector dimension be $n$ . Choose a random $n$ -dimensional row vector $\mathbf u^T$ , and find the shortest recurrence of the scalar sequence $\{\boldsymbol{u}^T\boldsymbol{v}_i\}$ . By the Schwartz–Zippel lemma, the two shortest recurrences are identical with probability at least $1 - \frac n p$ .

The matrix case is similar. For a matrix sequence $\{A_i\}$ of size $n \times m$ , choose a random $1 \times n$ row vector $\mathbf u^T$ and a random $m \times 1$ column vector $\boldsymbol{v}$ , then find the shortest recurrence of the scalar sequence $\{\boldsymbol{u}^T A_i \boldsymbol{v}\}$ . Again by the Schwartz–Zippel lemma, the probability that the two recurrences agree is at least $1 - \frac{n + m} p$ .

Optimizing matrix exponentiation

Let $\boldsymbol{f}_i$ be an $n$ -dimensional column vector satisfying $\boldsymbol{f}_i = A \boldsymbol{f}_{i - 1}$ . Then $\{\boldsymbol{f}_i\}$ is a vector sequence satisfying a linear recurrence of order at most $n$ . (Proof omitted.)

We can directly compute $\boldsymbol{f}_0 \dots \boldsymbol{f}_{2n - 1}$ , find the shortest recurrence of $\{\boldsymbol{f}_i\}$ using the method above, and then invoke Homogeneous linear recurrence with constant coefficients.

If $\boldsymbol{f}_m$ is required, the complexity is $O(n^3 + n\log n \log m)$ . If $A$ is sparse with only $k$ nonzero entries, this decreases to $O(nk + n\log n \log m)$ . However, since the algorithm already needs $O(nk)$ time for preprocessing, an $O(n^2 \log m)$ linear-recurrence algorithm is also acceptable when the constraints are not tight.

Finding the minimal polynomial of a matrix

The minimal polynomial of a square matrix $A$ is the polynomial $f$ of smallest degree satisfying $f(A) = 0$ .

The minimal polynomial is exactly the shortest recurrence of $\{A^i\}$ , so it can be found directly with the Berlekamp–Massey algorithm. If $A$ is an $n$ -th order square matrix, the degree of its minimal polynomial is clearly at most $n$ .

The bottleneck is obtaining $A^i$ , since multiplying matrices directly at every step would have complexity $O(n^4)$ . But when finding the shortest recurrence of a matrix sequence, we are actually finding the shortest recurrence of $\{\boldsymbol{u}^T A^i \boldsymbol{v}\}$ . Thus it is enough to compute $A^i \boldsymbol{v}$ .

If $A$ has $k$ nonzero entries, the complexity is $O(kn + n^2)$ .

Finding the determinant of a sparse matrix

If the characteristic polynomial of a square matrix $A$ is known, its constant term multiplied by $(-1)^n$ is the determinant. However, the minimal polynomial is not necessarily the characteristic polynomial.

If $A$ is multiplied by a random diagonal matrix $B$ , the minimal polynomial of $AB$ is the characteristic polynomial with probability at least $1 - \frac {2n^2 - n} p$ . It is then enough to divide by $\text{det}\;B$ .

For an $n$ -th order square matrix $A$ with $k$ nonzero entries, the complexity is $O(kn + n ^ 2)$ .

Finding the rank of a sparse matrix

For an $n\times m$ matrix $A$ , first choose a random $n\times n$ diagonal matrix $P$ and a random $m\times m$ diagonal matrix $Q$ , then find the minimal polynomial of $Q A P A^T Q$ .

There is no need to perform matrix multiplication explicitly. When finding the minimal polynomial, we only need to multiply $Q A P A^T Q$ by a vector, so we can apply these matrices to the vector one after another. The answer is the degree of the minimal polynomial after removing all factors of $x$ .

If $A$ has $k$ nonzero entries and $n \le m$ , the complexity is $O(kn + n ^ 2)$ .

Solving a sparse system of equations

Problem: Given $A \mathbf x = \mathbf b$ , where $A$ is a full-rank sparse $n \times n$ matrix and $\mathbf b$ and $\mathbf x$ are $1\times n$ column vectors, solve for $x$ in less than $n^\omega$ time. $A, \mathbf b$ are known.

Method: Clearly, $\mathbf x = A^{-1} \mathbf b$ . If we can find the shortest recurrence $\{r_0 \dots r_{m - 1}\}$ of $\{A^i \mathbf b\}$ for $i \ge 0$ , where $m \le n$ , then

$A^{-1} \mathbf b = -\frac 1 {r_{m - 1}} \sum_{i = 0} ^ {m - 2} A^i \mathbf b r_{m - 2 - i}$

(Proof omitted.)

Since $A$ is sparse, compute $\mathbf b \dots A^{2n - 1} \mathbf b$ directly from the definition of the recurrence.

If $A$ has $k$ nonzero entries, the complexity is likewise $O(kn + n^2)$ .