Berlekamp-Massey 算法

最新推荐文章于 2023-10-23 15:38:58 发布

FSYo

最新推荐文章于 2023-10-23 15:38:58 发布

阅读量620

点赞数 1

分类专栏：多项式 BM

本文链接：https://blog.csdn.net/sslz_fsy/article/details/102943807

版权

多项式同时被 2 个专栏收录

37 篇文章 0 订阅

订阅专栏

3 篇文章 0 订阅

订阅专栏

膜拜 ldxoi

用途：

首先我们有一种很快的方法线性递推可以求一个递推式的第 $n$ 项
但是我们不知道这个递推式，只知道这个序列的前面几项
这个时候可以通过 $B M$ 算法解出最短递推式

算法流程：

对于一个 $n$ 项的序列 $A=\{a_1,a_2,...,a_n\}$ ，和一个 $m$ 项的序列 $R=\{r_1,r_2,...,r_n\}(m\le n)$
当 $\forall k>m$ ，有 $a_k=\sum_{i=1}^mr_ia_{k-i}$ ，那么 $R$ 为 $A$ 的一个递推式
$m$ 最小的为最短递推式
现在要通过增量算法求出这个递推式

假设我们处理到当前位置 $n$ ， ${a_1,a_2,...,a_{n-1}\}$ 的最短递推式是 ${r_1,r_2,...,r_m\}$
记第 $i$ 次更改的最短递推式是 $R_i$ ， $R_0=\{\empty\}$
记 $R_{cur}$ 为当前递推式
现在要让第 $n$ 项合法，于是
设 $\Delta n=a_n-\sum_{i=1}^mr_ia_{n-i}$ ，分类讨论

$\Delta n=0$ ，那么 $R_{cur}$ 为当前递推式
$\Delta n != 0$ ，我们需要构造出 $R_{cur+1}$ 来满足 $n$

对 $c u r$ 分类讨论

$c u r = 0$ ，那么构造一个 $R=\{r_0=0,r_1=0,...,r_n=0\}$
考虑构造一个增量递推式 $R_{\Delta}$ ，满足 $R_{cur+1}=R_{cur}+R_{\Delta}$

如何求 $R_{\Delta}$
定义 $fail_i$ 为第 $i$ 个递推式的第一个出错的位置
$fail_{cur}=n$
设 $R_{\Delta}$ 有 $m^{'}$ 项，要满足如下条件，在把 $n$ 归位时不影响别人

$\Delta n=\sum_{i=1}^{m'}R_{\Delta,i}*a_{n-i}$
$\forall m'<k<n,\sum_{i=1}^{m'}R_{\Delta,i}a_{k-i}=0$

考虑到上一层的 $R_{cur-1}$ 与 $R_{\Delta}$ 长得真像！

$\Delta fail_{cur-1}=a_{fail_{cur-1}}-\sum_{i=1}^{m_{cur-1}}R_{cur-1,i}*a_{fail_{cur-1}-i}$
$\forall m_{cur-1}<k<fail_{cur-1},a_k-\sum_{i=1}^{m_{cur-1}}R_{cur-1,i}a_{k-i}=0$

设 $t=\frac{\Delta n}{\Delta fail_{cur-1}},b=\{t,-t*R_{cur-1,1},-t*R_{cur-1,2},...,-t*R_{cur-1,m_{cur-1}}\}$

然后就可以按如下方式构造 $R_{\Delta}$

$\forall fail_{cur-1}<k<n$ ，给 $a_{n-k}$ 分配 0 的系数
$\forall fail_{cur-1}-m_{cur-1}\le k\le fail_{cur-1}$ ，给 $a_{n-k}$ 分配 $b_{fail_{cur-1}-k+1}$ 的系数

这样构造出来的 $R_{\Delta}=\{0,0,...,0,t,-t*r_{cur-1,1},-t*r_{cur-1,2},...,-t*r_{cur-1,m_{cur-1}}\}$
满足条件

理解：
对 $\Delta fail_{cur-1}=a_{fail_{cur-1}}-\sum_{i=1}^{m_{cur-1}}R_{cur-1,i}*a_{fail_{cur-1}-i}$ 乘上 $t$
$\Delta n=a_{fail_{cur-1}}*t-\sum_{i=1}^{m_{cur-1}}a_{fail_{cur-1}-i}*(t*R_{cur-1,i})$
然后用 ${a_{fail_{cur-1}}-m_{cur-1},...,a_{fail_{cur-1}}\}$ 去构造 $n$
然后用 $0$ 去满足 $fail_{cur-1}，n-1]$ 不变

$u p t : 19 / 12 / 02$
新的递推式 $\Delta R$ 的项数是 $n - f a i l [b e s t] + l e n [b e s t]$ ，因为要求最短递推式，所以还要动态维护一个 $b e s t$

模板

#include<bits/stdc++.h>
#define cs const
using namespace std;
int read(){
	int cnt = 0, f = 1; char ch = 0;
	while(!isdigit(ch)){ ch = getchar(); if(ch == '-') f = -1; }
	while(isdigit(ch)) cnt = cnt*10 + (ch-'0'), ch = getchar();
	return cnt * f;
}
typedef long long ll;
const int Mod = 998244353, G = 3;
cs int N = 5e5 + 5;
int add(int a, int b){ return a + b >= Mod ? a + b - Mod : a + b;}
int mul(int a, int b){ return 1ll * a * b % Mod;}
int ksm(int a, int b){ int ans = 1; for(;b;b>>=1){if(b&1) ans = mul(ans, a); a = mul(a, a);} return ans;}
void Add(int &a, int b){ a = add(a, b); }
int dec(int a, int b){return a - b < 0 ? a - b + Mod : a - b; }
ll inv[N];
#define poly vector<ll>
#define C 20
poly w[C + 1];
int up, bit, rev[N]; 
void Init(int len){ up = 1, bit = 0;	
	while(up < len) up <<= 1, bit++;
	for(int i = 0; i < up; i++) rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit-1));
}
void prework(){
	inv[0] = inv[1] = 1;
	for(int i = 2; i <= N-5; i++) inv[i] = mul(Mod-Mod/i, inv[Mod%i]);
	for(int i = 1; i <= C; i++) w[i].resize(1<<(i-1));
	ll wn = ksm(G, (Mod-1)/(1<<C)); w[C][0] = 1; 
	for(int i = 1; i < (1<<(C-1)); i++) w[C][i] = mul(w[C][i-1], wn);
	for(int i = C-1; i; i--)
		for(int j = 0; j < (1<<(i-1)); j++)
			w[i][j] = w[i+1][j<<1];
}
void NTT(poly &a, int flag){
	for(int i = 0; i < up; i++) if(i < rev[i]) swap(a[i], a[rev[i]]);
	for(int i = 1, l = 1; i < up; i <<= 1, l++)
		for(int j = 0; j < up; j += (i<<1))
			for(int k = 0; k < i; k++){
				ll x = a[k + j], y = mul(w[l][k], a[k + j + i]);
				a[k + j] = add(x, y); a[k + j + i] = add(x, Mod-y);
			}
	if(flag == -1){
		reverse(a.begin() + 1, a.begin() + up);
		for(int i = 0; i < up; i++) a[i] = mul(a[i], inv[up]);
	}
}
poly operator * (poly a, poly b){
	int len = a.size() + b.size() - 1;
	Init(len); a.resize(up); b.resize(up);
	NTT(a, 1); NTT(b, 1);
	for(int i = 0; i < up; i++) a[i] = mul(a[i], b[i]); 
	NTT(a, -1); a.resize(len); return a;
}
poly Inv(poly a, int len){
	poly b(1, ksm(a[0], Mod - 2)), c;
	for(int lim = 4; lim < (len << 2); lim <<= 1){
		Init(lim); c = a; c.resize(lim >> 1);
		c.resize(up); b.resize(up); NTT(c, 1); NTT(b, 1);
		for(int i = 0; i < up; i++) b[i] = mul(b[i], add(2, Mod - mul(b[i], c[i])));
		NTT(b, -1); b.resize(lim >> 1);
	} b.resize(len); return b;
}
poly operator - (poly a, poly b){
	poly c; int len = max(a.size(), b.size()); c.resize(len);
	a.resize(len); b.resize(len);
	for(int i = 0; i < len; i++) c[i] = add(a[i], Mod - b[i]);
	return c;
}
poly operator + (poly a, poly b){
	poly c; int len = max(a.size(), b.size()); c.resize(len);
	a.resize(len); b.resize(len);
	for(int i = 0; i < len; i++) c[i] = add(a[i], b[i]);
	return c;
}
poly operator / (poly a, poly b){
	int lim = 1, len = a.size() - b.size() + 1;
	reverse(a.begin(), a.end()); reverse(b.begin(), b.end());
	while(lim <= len) lim <<= 1;
	b = Inv(b, lim); b.resize(len);
	a = a * b; a.resize(len);
	reverse(a.begin(), a.end()); 
	return a;
}
poly operator % (poly a, poly b){ 
	int len = a.size() - b.size() + 1;
	if(len < 0) return a;
	poly c = a - (a / b) * b;
	c.resize(b.size() - 1); return c;
}
int calc(poly coe, int *a, int k, int n){
	poly p(k + 1), f(2), g(1); f[1] = 1; g[0] = 1; p[k] = 1; 
	for(int i = 1; i <= k; i++) p[k - i] = coe[i] ? Mod - coe[i] : 0;
	for(;n;n>>=1, f=f*f%p) if(n&1) g = g*f % p;
	int ans = 0;
	for(int i = 0; i < k; i++) ans = add(ans, mul(a[i], g[i]));
	return ans;
}
int n, m, a[N];
namespace BM{
	poly best, cur; int fail = 0, val = 0;
	void upt(int n){
		int t = a[n];
		for(int i = 1; i < cur.size(); i++) Add(t, Mod-mul(cur[i], a[n - i]));
		if(!t) return;
		if(cur.empty()){ cur.resize(n + 1); fail = n; val = t; return; }
		int coef = mul(t, ksm(val, Mod-2));
		poly nxt(n - fail);
		nxt.push_back(coef);
		for(int i = 1; i < best.size(); i++) nxt.push_back(mul(coef, dec(0, best[i])));
		nxt = nxt + cur;
		if((int)cur.size() - n <= (int)best.size() - fail){
			best = cur; fail = n; val = t;
		} cur = nxt;
	}
	int solve(){
		for(int i = 1, up = cur.size() - 1; i <= up; i++) cout << cur[i] << " ";
		puts("");
		return calc(cur, a + 1, cur.size() - 1, m);
	}
}
int main(){
	n = read(), m = read(); prework();
	for(int i = 1; i <= n; i++) a[i] = read(), BM::upt(i);
	cout << BM::solve(); return 0;
}

FSYo

关注

1
点赞
踩
2

收藏

觉得还不错? 一键收藏
打赏
0
评论
Berlekamp-Massey 算法

膜拜 ldxoi用途：首先我们有一种很快的方法线性递推可以求一个递推式的第 nnn 项但是我们不知道这个递推式，只知道这个序列的前面几项这个时候可以通过 BMBMBM 算法解出最短递推式算法流程：对于一个 nnn 项的序列 A={a1,a2,...,an}A=\{a_1,a_2,...,a_n\}A={a1,a2,...,an}，和一个 mmm 项的序列 R={r1,r2...
复制链接

扫一扫