Berlekamp-Massey (src/Math/berlekamp_massey.hpp)

View this file on GitHub
View document part on GitHub
Last update: 2024-09-28 17:37:33+09:00
Include: #include "src/Math/berlekamp_massey.hpp"

数列$\lbrace a_i\rbrace$の最初の $n$ 項 $\lbrace a_i\rbrace_{i=0}^{n-1}$から、その数列を生成する最小の線形漸化式 $a_i = c_0 a_{i-1} + c_1 a_{i-2} + \cdots + c_{d-1}a_{i-d}$ を求める. ( $d\le n/2$ )

計算量

$O(n^2) $

Verify

天下一プログラマーコンテスト2015予選B C - 擬二等辺三角形

Required by

行列の最小多項式 (src/LinearAlgebra/MinimalPolynomial.hpp)

Verified with

Code

#pragma once
#include <vector>
// a[n] = c[0] * a[n-1] + c[1] * a[n-2] + ... + c[d-1] * a[n-d]
// return c
template <class K> std::vector<K> berlekamp_massey(const std::vector<K> &a) {
 size_t n= a.size(), d= 0, m= 0, i, j;
 if (n == 0) return {};
 std::vector<K> c(n), b(n), tmp;
 K x= 1, y, coef;
 for (c[0]= b[0]= 1, i= 0; i < n; ++i) {
  for (++m, y= a[i], j= 1; j <= d; ++j) y+= c[j] * a[i - j];
  if (y == K()) continue;
  for (tmp= c, coef= y / x, j= m; j < n; ++j) c[j]-= coef * b[j - m];
  if (2 * d <= i) d= i + 1 - d, b= tmp, x= y, m= 0;
 }
 c.resize(d + 1), c.erase(c.begin());
 for (auto &x: c) x= -x;
 return c;
}

#line 2 "src/Math/berlekamp_massey.hpp"
#include <vector>
// a[n] = c[0] * a[n-1] + c[1] * a[n-2] + ... + c[d-1] * a[n-d]
// return c
template <class K> std::vector<K> berlekamp_massey(const std::vector<K> &a) {
 size_t n= a.size(), d= 0, m= 0, i, j;
 if (n == 0) return {};
 std::vector<K> c(n), b(n), tmp;
 K x= 1, y, coef;
 for (c[0]= b[0]= 1, i= 0; i < n; ++i) {
  for (++m, y= a[i], j= 1; j <= d; ++j) y+= c[j] * a[i - j];
  if (y == K()) continue;
  for (tmp= c, coef= y / x, j= m; j < n; ++j) c[j]-= coef * b[j - m];
  if (2 * d <= i) d= i + 1 - d, b= tmp, x= y, m= 0;
 }
 c.resize(d + 1), c.erase(c.begin());
 for (auto &x: c) x= -x;
 return c;
}