【Learning】常系数线性齐次递推

给定数列前k项(h_0...h_{k-1})，其后的项满足：(h_i=sum_{i=1}^kh_{i-j}a_i)，其中(a_1...a_k)是给定的系数，求(h_n)
　　
　　
　　
数据范围小的时候：

​　　做法一：暴力(O(nk))的DP

​　　做法二：矩阵快速幂.

　　　　　　记(H_i=egin{bmatrix}h_i&h_{i+1}&...&h_{i+k-1}end{bmatrix}). 则(h_n)是(H_{n-k+1})的最后一项。

　　　　　　(H_{n-k+1}=H_0M^{n-k+1})

　　　　　　其中(M)是转移矩阵，如当(k=4)时是这么填的：

[M=egin{bmatrix} 0&0&0&a_4\ 1&0&0&a_3\ 0&1&0&a_2\ 0&0&1&a_1 end{bmatrix} ]

　　　　　　时间复杂度(O(k^3lg n))
　　
　　
　　
数据范围大一些的时候：
　　
　　(kleq2000,nleq10^9). 这时候矩阵快速幂也做不了了
　　
　　还是拿(k=4)时举例，(M)的特征多项式(f(lambda))为：

[f(lambda)=det(lambda I-M)=egin{bmatrix} lambda&0&0&0\ 0&lambda&0&0\ 0&0&lambda&0\ 0&0&0&lambda end{bmatrix} -egin{bmatrix} 0&0&0&a_4\ 1&0&0&a_3\ 0&1&0&a_2\ 0&0&1&a_1 end{bmatrix}=egin{bmatrix} lambda&0&0&-a_4\ -1&lambda&0&-a_3\ 0&-1&lambda&-a_2\ 0&0&-1&lambda-a_1 end{bmatrix} ]

　　用行列式的性质，将(f(lambda))按最后一列拉普拉斯展开，得到如下，其中((-1)^{i+j}f(x)_{i,j})即行列式定义里的代数余子式：

[egin{aligned} f(lambda)&=sum_{i=1}^ka_{k-i+1}(-1)^{i+j}f(lambda)_{i,j} &取j=k（按最后一列展开）\ &=sum_{i=1}^ka_{k-i+1}(-1)^{i+k}f(lambda)_{i,k} end{aligned} ]

　　化简得到如下式子（也可以按(k=4)带进去看看规律）

[f(lambda)=lambda^k-sum_{i=1}^ka_ilambda^{k-i} ]

　　
　　现在明确一个定义，(f(x))这个函数的自变量(x)可以是实数，也可以是矩阵等等。这个函数仅仅是表示如何将自变量组合起来。表达的意思也会多样化，比如多项式、矩阵的多项式...下文会随时切换自变量的种类，但是函数的本质不变。
　　
　　(lambda)是(M)的特征值，是一个数。但是根据Cayley-Hamilton定理，如果把(lambda)替换成(M)代入得到(f(M)=M^k-sum_{i=1}^ka_iM^{k-i})，结果为一个零矩阵，即(M^k-sum_{i=1}^ka_iM^{k-i}=0)
　　

　　
　　我们想要求(M)的(n)次方（这里的(n)只是代表(M)的(n)次方，题目中(n)应该用(n-k+1)替代），然而(M^n)直接快速幂求不现实，复杂度为(O(k^3lg n)).

　　首先退一步考虑，要求一个数字的n次方(x^n)，如果我们把(x^n)对(f(x))取模会发生什么？

​　　根据多项式取模的定义，(x^n ; ext{mod}; f(x)=f(x)g(x)+r(x))，其中(g(x))和(r(x))是两个多项式.

​　　将(x)看成(M)，那么(f(M))为0.

　　故(M^n ; ext{mod}; f(M)=r(M))，且(M^n=M^n ; ext{mod}; f(M))，那么(M_n=r(M))这个多项式

​　　根据多项式取模的特性，(r(x))的次数严格小于模数(f(x))的次数(k). 那么(r(x))所包含的(M)的指数一定小于(k)，到达了可以计算的范围。

　　要求(M^n)，就只需要求(M^n ; ext{mod}; f(M))的多项式(r(M))。如果两个多项式(A(x))和(B(x))对模数取模分别得到(C(x))和(D(x))，那么多项式(A(x)B(x))对模数取模结果就是(C(x)D(x))。

​　　那么就可以用快速幂来求解(M^n ; ext{mod}; f(M))的结果了，也就是求出了(r(x))的各项系数（记为(c_i)）。实际计算中，表面上是在计算(M^n)，实际上计算的是(M^n ; ext{mod}; f(M))的结果。
　　
　　
　　
　　至此求出(r(x)=sumlimits_{i=0}^{k-1}c_ix^i). 将它看成矩阵的多项式代入(M)，得(r(M)=sumlimits_{i=0}^{k-1}c_iM^i)

​　　所以(M^n=sumlimits_{i=0}^{k-1}c_iM^i)

​　　把(n)替换成题目所需要的(n-k+1)，最终答案(h_n)为(H_0M^{n-k+1})的最后一项。

​

[H_0M^{n-k+1}=H_0sum_{i=0}^{k-1}c_iM^i=sum_{i=0}^{k-1}c_iH_0M_i=sum_{i=0}^{k-1}c_iH_i ]

　　那么要求的是(H_0M^{n-k+1})的最后一项。记(last(H_i)=h_{k+i}) ，那么

[h_n=last(H_0M^{n-k+1})=sum_{i=0}^{k-1}c_ilast(H_i)=sum_{i=0}^{k-1}c_ih_{i+k} ]

　　发现(i+kin[k,2k-1])，所以暴力算出(h_k...h_{2k-1})，代入求解得到(h_n)，至此全部求完。

　　分析复杂度：多项式乘法此处用暴力算会比FFT快，耗时最多的集快速幂求(r(x)) ，复杂度为(O(k^2lgn))。

#include <cstdio>
using namespace std;
const int K=4005,mod=1e9+7;
int n,k;
int a[K],h[K];
int b[K],c[K],t[K],mo[K];
inline void add(int &x,int y){
	x+=y;
	if(x>=mod) x-=mod;
}
void mul(int *x,int *y,int *z){
	for(int i=0;i<=2*k-2;i++) t[i]=0;
	for(int i=0;i<k;i++)
		for(int j=0;j<k;j++)
			add(t[i+j],1LL*x[i]*y[j]%mod);
	for(int i=2*k-2;i>=k;i--){
		for(int j=k-1;j>=0;j--)
			add(t[i-k+j],mod-1LL*t[i]*mo[j]%mod);
		t[i]=0;
	}
	for(int i=0;i<k;i++) z[i]=t[i];
}
void ksm(int y){
	for(;y;mul(b,b,b),y>>=1)
		if(y&1)
			mul(c,b,c);
}
int main(){
	freopen("input.in","r",stdin);
	scanf("%d%d",&n,&k); n++;
	for(int i=1;i<=k;i++){
		scanf("%d",&a[i]);
		if(a[i]<0) a[i]+=mod;
	}
	for(int i=1;i<=k;i++){
		scanf("%d",&h[i]);
		if(h[i]<0) h[i]+=mod;
	}
	mo[k]=1;
	for(int i=1;i<=k;i++) mo[k-i]=mod-a[i];
	if(n<=k){printf("%d
",h[n]);return 0;}
	b[1]=1; c[0]=1;
	ksm(n-k);
	for(int i=k+1;i<=2*k;i++)
		for(int j=1;j<=k;j++)
			add(h[i],1LL*a[j]*h[i-j]%mod);
	int ans=0;
	for(int i=0;i<k;i++) 
		add(ans,1LL*c[i]*h[i+k]%mod);
	printf("%d
",ans);
	return 0;
}

　　
　　

EXT

　　　　
　　如果(k)也比较大，那么要上多项式全家桶来优化多项式计算了！复杂度(O(klog klog n))
　　
　　来啊
　　　

#include <cstdio>
#include <vector>
#include <algorithm>
using namespace std;
typedef long long ll;
typedef vector<int> vi;
const int K=200005,mod=998244353,G=3;
int n,k,a[K],h[K];
inline void swap(int &x,int &y){int t=x;x=y;y=t;}
inline int max(int x,int y){return x>y?x:y;}
inline int min(int x,int y){return x<y?x:y;}
inline void add(int &x,int y){
	y=(y%mod+mod)%mod;
	(x+=y)%=mod;
}
inline int pow(int x,int y){
	int ret=1;
	for(;y;x=1LL*x*x%mod,y>>=1)
		if(y&1) ret=1LL*ret*x%mod;
	return ret;
}
namespace NTT{/*{{{*/
	int n,invn,bit,rev[K*4],A[K*4],B[K*4],W[K*4][2];
	void build(){
		int bas=pow(G,mod-2);
		for(int i=0;i<=18;i++){
			W[1<<i][0]=pow(G,(mod-1)/(1<<i));
			W[1<<i][1]=pow(bas,(mod-1)/(1<<i));
		}
	}
	void init(int na,int nb,vi &a,vi &b,int fn=0){
		if(!fn) fn=na+nb;
		for(n=1,bit=0;n<fn;n<<=1,bit++);
		invn=pow(n,mod-2);
		for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
		for(int i=0;i<n;i++) A[i]=B[i]=0;
		for(int i=0;i<na;i++) A[i]=a[i];
		for(int i=0;i<nb;i++) B[i]=b[i];
	}
	void ntt(int *a,int f){
		for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
		int w_n,w,u,v;
		for(int i=2;i<=n;i<<=1){
			w_n=W[i][f==-1];
			for(int j=0;j<n;j+=i){
				w=1;	
				for(int k=0;k<i/2;k++){
					u=a[j+k]; v=1LL*a[j+i/2+k]*w%mod;
					a[j+k]=(u+v)%mod;
					a[j+i/2+k]=(u+mod-v)%mod;
					w=1LL*w*w_n%mod;
				}
			}
		}
		if(f==1) return;
		for(int i=0;i<n;i++) a[i]=1LL*a[i]*invn%mod;
	}
	void calc(){
		ntt(A,1);
		ntt(B,1);
		for(int i=0;i<n;i++) A[i]=1LL*A[i]*B[i]%mod;
		ntt(A,-1);
	}
	void calchh(){
		ntt(A,1);
		ntt(B,1);
		for(int i=0;i<n;i++) A[i]=(2LL*B[i]%mod+mod-1LL*A[i]*B[i]%mod*B[i]%mod)%mod;
		ntt(A,-1);
	}
}/*}}}*/
vi mop,b,c,T;
vi operator - (vi A,vi B){
	int n=A.size(),m=B.size(),fn=max(n,m);
	A.resize(fn);
	for(int i=0;i<m;i++) add(A[i],-B[i]);
	return A;
}
vi operator * (int a,vi A){
	int n=A.size();
	a=(a+mod)%mod;
	for(int i=0;i<n;i++) A[i]=1LL*a*A[i]%mod;
	return A;
}
vi operator * (vi &A,vi B){
	int n=A.size(),m=B.size();
	NTT::init(n,m,A,B);
	NTT::calc();
	A.resize(n+m-1);
	for(int i=0;i<n+m-1;i++) A[i]=NTT::A[i];
	return A;
}
vi inverse(vi A){
	int n=A.size();
	if(n==1){
		A[0]=pow(A[0],mod-2);
		return A;
	}
	vi B=A;
	B.resize((n+1)/2);
	B=inverse(B);

	int m=B.size();
	NTT::init(n,m,A,B,n+m-1+m-1);
	NTT::calchh();
	B.resize(NTT::n);
	for(int i=0;i<NTT::n;i++) B[i]=NTT::A[i];

	//B=(2*B)-((A*B)*B);
	B.resize(n);
	return B;
}
vi operator / (vi A,vi B){
	int n=A.size()-1,m=B.size()-1;
	vi C;
	if(n<m){
		C.resize(1); C[0]=0;
		return C;
	}
	reverse(A.begin(),A.end());
	reverse(B.begin(),B.end());
	B.resize(n-m+1);
	C=A*inverse(B);
	C.resize(n-m+1);	
	reverse(C.begin(),C.end());
	return C;
}
void module(vi &A,vi B){
	int n=A.size()-1,m=B.size()-1;
	if(n<m) return;
	vi D=A/B;
	A=A-(B*D);
	A.resize(m);
}
void ksm(int y){
	for(;y;y>>=1){
		if(y&1){
			c=c*b;
			module(c,mop);
		}
		b=b*b;
		module(b,mop);
	}
}
int main(){
	freopen("input.in","r",stdin);
	NTT::build();
	scanf("%d%d",&n,&k); n++;
	for(int i=1;i<=k;i++) scanf("%d",&h[i]),h[i]%=mod;
	for(int i=1;i<=k;i++) scanf("%d",&a[i]),a[i]%=mod;
	if(n<=k){printf("%d
",h[n]);return 0;}
	mop.resize(k+1);
	mop[k]=1;
	for(int i=1;i<=k;i++) mop[k-i]=(mod-a[i])%mod;
	b.resize(2); b[1]=1; 
	c.resize(1); c[0]=1;
	ksm(n-1);
	int ans=0;
	c.resize(k);
	for(int i=0;i<k;i++)
		add(ans,1LL*c[i]*h[i+1]%mod);
	printf("%d
",ans);
	return 0;
}

参考资料

http://blog.csdn.net/qq_33229466/article/details/78933309 "ORZ"