常系数齐次线性递推

最新推荐文章于 2024-08-06 21:39:55 发布

A1847225889

最新推荐文章于 2024-08-06 21:39:55 发布

阅读量151

点赞数

分类专栏：数学

本文链接：https://blog.csdn.net/A1847225889/article/details/107802751

版权

数学专栏收录该内容

53 篇文章 0 订阅

订阅专栏

常系数齐次线性递推

参考资料

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)$ ，满足 $A\vec v=(a_{k},a_{k-1},\dots,a_1)$

我们要求 $A^n\vec v$ 中的第 $k$ 项。

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

对于一个矩阵 $A$ ，如果能找到一个数 $\lambda$ 和一个列向量 $\overline v$ ，满足 $\lambda\vec v=A\vec v$ ，那么分别称它们是一组特征值和特征向量。

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

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

$\lambda E-A= \left( { \begin{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} } \right)$

用点基本的数学知识可以算出 $\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$ 的矩阵）

（感性理解： $A E - 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^n\vec v)_m$ ，即 $\sum_{i=0}^{k-1}q_i(A^i\vec v)_m=\sum_{i=0}^{k-1}q_ia_i$

于是做完了。总时间复杂度 $O(k\lg 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\n",ans);
	return 0;
}