\(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;
}