浅谈线性递推
农民
线性递推相关
常系数齐次线性递推
对于一般的递归式，我们有
\(\sum\limits_{j=0}^{k}A_{i-j}R_j=0,i\ge k\)
定义\(S=AR\),则\(S\)的最高次为\(k-1\),小于\(R\)的最高次项\(k\)
即求\([x^n]A=\frac{S}{R}\)
若已知\(A\)的前\(k\)项,如果直接与\(R\)卷积即为\(S\),即
\(S=AR(mod \ x^k)\)
现在问题在于求\([x^n]A=\frac{S}{R}\)
Bostan_Mori
对于\([x^n]\dfrac{F(x)}{G(x)}=[x^n]\dfrac{F(x)G(-x)}{G(-x)G(x)}\)
令\(U(x)=F(x)G(-x),V(x)=G(-x)G(x)\)
注意到\(V(x)=V(-x)\)则\(V(x)\)只有偶数次项,考虑对\(U\)奇数分类
\([x^n]\dfrac{F(x)}{G(x)}=[x^n]\dfrac{U_1(x^2)+xU_2(x^2)}{V'(x^2)}\)
若\(n\)为奇数，即为\([x^{\frac{n}{2}}]\dfrac{U_2(x^2)}{V'(x^2)}\)
若\(n\)为偶数，即为\([x^{\frac{n}{2}}]\dfrac{U_1(x^2)}{V'(x^2)}\)
int Bostan_Mori(Poly F,Poly G,int N)
{
// printf("fic\n");
if(F.size()>=G.size())
{
pair<Poly,Poly>DJNB=Mod(F,G);
int Fsd=0;
if(N<DJNB.first.size())
{
Fsd=DJNB.first.U[N];
}
return ((long long)Bostan_Mori(DJNB.second,G,N)+(Fsd))%MOD;
}
else
{
while(N)
{
// printf("fic\n");
Poly GT=G;
for(int i=1;i<GT.size();i+=2)
{
GT.U[i]=(MOD-GT.U[i]);
}
Poly V=Mul_NTT(G,GT);
Poly U=Mul_NTT(F,GT);
F.U.clear();
G.U.clear();
for(int i=N&1;i<U.size();i+=2)
{
F.U.push_back(U.U[i]);
}
for(int i=0;i<V.size();i+=2)
{
G.U.push_back(V.U[i]);
}
N>>=1;
}
if(F.size()==0)
{
return 0;
}
return ((long long)F.U[0]*inv(G.U[0],MOD))%MOD;
}
}
int Homo_recursion(Poly A,Poly R,int n)
{
int k=A.size();
Poly S=Mul_NTT(A,R);
S.U.resize(k);
printf("%d\n",Bostan_Mori(S,R,n));
}
BM
已知有限递推序列\(⟨a_0,a_1,...a_n⟩\),求最短的递推式\(⟨r_0,r_1,...r_m⟩\)
即求最短的\(R\)使得\(AR=S(mod \ x^n)\)且\(|S|<|R|\)
定义\(R_n\)为满足\(AR=S(mod \ x^{n+1})\)的\(R\),\(S_n\)同理
如果已知\(Rn,Sn\)
若\(AR_n=S_n(mod \ x^{n+1})\),直接令\(R_{n+1}=R_n,S_{n+1}=S_n\)
否则\(AR_n=S_n+a_nx^n(mod \ x^{n+1})\),称为失配
考虑上一个失配点\(p\)
\(AR_p=S_p+a_px^p(mod \ x^{p+1})\)
\(x^{n-p}\frac{a_n}{a_p}AR_p=x^{n-p}\frac{a_n}{a_p}S_p+a_nx^n(mod \ x^{n+1})\)
相减即为\(A(R_n-x^{n-p}\frac{a_n}{a_p}R_p)=Sn-x^{n-p}\frac{a_n}{a_p}S_p(mod \ x^{n+1})\)
令\(R_{n+1}=R_n-x^{n-p}\frac{a_n}{a_p}R_p\),\(S_{n+1}=S_n-x^{n-p}\frac{a_n}{a_p}S_p\)
若\(|R_n|>|S_n|,|R_p|>|S_p|\),则\(|R_{n+1}|>|S_{n+1}|\)
关于验证最短
有一个结论不会证
\(|R_{n+1}|\ge \max(|R_n|,n+1-|R_n|)\)
至于最短,即每次都能取等
假设前\(n\)项均满足
如果\(n+1\)没有失配,则\(|R_{n+1}|=|R_n|\)
如果失配,由递推式得\(|R_{n+1}|=max(|R_n|,|R_p|+n-p)=|R_p|+n-p\),
而对于\(i\in(p,n],|R_i|=|R_{p+1}|,|R_{p+1}|=p+1-|R_p|\)
则\(|R_{n+1}|=|R_p|+n-p=n-p+p+1-|R_{p+1}|=n+1-|R_n|\)
至于初始的\(R_0,S_0\)令\(R_0=⟨1⟩,S_0=⟨0⟩(|S_0|=-\infin)\)
找到第一个不满足的位置\(p\),即\(AR_p=S_p+a_px^p(mod \ x^{p+1})\),而\(a_0,a_1,a_2...a_{p-1}=0\)
由上结论得\(|R_{p+1}|=p+1-|R_p|=p+1\),则\(R_{p+1}\)的最高次为\(p+1\)
考虑构造\(R_{p+1}=1+x^{p+1},S_{p+1}=a_{p}x^p\),则
\(R_{p+1}A=a_px^p(1+x^{p+1})=a_px^p=S_{p+1}(mod\ x^{p+1})\)
至于为什么\(R_{p+1}=1+x^{p+1}\)还是在\(mod\ x^{p+1}\)意义下，因为如果要最短,则必须取等(感觉有点奇怪,好像也构造不出更优的解)
Poly BM(Poly A){
int N=A.size();
Poly R,Rp;
int Pa,Pap;
R.clear();
R.push_back(1);
Poly Nex;
Nex.clear();
int P=-1;
for(int i=0;i<N;i++)
{
Nex.push_back(A.U[i]);
Pa=0;
for(int j=0;j<R.U.size();j++)
{
Pa=((long long)Pa+((long long)Nex.U[(i)-j]*R.U[j])%MOD)%MOD;
}
if(Pa==0)
{
continue;
}
else
{
if(P==-1)
{
R.U.resize(i+2,0);
R.U[0]=1;
Rp.U.clear();
Rp.U.push_back(1);
Pap=Pa;
P=i;
}
else
{
int IDSY=((long long)Pa*inv(Pap,MOD))%MOD;
if(Rp.size()+i-P>R.size())
{
Poly WORNMDESBTIR;
WORNMDESBTIR.clear();
WORNMDESBTIR.U.resize(Rp.size()+i-P,0);
for(int j=0;j<Rp.size();j++)
{
WORNMDESBTIR.U[j+i-P]=(MOD-((long long)Rp.U[j]*IDSY)%MOD);
}
for(int j=0;j<R.size();j++)
{
WORNMDESBTIR.U[j]=((long long)WORNMDESBTIR.U[j]+R.U[j])%MOD;
}
Rp=R;
Pap=Pa;
P=i;
R=WORNMDESBTIR;
}
else
{
for(int j=0;j<Rp.size();j++)
{
R.U[j+i-P]=((long long)R.U[j+i-P]-((long long)Rp.U[j]*IDSY)%MOD+MOD)%MOD;
}
}
}
}
}
return R;
}
