【学习笔记】Berlekamp-Massey

orz zhenzhendong

之前贺过一边周指导博客但是弃疗了, 今天又贺了一次.

问题描述

给定一个长度为 \(n\) 的数列 \(\{a_i\}\), 求一个最短的齐次线性递推数列\(\{b_i\}\)（设长度为 \(m\)）,使得对于所有 \(m \leq k \leq n\), 有 \(a_k = \sum_{i = 1} ^ m a_{k - i} b_i\)

复杂度要求： \(O(n ^ 2)\)

一看就很适合配合常系数齐次线性递推食用.

算法流程

增量构造.

假设我们当前已经求出了 \(a_{0...i - 1}\) 的线性递推数列. 计算过程中, 我们曾经得出过 \(c\) 个递推式, 第 \(i\) 个递推式在 \(fail_i\) 的位置第一次失效了.

一开始 \(c = 0\), 我们有一个空的递推式.

现在我们加入数 \(a_i\).

设 \(R_c\) 的长度为 \(m\), \(delta_i = a_i - \sum_{k = 1} ^ {m} a_{i - k} R_c(k)\)

如果 \(delta_i = 0\), \(R_c\) 仍是一个合法的递推式.

否则我们要对 \(R_c\) 做出调整, 来得到一个新的符合条件的递推式.

若 \(c = 0\), 那么前 \(i - 1\) 个数都是 \(0\). 我们只需要构造一个包含 \(i\) 个 \(0\) 的递推式即可.

若 \(c \not= 0\), 只需要构造一个递推式 \(R'\), 当 \(|R'| + 1 \leq k < n\) 时, \(\sum_{i = 1} ^ {|R'|} a_{k - i} R'_i = 0\), \(\sum_{i = 1} ^ {|R'| } a_{n - i} R'_i = delta_n\), 那么 \(R_{c + 1} = R_c + R'\) 就符合条件.

我们随便找一个 $ 0 \leq id < c$, 它的前 \(fail_{id} - 1\) 个数都是 \(0\). 如果我们对它作一个位移, 即前面补上 \(i - fail_{id} - 1\) 个 \(0\), 后面跟个 \(1\), 然后接上 \(-R_{id}\), 我们就可以得到一个只有位置 \(i\) 是 \(delta_{fail_{id}}\),其余位置都是 \(0\) 的数组. 然后我们把它整个乘上 \(tmp = \frac{delta_{i}}{delta_{fail_{id}}}\), 就构造出了 \(R'\)

也就是 \(R'\) 为

\[\{0,0,...0,tmp,-tmp R_{id}(1),-tmp R_{id}(2),...\} \]

然后我们还要保证 \(R_c + R'\) 是最短的, 我们找到 \(i - fail_{id} + |R_{id}|\) 最小的即可. (不会严格证明)

模板

数据可以去周指导的博客上看.

#pragma GCC optimize("2,Ofast,inline")
#include<bits/stdc++.h>
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define LL long long
#define pii pair<int, int>
using namespace std;
const int mod = 1e9 + 7;

template <typename T> T read(T &x) {
	int f = 0;
	register char c = getchar();
	while (c > '9' || c < '0') f |= (c == '-'), c = getchar();
	for (x = 0; c >= '0' && c <= '9'; c = getchar())
		x = (x << 3) + (x << 1) + (c ^ 48);
	if (f) x = -x;
	return x;
}

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

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

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

inline int Qpow(int x, int p) {
	int ans = 1;
	for (; p; p >>= 1) {
		if (p & 1) ans = 1LL * ans * x % mod;
		x = 1LL * x * x % mod;
	}
	return ans;
}

inline int Inv(int x) {
	return Qpow(x, mod - 2);
}

namespace BM {

	const int Maxn = 5005;
	
	int n, c;
	int a[Maxn], del[Maxn], fail[Maxn];
	vector<int> R[Maxn];

	vector<int> solve() {
		c = 0;
		for (int i = 1; i <= n; ++i) {
			if (c == 0) {
				if (a[i]) {
					fail[0] = i;
					++c;
					del[i] = a[i];
					R[c].resize(i);
				}
				continue;
			}
			del[i] = a[i];
			for (int j = 0; j < R[c].size(); ++j) {
				del[i] = dec(del[i], 1LL * R[c][j] * a[i - j - 1] % mod);
			}
			if (del[i] == 0) continue;
			fail[c] = i;
			int id = c - 1, v = i - fail[id] + R[id].size();
			for (int j = c - 1; j >= 0; --j) {
				if (v > i - fail[j] + R[j].size()) {
					v = i - fail[j] + R[j].size();
					id = j;
				}
			}
			int p = i - fail[id];
			int tmp = 1LL * del[i] * Inv(del[fail[id]]) % mod;
			R[c + 1] = R[c];
			if (R[c + 1].size() < v) R[c + 1].resize(v);
			upd(R[c + 1][p - 1], tmp);
			for (int j = 0; j < R[id].size(); ++j) {
				upd(R[c + 1][p + j], -1LL * tmp * R[id][j] % mod + mod);
			}
			++c;
		}
		if (c == 0) return vector<int>(0);
		return R[c];
	}
}
using namespace BM;

int main() {
	read(n);
	for (int i = 1; i <= n; ++i) read(a[i]);
	vector<int> ans = BM::solve();
	cout << ans.size() << endl;
	for (int i = 0; i < ans.size(); ++i)
		cout << ans[i] << ' ';
	puts("");
}