Berlekamp-Massey算法
阿新 • • 發佈:2019-01-05
amp 以及 const sub pan https max 想要 可能
\(BM\) 算法
用處
它可以用來求常系數線性遞推的系數,並且可以求出最短的
求出來有什麽用呢?
你可以悶聲Cayley-Hamilton定理優化遞推矩陣快速冪
算法簡介
首先設一個數列 \(f\),我們想要試出其中滿足
\(f_n=\sum_{i=1}^{m}a_if_{n-i}(n>m)\)
的最小的 \(m\) 以及對應的系數 \(a\)
考慮增量法構造
- 首先因為要求 \(n>m\),所以 \(m=n\) 且 \(a\) 都為 \(0\) 顯然是滿足條件的,所以初始可以就是全 \(0\)
- 假設有一個長度為 \(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\)
這樣得到的長度為 \(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; }
Berlekamp-Massey算法