[Template][Luogu P5487]-Berlekamp-Massey algorithm

Portal

Linear Recursive Formula with Minimum Constant Coefficients Used to Solve Sequences of Numbers
That is to say, a sequence rrr with length of m m m satisfies for I [m+1,n] i in [m+1,n] I [m+1,n]
Satisfy A i=j=1mAi_j rj A_i=sum_{j=1}^{m} A_{i-j} r_jAi=j=1m Ai_j rj

The complexity of BMBMBM algorithm is O(n2)O(n^2)O(n2)

Firstly, r[i]r[i]r[i] is defined as the recursive formula of the second revision of iii.

When j=1m A i_j rj_sum_{j=1}^{m} A_{i-j} r_j_j=1m Ai_j RJ is not satisfied with the latest iii
Recursive modifications will be made

Define fail[i]fail[i]fail[i] to denote the second modification of iii
del[i]del[i]del[i] denotes the difference between a[fail[i]]a[fail[i]]a[fail[i]] a [fail [i]] before the second revision of iii and the value calculated by the recursive formula at that time.

Consider setting the number of modifications to cntcnt when an error occurs in iii
Consider modifying r[cnt]r[cnt]r[cnt] R [cnt] to r[cnt+1]r[cnt+1]r[cnt+1], so that it is also valid for iii

If cnt=0cnt=0cnt=0 CNT = 0, directly fill in iii 000, which is equivalent to the definition, without considering the current one.

Otherwise, let mul=del[i]del[fail[cnt_1]] mul =frac {del [i]} {del [fail [cnt-1]]} mul = del [fail [cnt_1]]del[i]

We just need to find a r'= R & x27; = r'={r1,r2,r3... rm'r_1,r_2,r_3... R_{m & x 27;} r1,r2,r3... rm '$$
Let it satisfy [m'+1,i_1], del[i]=0_for all i\ in [m '+1,i-1], del [i]=0[m'+1,i_1], del [i]=0, del [i]=0

And j=1m'A i_j rj'= del[i] sum {j=1} ^ {m & #x27;} A {i-j} R & #x27; {j}= del[i] J = 1m Ai j'= del[i]

Then you just need r CNT + 1 = rcnt + r'r {cnt + 1} = R {cnt} + R & # x27; rcnt+1=rcnt+r'

We found that we only need to make r'= R & x27; = r'={0,0,... 0,mul, mul rcnt 10, 0,... 0,mul,-mul*r_{cnt-1}0,0,... 0,mul,mul_rcnt_1}
Multiply i_fail cnt_1_1i-fail_{cnt-1}-1i_failcnt_1_1 0, one mulmulmulmulmulmulmulmul, and then multiply rcnt_1r_{cnt-1}rcnt_1 by_mul

You can find that this is satisfying (you can push it by hand, don't want to write it in detail)

The worst case scenario would be to modify rrr and O(n)O(n)O(n) O (n) times each time.
Complexity is O(n2)O(n^2)O(n2)

Code:

#include<bits/stdc++.h>
using namespace std;
#define gc getchar
inline int read(){
	char ch=gc();
	int res=0,f=1;
	while(!isdigit(ch))f^=ch=='-',ch=gc();
	while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
	return f?res:-res;
}
#define re register
#define pb push_back
#define cs const
#define pii pair<int,int>
#define fi first
#define se second
#define ll long long
#define poly vector<int>
#define bg begin
cs int mod=998244353,G=3;
inline int add(int a,int b){return (a+=b)>=mod?a-mod:a;}
inline void Add(int &a,int b){(a+=b)>=mod?(a-=mod):0;}
inline int dec(int a,int b){return (a-=b)<0?a+mod:a;}
inline void Dec(int &a,int b){(a-=b)<0?(a+=mod):0;}
inline int mul(int a,int b){return 1ll*a*b>=mod?1ll*a*b%mod:a*b;}
inline void Mul(int &a,int b){a=mul(a,b);}
inline int ksm(int a,int b,int res=1){
	for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));return res;
}
inline void chemx(ll &a,ll b){a<b?a=b:0;}
inline void chemn(int &a,int b){a>b?a=b:0;}
cs int N=(1<<16)|7,C=16;
poly w[C+1];
inline void init_w(){
	for(int i=1;i<=C;i++)w[i].resize(1<<(i-1));
	int 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];
}
int rev[N<<2];
inline void init_rev(int lim){
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)*(lim>>1));
}
inline void ntt(poly &f,int lim,int kd){
	for(int i=0;i<lim;i++)if(i>rev[i])swap(f[i],f[rev[i]]);
	for(int mid=1,l=1;mid<lim;mid<<=1,l++)
	for(int i=0;i<lim;i+=(mid<<1))
		for(int j=0,a0,a1;j<mid;j++)
			a0=f[i+j],a1=mul(f[i+j+mid],w[l][j]),
			f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
	if(kd==-1&&(reverse(f.bg()+1,f.bg()+lim),1))
		for(int inv=ksm(lim,mod-2),i=0;i<lim;i++)Mul(f[i],inv);
}
inline poly operator +(poly a,poly b){
	poly c;int lim=max(a.size(),b.size());c.resize(lim);
	a.resize(lim),b.resize(lim);
	for(int i=0;i<lim;i++)c[i]=add(a[i],b[i]);return c;
}
inline poly operator -(poly a,poly b){
	poly c;int lim=max(a.size(),b.size());c.resize(lim);
	a.resize(lim),b.resize(lim);
	for(int i=0;i<lim;i++)c[i]=dec(a[i],b[i]);return c;
}
inline poly operator *(poly a,int b){
	for(int i=0;i<a.size();i++)Mul(a[i],b);return a;
}
inline poly operator /(poly a,int b){
	for(int i=0,inv=ksm(b,mod-2);i<a.size();i++)Mul(a[i],inv);
	return a;
}
inline poly operator *(poly a,poly b){
	int deg=a.size()+b.size()-1,lim=1;
	if(deg<=128){
		poly c(deg,0);
		for(int i=0;i<a.size();i++)
		for(int j=0;j<b.size();j++)
			Add(c[i+j],mul(a[i],b[j]));
		return c;
	}
	while(lim<deg)lim<<=1;init_rev(lim);
	a.resize(lim),ntt(a,lim,1);
	b.resize(lim),ntt(b,lim,1);
	for(int i=0;i<lim;i++)Mul(a[i],b[i]);
	ntt(a,lim,-1),a.resize(deg);
	return a;
}
inline poly Inv(poly a,int deg){
	poly c,b(1,ksm(a[0],mod-2));
	for(int lim=4;lim<(deg<<2);lim<<=1){
		init_rev(lim);
		c=a,c.resize(lim>>1);
		c.resize(lim),ntt(c,lim,1);
		b.resize(lim),ntt(b,lim,1);
		for(int i=0;i<lim;i++)Mul(b[i],dec(2,mul(b[i],c[i])));
		ntt(b,lim,-1),b.resize(lim>>1);
	}b.resize(deg);return b;
}
inline poly operator /(poly a,poly b){
	int lim=1,deg=(int)a.size()-(int)b.size()+1;
	reverse(a.bg(),a.end());
	reverse(b.bg(),b.end());
	while(lim<=deg)lim<<=1;
	b=Inv(b,lim),b.resize(deg);
	a=a*b,a.resize(deg);
	reverse(a.bg(),a.end());
	return a;
}
inline poly operator %(poly a,poly b){
	int deg=(int)a.size()-(int)b.size()+1;
	if(deg<0)return a;
	poly c=a-(a/b)*b;
	c.resize(b.size()-1);
	return c;
}
inline poly ksm(poly a,int b,poly res,cs poly &mod){
	for(;b;b>>=1,a=a*a%mod)if(b&1)res=res*a%mod;
	return res;
}
namespace Cayley_Hamilton{
	int n;
	inline int solve(poly coef,int *a,int k){
		n=coef.size(),init_w();
		poly f(n+1),g(2),res(1);g[1]=res[0]=1;
		for(int i=1; i<n;i++)f[n-i]=dec(0,coef[i]);
		f[n]=1;
		res=ksm(g,k,res,f);
		int anc=0;
		for(int i=0;i<n;i++)Add(anc,mul(res[i],a[i+1]));
		return anc;
	}
}
namespace Berlekamp_Massey{
	poly r[N];
	int n,cnt,m,a[N],fail[N],del[N];
	inline void update(int i){
		++cnt;
		int MUL=mul(dec(a[i],del[i]),ksm(dec(a[fail[cnt-2]],del[fail[cnt-2]]),mod-2));
		r[cnt].resize(i-fail[cnt-2],0);
		r[cnt].pb(MUL);
		for(int j=1;j<r[cnt-2].size();j++){
			r[cnt].pb(mul(mod-r[cnt-2][j],MUL));
		}
		r[cnt]=r[cnt]+r[cnt-1];
	}
	inline void BM(){
		for(int i=1;i<=n;i++){
			for(int j=1;j<r[cnt].size();j++){
				Add(del[i],mul(a[i-j],r[cnt][j]));
			}
			if(a[i]!=del[i]){
				fail[cnt]=i;
				if(!cnt)r[++cnt].resize(i+1);
				else update(i);
			}
		}
	}
	inline void solve(){
		n=read(),m=read();
		for(int i=1;i<=n;i++)
		a[i]=read();
		BM();
		for(int i=1;i<r[cnt].size();i++)
		cout<<r[cnt][i]<<" ";puts("");
		cout<<Cayley_Hamilton::solve(r[cnt],a,m)<<'\n';
	}
}
int main(){
	Berlekamp_Massey::solve();
}

Posted on Tue, 08 Oct 2019 03:36:11 -0700 by Liodel