Berlekamp-Massey算法学习小记

介绍：

BM算法的作用是求一个序列的最短递推式。

设原序列为$a$，即求出一个$f$，$f[1-len]$有值。

满足：$a[i]={\sum }_{j=1}^{len}a[i-j]\ast f[j](i>len)$

$len$最小的$f$即为所求。

BM算法本质上是一个构造算法，证明我并不会。

---

算法流程：

假设现在已经求出了对于序列$a$的前$i-1$项的最短递推式$f$

现在考虑第$i$项。

设$delta={\sum }_{j=1}^{len}a[i-j]\ast f[j]-a[i]$

若$delta=0$，说明$f$可以继续适用于序列前$i$项，则不用找了。

不然的话，我们可以构造一个递推式$f^{\prime }$，使得$f^{\prime }$代进前$i-1$项都是0，代入第$i$项是1。

那么对于前$i$项的递推式就是$f-f^{\prime }\ast delta$

我们记录下历史中失配的那些递推式，先随便取一个递推式，假设失配的位置是$fail$，这个递推式是$g$。

假设在$g$的前面添加了$i-fail$个0，设为$g^{\prime }$，即$g^{\prime }[x]=g[x-(i-fail)]$。

发现把$g^{\prime }$代入$i$位置，会得到把$g$代入$fail$位置一样的结果，代入$i-1$位置，会得到和$fail-1$一样的结果，……

此时不妨让$g^{\prime }[i-fail]-=1$，这时代入$i-1$及更小的位置就是0了，而代入$i$位置会得到和先前$g$得到的$delt{a}^{\prime }$一样的值。

那么$f^{\prime }=g^{\prime }/delt{a}^{\prime }\ast delta$

最后一个问题，我们这样只是随便构造了一组解，并不一定是最短的。

这个很简单，我们只要使g’的长度最短就好了，可以感性理解。

注意g’的长度是$i-fail+len(g)$，所以在每一次失配时，比较$fail-len(now)$，更新一下即可。

直接存也许需要$O(n^2)$的空间，所以要滚动。

---

模板：

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, B = y; i <= B; i ++)
#define ff(i, x, y) for(int i = x, B = y; i <  B; i ++)
#define fd(i, x, y) for(int i = x, B = y; i >= B; i --)
#define ll long long
#define pp printf
#define hh pp("\n")
#define db double
using namespace std;

const int N = 1e5 + 5;

const db eps = 1e-8;

int n;
db a[N];
db f[3][N], dt[3]; int len[3], fail[3];

void cp(int x, int y) {
	fo(i, 1, len[x]) f[y][i] = f[x][i];
	len[y] = len[x]; dt[y] = dt[x]; fail[y] = fail[x];
}

void solve(db *a, int n) {
	fo(i, 1, n) f[0][i] = 0; len[0] = 0;
	fo(i, 1, n) {
		db tmp = -a[i];
		fo(j, 1, len[0]) tmp += a[i - j] * f[0][j];
		dt[0] = tmp;
		if(abs(tmp) < eps) continue;
		fail[0] = i;
		if(!len[0]) {
			cp(0, 1); f[0][len[0] = 1] = 0;
			continue;
		}
		cp(0, 2);
		db mul = -dt[0] / dt[1];
		int st = i - fail[1]; len[0] = max(len[0], st + len[1]);
		f[0][st] -= mul;
		fo(j, 1, len[1]) f[0][st + j] += mul * f[1][j];
		if(len[2] - fail[2] < len[1] - fail[1]) cp(2, 1);
	}
}

int main() {
	freopen("a.in", "r", stdin);
	scanf("%d", &n);
	fo(i, 1, n) scanf("%lf", &a[i]);
	solve(a, n);
	pp("len = %d ", len[0]);
	fo(i, 1, len[0]) pp("%lf ", f[0][i]); hh;
	fo(i, 1, n) {
		db s = 0;
		fo(j, 1, min(i - 1, len[0])) s += a[i - j] * f[0][j];
		pp("%lf %lf\n", a[i], s);
	}
}

---

例题：

JZOJ6165

大概就是求稀疏图游走期望次数。
$n<=5000$

可以设$f[i]$表示在$i$回合后还没有结束的期望，$Ans=\sum f$

不会证明$f$是n阶的线性递归。

那么用$O(nm)$的时间求$f[1-n]$，然后套个BM，再解个方程就好了。

Code：

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, B = y; i <= B; i ++)
#define ff(i, x, y) for(int i = x, B = y; i <  B; i ++)
#define fd(i, x, y) for(int i = x, B = y; i >= B; i --)
#define ll long long
#define pp printf
#define db double
using namespace std;

const int mo = 998244353;

ll ksm(ll x, ll y) {
	ll s = 1;
	for(; y; y /= 2, x = x * x % mo)
		if(y & 1) s = s * x % mo;
	return s;
}

const int N = 1e4 + 5;

ll a[N];
ll f[3][N], len[3], fail[3], delta[3];

void cp(int x, int y) {
	fo(i, 1, len[x]) f[y][i] = f[x][i];
	len[y] = len[x]; fail[y] = fail[x]; delta[y] = delta[x];
}

void solve(ll *a, int n) {
	fo(i, 1, n) f[0][i] = 0; len[0] = 0;
	fo(i, 1, n) {
		ll tmp = a[i];
		fo(j, 1, len[0]) tmp -= f[0][j] * a[i - j] % mo;
		tmp = (tmp % mo + mo) % mo;
		delta[0] = tmp;
		if(!tmp) continue;
		fail[0] = i;
		if(!len[0]) {
			cp(0, 1); f[0][len[0] = 1] = 0;
			continue;
		}
		cp(0, 2);
		ll mul = delta[0] * ksm(delta[1], mo - 2) % mo;
		int st = i - fail[1]; len[0] = max(len[0], st + len[1]);
		f[0][st] = (f[0][st] + mul) % mo;
		fo(j, 1, len[1]) f[0][st + j] = (f[0][st + j] + (mo - mul) * f[1][j]) % mo;
		if(len[2] - fail[2] < len[1] - fail[1]) cp(2, 1);
	}
}

int T, n, m, x, y, b0[5005];
ll p[2][N][2], o, ni[N];
vector<int> b[N];

#define mem(a) memset(a, 0, sizeof a)

ll r[N];

void link(int x, int y) {
	b[x].push_back(y);
}

int main() {
	freopen("c.in", "r", stdin);
	freopen("c.out", "w", stdout);
	fo(i, 1, 5000) ni[i] = ksm(i, mo - 2);
	for(scanf("%d", &T); T; T --) {
		scanf("%d %d", &n, &m);
		fo(i, 1, n) {
			scanf("%d %d", &x, &y);
			b0[i] = 0;
			b[i].clear();
		}
		fo(i, 1, m) {
			scanf("%d %d", &x, &y);
			link(x, y); link(y, x);
			b0[x] ++; b0[y] ++;
		}
		mem(p);
		p[o][1][0] = 1;
		fo(i, 1, 2 * n) {
			mem(p[!o]);
			fo(j, 1, n - 1) {
				ll x = (p[o][j][0] + p[o][j][1]) % mo;
				ff(k, 0, b[j].size()) p[!o][b[j][k]][0] += x * ni[b0[j]] % mo;
			}
			fo(j, 1, n) {
				p[!o][j][0] = p[!o][j][0] % mo * ni[2] % mo;
				ff(k, 0, b[j].size()) p[!o][b[j][k]][1] += p[!o][j][0] * ni[b0[j]] % mo;
			}
			o = !o;
			a[i] = 0;
			fo(j, 1, n - 1) a[i] = (a[i] + p[o][j][0] + p[o][j][1]) % mo;
		}
		solve(a, 2 * n);
		ll s1 = 0, s2 = 0;
		fo(i, 1, len[0] + 1) {
			r[i] = a[i];
			fo(j, 1, i - 1) r[i] -= a[i - j] * f[0][j] % mo;
			s1 = (s1 + r[i] % mo + mo) % mo;
		}
		s2 = 1; fo(i, 1, len[0]) s2 = (s2 - f[0][i] + mo) % mo;
		pp("%lld\n", (s1 * ksm(s2, mo - 2) % mo + 1) % mo);
	}
}