• 【模板】常系数齐次线性递推


    XXVII.【模板】常系数齐次线性递推

    题意:已知\(f_0,\dots,f_{m-1}\),且对于\(k\geq m\),有

    \[f_k=\sum\limits_{i=1}^ma_if_{k-i} \]

    其中\(a_1,\dots,a_m\)是给定的系数。

    \(f_n\)

    我们一个naive的思路就是矩阵快速幂。

    考虑设\(m\)阶向量

    \[\mathbf{V}=\left[\begin{matrix}f_0\\\vdots\\f_{m-1}\end{matrix}\right] \]

    再设矩阵\(\mathrm{A}\)为构造出来的转移矩阵。

    则有

    \[f_n=(\mathbf{V}\mathrm{A}^n)_0 \]

    这个复杂度是\(O(m^3\log n)\),无法承受。

    我们假设我们已经通过魔法找到了一组\(c_0,\dots,c_{m-1}\)使得

    \[\mathrm{A}^n=\sum\limits_{i=0}^{m-1}c_i\mathrm{A}^i \]

    则两边同时乘上\(\mathbf{V}\),得到

    \[\mathbf{V}\mathrm{A}^n=\sum\limits_{i=0}^{m-1}c_i\mathbf{V}\mathrm{A}^i \]

    因为我们只考虑求出\(f_n\)\((\mathbf{V}\mathrm{A}^n)_0\),所以我们在上式中也只需考虑\(0\)一位即可。

    于是

    \[(\mathbf{V}\mathrm{A}^n)_0=\sum\limits_{i=0}^{m-1}c_i(\mathbf{V}\mathrm{A}^i)_0 \]

    \[f_n=\sum\limits_{i=0}^{m-1}c_if_i \]

    则我们现在的目标就变成求出这样一组\(c\)出来。

    我们设一个矩阵函数\(C(\mathrm{A})=\sum\limits_{i=0}^{m-1}c_i\mathrm{A}^i=\mathrm{A}^n\)

    显然两边次数不等,不好求出。我们考虑另外找三个多项式\(Q,G,R\),满足

    \[\mathrm{A}^n=(QG)(\mathrm{A})+R(\mathrm{A}) \]

    即多项式除法的形式。则应有\(R\)的次数小于\(G\)的次数。

    如果\(G\)的次数是\(m\)的话,就可以变成如下形式

    \[\mathrm{A}^n=(QG)(\mathrm{A})+\sum\limits_{i=0}^{m-1}r_i\mathrm{A}^i \]

    此时,如果我们可以幸运地发现有\(G(\mathrm{A})=\mathrm{0}\)的话,最终就有

    \[\mathrm{A}^n=R(\mathrm{A})=C(\mathrm{A}) \]

    但真有这种好事吗?

    我们考虑,有

    \[\mathrm{A}=\begin{bmatrix}0,1,0,0,\dots,0\\0,0,1,0,\dots,0\\0,0,0,1,\dots,0\\\vdots\\0,0,0,0,\dots,1\\a_m,a_{m-1}\dots,a_1\end{bmatrix} \]

    我们发现,假如

    \[g_i=-a_{m-i},g_m=1 \]

    的话,恰有\(G(\mathrm{A})=0\)

    故我们只需要得到\(G\)后,计算

    \[\mathrm{A}^n\bmod{G} \]

    即可得到\(C\)

    多项式取模可以用多项式除法魔改得到。然后套上一个快速幂即可。

    时间复杂度\(O(m\log m\log n)\)

    代码:

    #include<bits/stdc++.h>
    using namespace std;
    #include<bits/stdc++.h>
    using namespace std;
    const int mod=998244353;
    const int N=1<<20;
    int ksm(int x,int y){
    	int rt=1;
    	for(;y;x=(1ll*x*x)%mod,y>>=1)if(y&1)rt=(1ll*rt*x)%mod;
    	return rt;
    }
    int n,m,f[N],g[N],invrevg[N],all,res;
    namespace Poly{
    	const int G=3;
    	int rev[N],lim,invlim;
    	void NTT(int *a,int tp){
    		for(int i=0;i<lim;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    		for(int md=1;md<lim;md<<=1){
    			int rt=ksm(G,(mod-1)/(md<<1));
    			if(tp==-1)rt=ksm(rt,mod-2);
    			for(int stp=md<<1,pos=0;pos<lim;pos+=stp){
    				int w=1;
    				for(int i=0;i<md;i++,w=(1ll*w*rt)%mod){
    					int x=a[pos+i],y=(1ll*w*a[pos+md+i])%mod;
    					a[pos+i]=(x+y)%mod;
    					a[pos+md+i]=(x-y+mod)%mod;
    				}
    			}
    		}
    		if(tp==-1)for(int i=0;i<lim;i++)a[i]=(1ll*a[i]*invlim)%mod;
    	}
    	int A[N],B[N],C[N],D[N],E[N];
    	void mul(int *a,int *b,int *c,int LG=-1){//using: Array A and B
    		if(LG!=-1){
    			lim=(1<<LG),invlim=ksm(lim,mod-2);
    			for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(LG-1));			
    		}
    		for(int i=0;i<lim;i++)A[i]=B[i]=0;
    		for(int i=0;i<(lim>>1);i++)A[i]=a[i],B[i]=b[i];
    		NTT(A,1),NTT(B,1);
    		for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
    		NTT(A,-1);
    		for(int i=0;i<lim;i++)c[i]=A[i];
    	}
    	void inv(int *a,int *b,int LG){//using: Array C
    		b[0]=ksm(a[0],mod-2);
    		for(int k=1;k<=LG+1;k++){
    			mul(b,a,C,k);
    			for(int i=0;i<(1<<k);i++)C[i]=(mod-C[i])%mod;
    			(C[0]+=2)%=mod;
    			mul(C,b,b,k);
    		}
    	}
    	void rever(int *a,int *b,int lim){
    		for(int i=0;i<lim;i++)b[i]=a[i];
    		reverse(b,b+lim);
    	}
    	void modulo(int *a){//using:Array D
    		rever(a,D,2*m);
    		for(int i=m;i<lim;i++)D[i]=0;
    		mul(D,invrevg,E);
    		reverse(E,E+m);
    		for(int i=m;i<lim;i++)E[i]=0;
    		mul(E,g,D);
    		for(int i=0;i<m;i++)a[i]=(a[i]-D[i]+mod)%mod;
    		for(int i=m;i<lim;i++)a[i]=0;
    	}
    }
    using namespace Poly;
    int c[N],d[N];
    void KSM(){
    	c[1]=1,d[0]=1;
    	for(;n;n>>=1,mul(c,c,c),modulo(c))if(n&1)mul(d,c,d),modulo(d);
    }
    int main(){
    	scanf("%d%d",&n,&m);
    	for(int i=1;i<=m;i++)scanf("%d",&g[m-i]),g[m-i]=-g[m-i];
    	g[m]=1;
    	for(int i=0;i<=m;i++)(g[i]+=mod)%=mod;
    	for(int i=0;i<m;i++)scanf("%d",&f[i]),(f[i]+=mod)%=mod;
    	while((1<<all)<(m<<1))all++;
    	reverse(g,g+m+1);
    	inv(g,invrevg,all);
    	reverse(g,g+m+1);
    	for(int i=m+1;i<(1<<all);i++)invrevg[i]=0;
    	lim=(1<<all),invlim=ksm(lim,mod-2);
    	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(all-1)); 
    	KSM();
    	for(int i=0;i<m;i++)(res+=1ll*d[i]*f[i]%mod)%=mod;
    	printf("%d\n",res);
    	return 0;
    } 
    

  • 相关阅读:
    iOS iPad开发之Modal
    iOS SVN终端指令
    iOS iPad开发之UIPopoverController的使用
    iOS Xcode6和Xcode5的区别?
    算法 查找算法--二分查找
    算法 排序算法--快速排序
    算法 排序算法--选择排序
    XCode签名证书死活不能选
    IOS7 适配时导航栏变黑
    Xcode Provisioning 路径
  • 原文地址:https://www.cnblogs.com/Troverld/p/14608069.html
Copyright © 2020-2023  润新知