标签：Massey ch 数列 cdots 算法 Berlekamp mathrm 递推 define
都 2202 年了，现代 OIer 早该会会了！参考了 此博客。
引入
Berlekamp-Massey 算法，又称为 BM 算法，其可以在 \(O(n^2)\) 时间内求解一个长度为 \(n\) 的数列的最短线性递推式。
在当今 OI 界，尚没有很多 BM 算法的应用，但在一些输入的数很少的题目中，BM 能够成为发掘题目性质（找规律）的一大助力，甚至有可能直接解出答案的线性递推式，不失为一种有效的工具。
算法流程
对于数列 \({a_{1 \cdots n}}\)，我们称数列 \({r_{1 \cdots m}}\) 为其线性递推式当且仅当 \(a_i = \sum_{j=1}^m r_j \cdot a_{i-j}\) 对于任意 \(m+1 \leq i \leq n\) 均成立。若数列 \({a_{1 \cdots n}}\) 的线性递推式 \({r_{1 \cdots m}}\) 还满足 \(m\) 是所有该数列的线性递推式中最小的，则称 \({r_{1 \cdots m}}\) 为数列 \({a_{1 \cdots n}}\) 的最短线性递推式。我们的问题是，在数列 \(a\) 已知的情况下，如何求出其最短线性递推式 \(r\)。
考虑增量法，假设我们已经求出了 \(a_{1 \cdots i-1}\) 的最短线性递推式 \(r_{1 \cdots m}\)，如何求出 \(a_{1 \cdots i}\) 的最短线性递推式。
定义 \({a_{1 \cdots i-1}}\) 的最短线性递推式 \(r_{1 \cdots m}\) 为当前递推式，记递推式被更改的次数为 \(c\)，第 \(i\) 次更改后的递推式为 \(R_i\)，特别地，定义 \(R_0\) 为空，那么当前递推式应当为 \(R_c\)。
记 \(\Delta_i=a_i−\sum_{j=1}^m r_j∗a_{i−j}\)，其中 \(r_{1 \cdots m}\) 为当前递推式，显然若 \(\Delta_i = 0\)，那么当前递推式就是 \(a_{1 \cdots i}\) 的最短线性递推式。否则，我们认为 \(R_c\) 在 \(a_i\) 处出错了，定义 \(\mathrm{fail}_i\) 为 \(R_i\) 最早的出错位置，则有 \(\mathrm{fail}_c=i\)。
若 \(c=0\)，这意味着 \(a_i\) 是序列中第一个非零元素，我们可以令 \(R_{c+1}= \{ 0,0,0,...,0 \}\)，即用 \(i\) 个 \(0\) 填充线性递推式，此时由于不存在 \(j\) 使得 \(m+1 \le j \leq i\)，因此 \(R_{c+1}\) 显然为 \(a_{1 \cdots i}\) 的线性递推式，并且由于 \(a_i\) 是序列中第一个非零元素，不难证明 \(R_{c+1}\) 也是 \(a_{1 \cdots i}\) 的最短线性递推式。
否则，即 \(c>0\)，考虑 \(R_{c-1}\) 出错的位置 \(\mathrm{fail}_{c-1}\)，记 \(\mathrm{mul} = \frac{\Delta_i}{\Delta_{\mathrm{fail}_{c-1}}}\)。我们希望得到数列 \(R'=r'_{1 \cdots m'}\)，使得 \(\sum_{j=1}^{m'} r'_j \cdot a_{k-j} = 0\) 对于任意 \(m'+1≤k≤i−1\) 均成立，并且 \(\sum_{j=1}^{m'} r'_j \cdot a_{i-j} = \Delta_i\)。如果能够找到这样的数列 \(R'\)，那么令 \(R_{c+1}=R_c+R'\) 即可（其中 \(+\) 定义为各位分别相加）。
构造数列 \(R'\) 如下：\(\{0,0,0,...,0,\mathrm{mul},−\mathrm{mul} \cdot R_{c-1} \}\)，即填充 \(i−\mathrm{fail}_{c−1}−1\) 个零，然后将数列 \(\{1,−R_{c−1} \}\) 的 \(\mathrm{mul}\) 倍放在后面。容易验证其合法性，故令 \(R_{c+1}=R_c+R'\) 即可。在最坏情况下，我们可能需要对数列进行 \(O(n)\) 次修改，因此该算法的时间复杂度为 \(O(n^2)\)。
代码
以下代码实现了求解给定 \(n\) 元数列在实数域上的最短线性递推式。显然，BM 算法只需要数域中每个非零元素均存在乘法逆元即可实现，读者不妨自行实现一下在模质数意义下的 BM 算法。
Code
/*
最黯淡的一个 梦最为炽热
万千孤单焰火 让这虚构灵魂鲜活
至少在这一刻 热爱不问为何
存在为将心声响彻
*/
#include <bits/stdc++.h>
#define pii pair<int, int>
#define mp(x, y) make_pair(x, y)
#define pb push_back
#define eb emplace_back
#define fi first
#define se second
#define int long long
#define mem(x, v) memset(x, v, sizeof(x))
#define mcpy(x, y, n) memcpy(x, y, sizeof(int) * (n))
#define lob lower_bound
#define upb upper_bound
using namespace std;
inline int read() {
int x = 0, w = 1;char ch = getchar();
while (ch > '9' || ch < '0') { if (ch == '-')w = -1;ch = getchar(); }
while (ch >= '0' && ch <= '9') x = x * 10 + ch - '0', ch = getchar();
return x * w;
}
const int MN = 2e3 + 5;
const int Mod = 998244353;
const int inf = 1e9;
const double eps = 1e-8;
inline int qPow(int a, int b = Mod - 2, int ret = 1) {
while (b) {
if (b & 1) ret = ret * a % Mod;
a = a * a % Mod, b >>= 1;
}
return ret;
}
#define dbg
int N, c, fail[MN];
double val[MN], delta[MN];
vector <double> ans[MN];
signed main(void) {
N = read();
for (int i = 1; i <= N; i++)
scanf("%lf", &val[i]);
for (int i = 1; i <= N; i++) {
double tmp = val[i];
for (int j = 0; j < ans[c].size(); j++)
tmp -= ans[c][j] * val[i - j - 1];
delta[i] = tmp;
if (fabs(tmp) <= eps) continue;
fail[c] = i;
if (!c) {
ans[++c].resize(i);
continue;
}
double mul = delta[i] / delta[fail[c - 1]];
++c, ans[c].resize(i - fail[c - 2] - 1);
ans[c].pb(mul);
for (int j = 0; j < ans[c - 2].size(); j++)
ans[c].pb(ans[c - 2][j] * -mul);
if (ans[c].size() < ans[c - 1].size()) ans[c].resize(ans[c - 1].size());
for (int j = 0; j < ans[c - 1].size(); j++)
ans[c][j] += ans[c - 1][j];
}
for (int i = 0; i < ans[c].size(); i++)
printf("%.lf ", ans[c][i]);
return 0;
}
标签：Massey,ch,数列,cdots,算法,Berlekamp,mathrm,递推,define 来源： https://www.cnblogs.com/came11ia/p/16597854.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。