可以用来求解数列的最短常系数线性递推式
即一个长度为
m
m
m的数列
r
r
r满足对于
i
∈
[
m
+
1
,
n
]
i\in[m+1,n]
i∈[m+1,n]
满足
A
i
=
∑
j
=
1
m
A
i
−
j
r
j
A_i=\sum_{j=1}^{m}A_{i-j}r_j
Ai=∑j=1mAi−jrj
B M BM BM算法复杂度是 O ( n 2 ) O(n^2) O(n2)的
首先有个东西叫做常系数齐次线性递推
就是对于一个长度为
k
k
k的递推式
可以在
O
(
k
l
o
g
n
l
o
g
k
)
O(klognlogk)
O(klognlogk)的时间求出第
n
n
n项
首先定义 r [ i ] r[i] r[i]为第 i i i次修改的递推式
当对于最新的一个
i
i
i不满足
∑
j
=
1
m
A
i
−
j
r
j
\sum_{j=1}^{m}A_{i-j}r_j
∑j=1mAi−jrj时
就会对递推式修改
定义
f
a
i
l
[
i
]
fail[i]
fail[i]表示第
i
i
i次修改
d
e
l
[
i
]
del[i]
del[i]表示第
i
i
i次修改前
a
[
f
a
i
l
[
i
]
]
a[fail[i]]
a[fail[i]]和用当时的递推式算出来的值之差
考虑当在
i
i
i出错时,设修改次数为
c
n
t
cnt
cnt
考虑对
r
[
c
n
t
]
r[cnt]
r[cnt]修改为
r
[
c
n
t
+
1
]
r[cnt+1]
r[cnt+1],使其对
i
i
i同样成立
如果当 c n t = 0 cnt=0 cnt=0,直接填入 i i i个 0 0 0,这样就相当于作为定义式,不用考虑当前这个了
否则设 m u l = d e l [ i ] d e l [ f a i l [ c n t − 1 ] ] mul=\frac{del[i]}{del[fail[cnt-1]]} mul=del[fail[cnt−1]]del[i]
我们只需要求出一个
r
′
=
r'=
r′={
r
1
,
r
2
,
r
3
…
…
r
m
′
r_1,r_2,r_3……r_{m'}
r1,r2,r3……rm′}$
使其满足
∀
i
∈
[
m
′
+
1
,
i
−
1
]
,
d
e
l
[
i
]
=
0
\forall i\in[m'+1,i-1] ,del[i]=0
∀i∈[m′+1,i−1],del[i]=0
并且 ∑ j = 1 m ′ A i − j r j ′ = d e l [ i ] \sum_{j=1}^{m'}A_{i-j}r'_{j}=del[i] ∑j=1m′Ai−jrj′=del[i]
然后只需要 r c n t + 1 = r c n t + r ′ r_{cnt+1}=r_{cnt}+r' rcnt+1=rcnt+r′即可
我们发现只需要令
r
′
=
r'=
r′={
0
,
0
,
…
…
0
,
m
u
l
,
−
m
u
l
∗
r
c
n
t
−
1
0,0,……0,mul,-mul*r_{cnt-1}
0,0,……0,mul,−mul∗rcnt−1}
将
i
−
f
a
i
l
c
n
t
−
1
−
1
i-fail_{cnt-1}-1
i−failcnt−1−1个0,一个
m
u
l
mul
mul,再把
r
c
n
t
−
1
r_{cnt-1}
rcnt−1乘上
−
m
u
l
-mul
−mul接在后面(注意是
−
-
−的)
可以发现这样是满足的(可以手推一下,不想具体写了)
UPD:
但是我们发现这样并不一定是最短的递推式,问题在于前面我们用的
c
n
t
−
1
cnt-1
cnt−1来构造
而
c
n
t
−
1
cnt-1
cnt−1的递推式不一定是最短的,这也是xyx神仙博客里的问题
这对一般求答案不会有影响,但是如果在写成生成函数递推算的时候不是最短就可能会出事
观察上面对递推式长度影响的,我们要找的是
[
1
,
c
n
t
−
1
]
[1,cnt-1]
[1,cnt−1]内
r
[
i
]
−
f
a
i
l
[
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];
}
}
这样最坏情况每次都会修改一次
r
r
r,修改
O
(
n
)
O(n)
O(n)次
复杂度就是
O
(
n
2
)
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();
}
例题: