cp-library

C++ Library for Competitive Programming

View the Project on GitHub emthrm/cp-library

:heavy_check_mark: Berlekamp–Massey algorithm
(include/emthrm/math/formal_power_series/berlekamp-massey.hpp)

任意の体上で線形回帰数列 (linear recurrence sequence) の特性多項式を求められるアルゴリズムである。

時間計算量

$O(N^2)$

仕様

名前 戻り値 要件
template <typename T>
std::vector<T> berlekamp_massey(const std::vector<T>& s); 		$S(x) = \frac{P(x)}{Q(x)} \bmod{x^n}$ を満たす最小次数の $Q(x)$？ ${\lbrack x^0 \rbrack}Q = 1$
$\mathrm{deg}(P) < \mathrm{deg}(Q) = d$ を満たす $\frac{P(x)}{Q(x)}$ を想定するならば $n \geq 2d$ を満たす。

実装

入力を $S \mathrel{:=} (S_0, S_1, \ldots)$ とおく。数列 $C$ に対して

\[C(i) \mathrel{:=} \begin{cases} \sum_{k = 0}^{\lvert C \rvert - 1} C_k S_{i - k} & (i \geq \lvert C \rvert - 1), \\ 0 & (\text{otherwise}) \end{cases}\]

とおく。任意の $i \in \mathbb{N}$ に対して $C(i) = 0$ が成り立つ $C$ を求めたい。

初期値 $C \gets (1)$ とする。$n$ 番目 (0-based) のイテレーションを考える。$\lvert C \rvert \mathrel{:=} l$ とおき、任意の $i \in \lbrace 0, 1, \ldots, n - 1 \rbrace$ に対して $C(i) = 0$ が成り立つとする。

$\Delta_n \mathrel{:=} C(n)$ とおく。$\Delta_n \neq 0$ を仮定する。

\[D(i) = \begin{cases} 0 & (i = 0, 1, \ldots, n - 1), \\ -\Delta_n & (i = n) \end{cases}\]

を満たす $D$ を求めたい。

$\Delta_f \neq 0$ を満たす $f < n$ を一つとる（$n = \min \lbrace i \in \mathbb{N} \mid S_i \neq 0 \rbrace$ のときは後述する）。$f$ 番目のイテレーションのときの $C$ を $B$ とおく。

\[B(i) = \begin{cases} 0 & (i = 0, 1, \ldots, f - 1), \\ \Delta_f & (i = f) \end{cases}\]

が成り立つ。

\[D_i \mathrel{:=} \begin{cases} 0 & (i = 0, 1, \ldots, n - f - 1), \\ -\frac{\Delta_n}{\Delta_f} B_{i - (n - f)} & (i = n - f, n - f + 1, \ldots, n - f + \lvert B \rvert - 1) \end{cases}\]

とおくと

\[\begin{split} D(i) &= \sum_{k = 0}^{n - f + \lvert B \rvert - 1} D_k C_{i - k} \\ &= \sum_{k = n - f}^{n - f + \lvert B \rvert - 1} D_k C_{i - k} \\ &= \sum_{k = 0}^{\lvert B \rvert - 1} -\frac{\Delta_n}{\Delta_f} B_k C_{i - (k + n - f)} \\ &= -\frac{\Delta_n}{\Delta_f} \sum_{k = 0}^{\lvert B \rvert - 1} B_k C_{(i - (n - f)) - k} \\ &= -\frac{\Delta_n}{\Delta_f} B(i - (n - f)) \\ &= \begin{cases} 0 & (i = 0, 1, \ldots, n - f - 1) \\ 0 & (i = n - f, n - f + 1, \ldots, n - 1) \\ -\frac{\Delta_n}{\Delta_f} \cdot \Delta_f = -\Delta_n & (i = n) \end{cases} \end{split}\]

となり、条件を満たす。

$m \mathrel{:=} \lvert D \rvert = n - f + \lvert B \rvert$ とおく。

$n = \min \lbrace S_i \neq 0 \mid i \in \mathbb{N} \rbrace$ のとき、$C^\prime_0 \neq 0$ かつ $\lvert C^\prime \rvert \leq n + 1$ を満たす任意の数列 $C^\prime$ に対して $C^\prime(n) \neq 0$ が成り立つ。$\lvert C \rvert > n + 1$ を満たすように更新しなければならない。実装では $C \gets (1, 0, 0, \ldots, 0, -\Delta_n)$ と更新するように初期値を設定している。

参考文献

Submissons

https://judge.yosupo.jp/submission/80147

Verified with

Code

#ifndef EMTHRM_MATH_FORMAL_POWER_SERIES_BERLEKAMP_MASSEY_HPP_
#define EMTHRM_MATH_FORMAL_POWER_SERIES_BERLEKAMP_MASSEY_HPP_

#include <vector>

namespace emthrm {

template <typename T>
std::vector<T> berlekamp_massey(const std::vector<T>& s) {
  const int n = s.size();
  std::vector<T> b{1}, c{1};
  b.reserve(n);
  c.reserve(n + 1);
  int m = b.size(), l = c.size(), f = -1;
  T prv_delta = 1;
  for (int i = 0; i < n; ++i) {
    T delta = s[i];
    for (int j = 1; j < l; ++j) {
      delta += c[j] * s[i - j];
    }
    if (delta == 0) continue;
    const T mul = -delta / prv_delta;
    const int shift = i - f;
    if (m + shift > l) {
      l = m + shift;
      const std::vector<T> nxt_b = c;
      c.resize(l, 0);
      for (int j = 0; j < m; ++j) {
        c[l - 1 - j] += mul * b[m - 1 - j];
      }
      b = nxt_b;
      m = b.size();
      f = i;
      prv_delta = delta;
    } else {
      for (int j = 0; j < m; ++j) {
        c[shift + j] += mul * b[j];
      }
    }
  }
  return c;
}

}  // namespace emthrm

#endif  // EMTHRM_MATH_FORMAL_POWER_SERIES_BERLEKAMP_MASSEY_HPP_
#line 1 "include/emthrm/math/formal_power_series/berlekamp-massey.hpp"



#include <vector>

namespace emthrm {

template <typename T>
std::vector<T> berlekamp_massey(const std::vector<T>& s) {
  const int n = s.size();
  std::vector<T> b{1}, c{1};
  b.reserve(n);
  c.reserve(n + 1);
  int m = b.size(), l = c.size(), f = -1;
  T prv_delta = 1;
  for (int i = 0; i < n; ++i) {
    T delta = s[i];
    for (int j = 1; j < l; ++j) {
      delta += c[j] * s[i - j];
    }
    if (delta == 0) continue;
    const T mul = -delta / prv_delta;
    const int shift = i - f;
    if (m + shift > l) {
      l = m + shift;
      const std::vector<T> nxt_b = c;
      c.resize(l, 0);
      for (int j = 0; j < m; ++j) {
        c[l - 1 - j] += mul * b[m - 1 - j];
      }
      b = nxt_b;
      m = b.size();
      f = i;
      prv_delta = delta;
    } else {
      for (int j = 0; j < m; ++j) {
        c[shift + j] += mul * b[j];
      }
    }
  }
  return c;
}

}  // namespace emthrm
Back to top page