zoukankan      html  css  js  c++  java
  • 常系数齐次线性递推

    常系数齐次线性递推

    参考资料

    https://blog.csdn.net/hzj1054689699/article/details/85683342

    https://www.luogu.com.cn/blog/Zhang-RQ/chang-ji-shuo-ji-ci-xian-xing-di-tui-chu-tan

    问题描述

    见洛谷模板题。

    有个数列({a}),给出前(k)项,即(a_0,a_1,dots,a_{m-1})

    对于后面的所有(a_n)(a_n=sum_{i=1}^{k}f_ia_{n-i})

    (f)给出。

    暴力?

    小学生:直接干。

    初中生:矩阵乘法。

    更好的做法

    现在设(A)为转移矩阵。列向量(vec v=(a_{k-1},a_{k-2},dots,a_0)),满足(Avec v=(a_{k},a_{k-1},dots,a_1))

    我们要求(A^nvec v)中的第(k)项。

    先介绍一些奇奇怪怪的东西(以下的东西不带严谨证明)

    对于一个矩阵(A),如果能找到一个数(lambda)和一个列向量(overline v),满足(lambdavec v=Avec v),那么分别称它们是一组特征值特征向量

    移项得((lambda E-A)vec v=0)

    这个方程满足当且仅当(det(lambda E -A)=0)(不会证)

    [lambda E-A= left( { egin{matrix} lambda-a_1 & -a_2 & cdots &-a_{m-1} & -a_m \ -1 & lambda & cdots & 0 &0 \ 0 & -1 &cdots & 0 & 0\ vdots & vdots & ddots & vdots &vdots \ 0 & 0 & cdots & -1 & lambda end{matrix} } ight)]

    用点基本的数学知识可以算出(det(lambda E-A)=lambda^k-sum_{i=0}^{k-1}a_{k-i}lambda^i)

    我们可以将其看作一个关于(lambda)的多项式(叫特征多项式),记作(g(lambda))

    Cayley-Hamilton定理(g(A)=O)(O)表示全(0)的矩阵)

    (感性理解:(AE-A=0),所以这个东西成立。当然这样证明显然是不对的)

    (A)替代(lambda),可以看做一个矩阵的多项式。并且显然每一项底数为(A)的两个多项式相乘符合交换律。

    考虑除法,形如(A^n=p(A)g(A)+q(A))

    由于(g(A)=O),所以(A^n=q(A))

    现在我们就真把它当多项式来算。我们计算(q(x)equiv x^n pmod {g(x)})

    用倍增+多项式取模算出(q(x))的每一项的系数。

    现在我们要求((A^nvec v)_m),即(sum_{i=0}^{k-1}q_i(A^ivec v)_m=sum_{i=0}^{k-1}q_ia_i)

    于是做完了。总时间复杂度(O(klg k lg n))

    using namespace std;
    #include <cstdio>
    #include <cstring>
    #include <algorithm>
    #define N 65536
    #define mo 998244353
    #define ll long long
    ll qpow(ll x,ll y=mo-2){
    	ll r=1;
    	for (;y;y>>=1,x=x*x%mo)
    		if (y&1)
    			r=r*x%mo;
    	return r;
    }
    int n,k;
    int a[N],f[N];
    int nN,re[N];
    void setlen(int n){
    	int bit=0;
    	for (nN=1;nN<=n;nN<<=1,++bit);
    	for (int i=1;i<nN;++i)
    		re[i]=re[i>>1]>>1|(i&1)<<bit-1;
    }
    void clear(int A[],int n){memset(A,0,sizeof(int)*n);}
    void copy(int A[],int a[],int n){
    	clear(A,nN);
    	for (int i=0;i<=n;++i)
    		A[i]=a[i];
    }
    void dft(int A[],int flag){
    	for (int i=0;i<nN;++i)
    		if (i<re[i])
    			swap(A[i],A[re[i]]);
    	static int wnk[N];
    	for (int i=1;i<nN;i<<=1){
    		ll wn=qpow(3,flag==1?(mo-1)/(2*i):mo-1-(mo-1)/(2*i));
    		wnk[0]=1;
    		for (int k=1;k<i;++k)
    			wnk[k]=wnk[k-1]*wn%mo;
    		for (int j=0;j<nN;j+=i<<1)
    			for (int k=0;k<i;++k){
    				ll x=A[j+k],y=(ll)A[j+k+i]*wnk[k];
    				A[j+k]=(x+y)%mo;
    				A[j+k+i]=(x-y)%mo;
    			}
    	}
    	if (flag==-1)
    		for (int i=0,invn=qpow(nN);i<nN;++i)
    			A[i]=(ll)A[i]*invn%mo;
    	for (int i=0;i<nN;++i)
    		A[i]=(A[i]+mo)%mo;
    }
    void multi(int c[],int a[],int b[],int n,int an=-1,int bn=-1){
    	if (an==-1) an=n-1;
    	if (bn==-1) bn=n-1;
    	static int A[N],B[N],C[N];
    	setlen(an+bn);
    	copy(A,a,an),dft(A,1);
    	if (a==b)
    		for (int i=0;i<nN;++i)
    			C[i]=(ll)A[i]*A[i]%mo;
    	else{
    		copy(B,b,bn),dft(B,1);
    		for (int i=0;i<nN;++i)
    			C[i]=(ll)A[i]*B[i]%mo;
    	}
    	dft(C,-1);
    	for (int i=0;i<=min(n-1,an+bn);++i)
    		c[i]=C[i];
    }
    void getinv(int c[],int a[],int n,int an=-1){
    	if (an==-1) an=n-1;
    	static int b[N],g[N];
    	int nn=1;for (;nn<n;nn<<=1);
    	clear(b,nn);
    	b[0]=qpow(a[0]);
    	for (int i=1;i<n;i<<=1){
    		clear(g,i<<1);
    		multi(g,b,b,2*i,i-1,i-1);
    		multi(g,g,a,2*i,2*i-1,min(2*i-1,an));
    		for (int j=0;j<2*i;++j)
    			b[j]=(2*b[j]-g[j]+mo)%mo;
    	}
    	for (int i=0;i<n;++i)
    		c[i]=b[i];
    }
    void getrev(int A[],int a[],int n){for (int i=0;i<=n;++i) A[i]=a[n-i];}
    void getdiv(int c[],int a[],int b[],int an,int bn){
    	static int A[N],B[N],C[N];
    	clear(A,an-bn+1),clear(B,an-bn+1);
    	getrev(A,a,an),getrev(B,b,bn);
    	getinv(B,B,an-bn+1,an-bn+1);
    	multi(C,A,B,an-bn+1,an-bn,an-bn);
    	getrev(c,C,an-bn);	
    }
    void getmod(int c[],int a[],int b[],int an,int bn){
    	static int D[N];
    	getdiv(D,a,b,an,bn);
    	multi(D,D,b,an+1,an-bn,bn);
    	for (int i=0;i<bn;++i)
    		c[i]=(a[i]-D[i]+mo)%mo;
    }
    int g[N],q[N],mx;
    void dfs(int n){
    	if (n==0){
    		q[0]=1;
    		mx=0;
    		return;
    	}
    	if (n&1){
    		dfs(n-1);
    		for (int i=mx;i>=0;--i)
    			q[i+1]=q[i];
    		q[0]=0;
    		if (mx+1<k)
    			mx++;
    		else{
    			getmod(q,q,g,mx+1,k);
    			mx=k-1;
    		}
    	}
    	else{
    		dfs(n>>1);
    		multi(q,q,q,2*mx+1,mx,mx);
    		if (2*mx<k)
    			mx*=2;
    		else{
    			getmod(q,q,g,2*mx,k);
    			mx=k-1;
    		}
    	}
    }
    int main(){
    //	freopen("in.txt","r",stdin); 
    	scanf("%d%d",&n,&k);
    	for (int i=1;i<=k;++i)
    		scanf("%d",&f[i]);
    	for (int i=0;i<k;++i)
    		scanf("%d",&a[i]),(a[i]+=mo)%=mo;
    	for (int i=0;i<k;++i)
    		g[i]=(mo-f[k-i])%mo;
    	g[k]=1;
    	dfs(n);
    	ll ans=0;
    	for (int i=0;i<k;++i)
    		(ans+=(ll)q[i]*a[i])%=mo;
    	printf("%lld
    ",ans);
    	return 0;
    }
    
  • 相关阅读:
    ReactJS入门学习一
    js控制html5 【video】标签中视频的播放和停止
    CentOS中vsftp安装与配置
    linux 添加多个网段
    js图片预加载后触发操作
    nodejs在后台运行
    asp.net环境搭建
    aspx aspx.cs
    linux 添加静态ip dns
    kali ssh服务开启登录
  • 原文地址:https://www.cnblogs.com/jz-597/p/13436814.html
Copyright © 2011-2022 走看看