zoukankan      html  css  js  c++  java
  • 【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"

  • 相关阅读:
    Understanding about Baire Category Theorem
    Isometric embedding of metric space
    Convergence theorems for measurable functions
    Mindmap for "Principles of boundary element methods"
    Various formulations of Maxwell equations
    Existence and uniqueness theorems for variational problems
    Kernels and image sets for an operator and its dual
    [loj6498]农民
    [luogu3781]切树游戏
    [atAGC051B]Three Coins
  • 原文地址:https://www.cnblogs.com/RogerDTZ/p/8503773.html
Copyright © 2011-2022 走看看