Berlekamp–Massey 算法

线性递推的求解

求解序列的最短线性递推: Berlekamp-Massey 算法

算法过程

求解向量序列的线性递推

求解矩阵序列的线性递推

【模板】Berlekamp–Massey 算法

#include <cstdio>
#include <cstring>
#include <iostream>
#include <vector>

using namespace std;

typedef unsigned long long u64;

const int N = 10005, mod = 998244353;
int n, m, k, a[N];
vector<int> f;

int add(int x, int y) { return x + y < mod ? x + y : x + y - mod; }
int sub(int x, int y) { return x < y ? x + mod - y : x - y; }

int power(int a, int n) {
  int tp = 1;
  while (n) {
    if (n & 1)
      tp = 1ll * tp * a % mod;
    a = 1ll * a * a % mod, n >>= 1;
  }
  return tp;
}

vector<int> bm(int n, int *a) {
  int ld = 0, p = 0;
  vector<int> f(1, 0), g(1, 0);

  for (int i = 0; i < n; i++) {
    u64 sum = 0;
    int m = f[0];
    for (int j = 1; j <= m; j++) {
      sum += 1ll * f[j] * a[i - j];
      if (j % 18 == 0)
        sum %= mod;
    }
    int s = sum %= mod;
    if (s == a[i])
      continue;

    if (!m) {
      m = i + 1, ld = a[i], p = i;
      f[0] = i + 1;
      f.resize(m + 1, 0);
      continue;
    }

    int d = sub(a[i], s);
    int x = 1ll * d * power(ld, mod - 2) % mod;
    int len = g[0];
    vector<int> h = f;
    f.resize(i - p + len + 1);
    f[0] = max(f[0], i - p + len);
    f[i - p] = add(f[i - p], x);
    for (int j = 1; j <= len; j++)
      f[i - p + j] = sub(f[i - p + j], 1ll * x * g[j] % mod);
    ld = d, p = i, g = h;
  }

  return f;
}

void polymul(int *a, int n, int *b, int m, int *c) {
  static u64 h[N];
  for (int i = 0; i <= n; i++) {
    for (int j = 0; j <= m; j++)
      h[i + j] += 1ll * a[i] * b[j];
    if (i % 18 == 0)
      for (int i = 0; i <= n + m; i++)
        h[i] %= mod;
  }

  for (int i = 0; i <= n + m; i++)
    c[i] = h[i] % mod;
  memset(h, 0, sizeof h);
}

void polymod(int *a, int n, int *b, int m) {
  static u64 h[N];
  for (int i = n; i >= m; i--) {
    int x = a[i] = sub(a[i], h[i] % mod);
    h[i] = 0;
    for (int j = 1; j <= m; j++)
      h[i - j] += 1ll * x * b[m - j];
    if (i % 18 == 0)
      for (int j = 0; j < i; j++)
        h[j] %= mod;
  }

  for (int i = m - 1; ~i; i--)
    a[i] = sub(a[i], h[i] % mod), h[i] = 0;
}

void fpower(int n, int *f, int m, int *a) {
  if (n == 1) {
    f[1] = 1;
    return;
  }

  fpower(n >> 1, f, m, a);
  polymul(f, m - 1, f, m - 1, f);
  polymod(f, 2 * m - 2, a, m);

  if (n & 1) {
    for (int i = m; i; i--)
      f[i] = f[i - 1];
    f[0] = 0;
    polymod(f, m, a, m);
  }
}

int query(int n, int *a, vector<int> &f) {
  static int g[N], _m[N];

  int m = f[0];
  if (n < m)
    return a[n];

  _m[m] = 1;
  for (int i = 1; i <= m; i++)
    _m[m - i] = sub(0, f[i]);

  fpower(n, g, m, _m);

  int ans = 0;
  for (int i = 0; i < m; i++)
    ans = (ans + 1ll * g[i] * a[i]) % mod;
  return ans;
}

int main() {
  ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);

  cin >> n >> m;
  for (int i = 0; i < n; i++)
    cin >> a[i];

  f = bm(n, a), k = f[0];
  for (int i = 1; i <= k; i++)
    cout << f[i] << " \n"[i == k];

  cout << query(m, a, f);
}

标签：Massey,return,int,tp,算法,vector,Berlekamp,include,mod
来源： https://www.cnblogs.com/spaceswalker/p/16290721.html