Berlekamp-Massey算法

Posted 2021-02-01 cjoieryl

(BM) 算法

用处

它可以用来求常系数线性递推的系数，并且可以求出最短的
求出来有什么用呢？
你可以闷声Cayley-Hamilton定理优化递推矩阵快速幂

算法简介

首先设一个数列 (f)，我们想要试出其中满足
(f_n=sum_{i=1}^{m}a_if_{n-i}(n>m))
的最小的 (m) 以及对应的系数 (a)
考虑增量法构造

1. 首先因为要求 (n>m)，所以 (m=n) 且 (a) 都为 (0) 显然是满足条件的，所以初始可以就是全 (0)
2. 假设有一个长度为 (m) 的 (a) 在 (f_{1...n-1}) 都满足条件，并且 (f_n) 不满足了
   设 (delta_i=sum_{i=1}^{m}f_{n-i}a_i-f_n)
   我们只要构造出一个长度为 (m') 最短的 (a')
   使得 (sum_{i=1}^{m'}f_{n-i}a'_i=-delta_i) 然后 (a,a') 按位相加就好了
   怎么找到呢，实际上我们之前已经存在有一些不满足条件的情况
   假设有个 (x)
   (delta_x=sum_{i=1}^{m'}f_{x-i}a'_i-f_x)
   把 (a') 向后移动 (n-x) 位，前面补 (n-x-1) 个 (0)，第 (n-x) 位搞个 (-1)
   这样得到的长度为 (m'+n-x) 的 (b) 再搞个 (frac{-delta_i}{delta_x}) 乘起来就好了
   搞出来的 (b) 显然就是我们要求的，但是可能不是最短的
   万物皆可持久化把之前所有求过的 (a) 全部记录下来
   然后又搞个 (fail_i) 表示第 (i) 个 (a) 挂了的位置
   最后弄个变量记录一下最短的就好了

代码

可能是对的
可以去 zzq 的博客里面搞个数据测一下正确性

# include <bits/stdc++.h>
using namespace std;
typedef long long ll;

const int maxn(3005);
const int mod(1e9 + 7);

inline void Inc(int &x, int y) {
    x = x + y >= mod ? x + y - mod : x + y;
}

inline void Dec(int &x, int y) {
    x = x - y < 0 ? x - y + mod : x - y;
}

inline int Add(int x, int y) {
    return x + y >= mod ? x + y - mod : x + y;
}

inline int Sub(int x, int y) {
    return x - y < 0 ? x - y + mod : x - y;
}

inline int Pow(ll x, int y) {
    register ll ret = 1;
    for (; y; y >>= 1, x = x * x % mod)
        if (y & 1) ret = ret * x % mod;
    return ret;
}

int n, f[maxn], dt[maxn], fail[maxn], cnt, inv, mn;
vector <int> cur, vc[maxn];

int main() {
    freopen("BM-in.txt", "r", stdin);
    register int i, j, l;
    scanf("%d", &n), mn = 0;
    for (i = 1; i <= n; ++i) scanf("%d", &f[i]);
    for (i = 1; i <= n; ++i) {
        dt[i] = mod - f[i], l = vc[cnt].size();
        for (j = 0; j < l; ++j) Inc(dt[i], (ll)f[i - j - 1] * vc[cnt][j] % mod);
        if (!dt[i]) continue;
        fail[cnt] = i;
        if (!cnt) {
            vc[++cnt].resize(i);
            continue;
        }
        inv = mod - (ll)dt[i] * Pow(dt[fail[mn]], mod - 2) % mod, l = vc[mn].size();
        cur.clear(), cur.resize(i - fail[mn] - 1), cur.push_back(mod - inv);
        for (j = 0; j < l; ++j) cur.push_back((ll)inv * vc[mn][j] % mod);
        if (vc[cnt].size() > cur.size()) cur.resize(vc[cnt].size());
        for (l = vc[cnt].size(), j = 0; j < l; ++j) Inc(cur[j], vc[cnt][j]);
        if (vc[cnt].size() - i < vc[mn].size() - fail[mn]) mn = cnt;
        vc[++cnt] = cur;
    }
    cout << cur.size() << endl;
    return 0;
}