【模板】【洛谷P5487】—Berlekamp-Massey算法

传送门

可以用来求解数列的最短常系数线性递推式
即一个长度为mmm的数列rrr满足对于i∈[m+1,n]i\in[m+1,n]i∈[m+1,n]
满足Ai=∑j=1mAi−jrjA_i=\sum_{j=1}^{m}A_{i-j}r_jAi=∑j=1mAi−jrj

BMBMBM算法复杂度是O(n2)O(n^2)O(n2)的

首先有个东西叫做常系数齐次线性递推
就是对于一个长度为kkk的递推式
可以在O(klognlogk)O(klognlogk)O(klognlogk)的时间求出第nnn项

传送门

首先定义r[i]r[i]r[i]为第iii次修改的递推式

当对于最新的一个iii不满足∑j=1mAi−jrj\sum_{j=1}^{m}A_{i-j}r_j∑j=1mAi−jrj时
就会对递推式修改

定义fail[i]fail[i]fail[i]表示第iii次修改
del[i]del[i]del[i]表示第iii次修改前a[fail[i]]a[fail[i]]a[fail[i]]和用当时的递推式算出来的值之差

考虑当在iii出错时，设修改次数为cntcntcnt
考虑对r[cnt]r[cnt]r[cnt]修改为r[cnt+1]r[cnt+1]r[cnt+1]，使其对iii同样成立

如果当cnt=0cnt=0cnt=0，直接填入iii个000，这样就相当于作为定义式，不用考虑当前这个了

否则设mul=del[i]del[fail[cnt−1]]mul=\frac{del[i]}{del[fail[cnt-1]]}mul=del[fail[cnt−1]]del[i]

我们只需要求出一个r′=r'=r′={r1,r2,r3……rm′r_1,r_2,r_3……r_{m'}r1,r2,r3……rm′}$
使其满足∀i∈[m′+1,i−1]，del[i]=0\forall i\in[m'+1,i-1] ，del[i]=0∀i∈[m′+1,i−1]，del[i]=0

并且∑j=1m′Ai−jrj′=del[i]\sum_{j=1}^{m'}A_{i-j}r'_{j}=del[i]∑j=1m′Ai−jrj′=del[i]

然后只需要rcnt+1=rcnt+r′r_{cnt+1}=r_{cnt}+r'rcnt+1=rcnt+r′即可

我们发现只需要令r′=r'=r′={0,0,……0,mul,−mul∗rcnt−10,0,……0,mul,-mul*r_{cnt-1}0,0,……0,mul,−mul∗rcnt−1}
将i−failcnt−1−1i-fail_{cnt-1}-1i−failcnt−1−1个0，一个mulmulmul,再把rcnt−1r_{cnt-1}rcnt−1乘上−mul-mul−mul接在后面（注意是−-−的）

可以发现这样是满足的（可以手推一下，不想具体写了）

UPD：

但是我们发现这样并不一定是最短的递推式，问题在于前面我们用的cnt−1cnt-1cnt−1来构造
而cnt−1cnt-1cnt−1的递推式不一定是最短的，这也是xyx神仙博客里的问题
这对一般求答案不会有影响，但是如果在写成生成函数递推算的时候不是最短就可能会出事

观察上面对递推式长度影响的，我们要找的是[1,cnt−1][1,cnt-1][1,cnt−1]内r[i]−fail[i]r[i]-fail[i]r[i]−fail[i]最小的来更新
维护一个最短的即可
代码（后面的就没改了）：

namespace BM{int fail[N],del[N];
vector<int> r[N],t;
inline poly work(poly a,int n){int cnt=0,mx=0;for(int i=1;i<=n;i++){int w=a[i];for(int j=0;j<r[cnt].size();j++)Dec(w,mul(r[cnt][j],a[i-j-1]));if(w==0)continue;del[cnt]=w,fail[cnt]=i;if(!cnt){r[++cnt].resize(i);continue;}int ml=mul(del[cnt],Inv(del[mx]));r[++cnt].resize(i-fail[mx]-1),r[cnt].pb(ml);for(int j=0;j<r[mx].size();j++)r[cnt].pb(dec(0,mul(r[mx][j],ml)));r[cnt]=r[cnt]+r[cnt-1];if(r[cnt-1].size()-i<r[mx].size()-fail[mx])mx=cnt-1;}return r[cnt];
}}

这样最坏情况每次都会修改一次rrr，修改O(n)O(n)O(n)次
复杂度就是O(n2)O(n^2)O(n2)的

代码：

// luogu-judger-enable-o2
// luogu-judger-enable-o2
#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,ll k){n=coef.size(),init_w();poly f(n),g(2),res(1,0);g[1]=res[0]=1;for(int i=1; i<n;i++)f[n-i-1]=dec(0,coef[i]);f[n-1]=1;res=ksm(g,k,res,f);int anc=0;for(int i=0;i<res.size();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();
}

例题：

【HDU6172】

【Codechef walk on tree】

【模板】【洛谷P5487】—Berlekamp-Massey算法相关推荐

洛谷专题训练 ——【算法1-1】模拟与高精度
洛谷题单[算法1-1]模拟与高精度 ACM-ICPC在线模板题单链接: [算法1-1]模拟与高精度下面的这一坨都是洛谷题单上的东东题单简介恭喜大家完成了第一部分语言入门,相信大家已经可以使用 ...
Berlekamp–Massey算法简要介绍
这是一篇翻译向的文章,笔者整理了一些有关Berlekamp–Massey算法的笔记,还增加了一些自己的理解. 下面列出了笔者写此文时所参考的一些资料: wikipedia fjzzq2002 别人的博 ...
高精度模板洛谷Luogu P1932 A+B A-B A*B A/B Problem
P1932 A+B & A-B & A*B & A/B Problem 题目背景这个题目很新颖吧!!! 题目描述求A.B的和差积商余! 输入输出格式输入格式: 两个数两行 ...
Berlekamp–Massey 算法
学习博客学习博客 BM算法用于求解常系数线性递推式. 它可以在 O(n2)O(n^2)O(n2) 的时间复杂度内解决问题. 先打个比赛,一会再更学习算法看上边的两个学习链接就行,板子可以参考我下边 ...
字典树模板+洛谷P2580 于是他错误的点名开始了
题目: 题目背景 XS中学化学竞赛组教练是一个酷爱炉石的人. 他会一边搓炉石一边点名以至于有一天他连续点到了某个同学两次,然后正好被路过的校长发现了然后就是一顿欧拉欧拉欧拉(详情请见已结束比赛 CON ...
扩展 KMP(模板) 洛谷P5410
题目:https://www.luogu.org/problemnew/show/P5410 博客 #include<bits/stdc++.h>using namespace std;# ...
[模板]洛谷T3379 最近公共祖先(LCA) 倍增+邻接表
一年前听说的这东西...现在终于会了... 1 #include<cstdio> 2 #include<iostream> 3 #include<cstring> ...
Berlekamp Massey算法(HDU6172)
用途:常系数k阶递推式求值 HDU6172 #include <bits/stdc++.h>using namespace std; #define rep(i, a, n) for (l ...
洛谷or牛客数据结构+算法
栈思想:先进后出 tips:栈里能放下标就放下标 (牛客)小c的计事本(直接用stack可以简化代码,且不会被自己绕晕,当时没意识到) (牛客)吐泡泡(没意识到用栈),(牛客)好串 1.后缀表达式(栈 ...

【模板】【洛谷P5487】—Berlekamp-Massey算法

UPD：

【模板】【洛谷P5487】—Berlekamp-Massey算法相关推荐

最新文章

热门文章