常系数齐次线性递推快速算法学习笔记

今天集训被线代狠狠的虐了一发。

不过还有一点收获的，比如这个。

---

数列 (f) 满足 (f_n=sumlimits_{i=1}^ka_if_{n-i}(nge k))，其中 (a_1dots a_k,f_0dots f_{k-1}) 均给出。求 (f_n)。

(nle 10^9,kle 30000)。

---

先要弄懂一些基本定义，

矩阵，行列式，高斯消元这些基本的东西就自己看别的东西去吧，我也不知道有什么好的资料，不妨去洛谷搜模板看题解，那里面的都不错。

然后讲一下特征值，特征多项式和 Hamilton-Cayley 定理，就能做这题了。

对于 (n imes n) 的矩阵 (A)，如果存在数 (lambda) 和非零列向量 (x) 满足 (Ax=lambda x)，那么 (lambda) 是 (A) 的特征值，(x) 是 (A) 的特征向量。

那么 (Ax=lambda I x)。((lambda I-A)x=0)。因为 (x) 不为零向量，所以 (det(lambda I-A)=0)，也就是 (lambda I-A) 不满秩。

我们称 (det(lambda I-A)=0) 为 (A) 的特征多项式，元是 (lambda)。特征多项式是 (n) 次的，他的 (n) 个根就是 (A) 的所有特征值。（可能有相等的根）

对于上三角矩阵，所有的特征值就是主对角线上的所有值。

如果有 (n) 个线性无关的特征向量 (x_i)（当且仅当 (A) 满秩），那么有 (Aegin{bmatrix}x_1&x_2&cdots&x_nend{bmatrix}=egin{bmatrix}x_1&x_2&cdots&x_nend{bmatrix}egin{bmatrix}lambda_1&0&0&cdots&0\0&lambda_2&0&cdots&0\cdots\0&0&0&cdots&lambda_nend{bmatrix})。

Hamilton-Cayley 定理：对于矩阵 (A) 的特征多项式 (f(lambda)=sumlimits_{i=0}^nc_ilambda^i)，有 (f(A)=0)，即 (sumlimits_{i=0}^nc_iA^i=0)。

会了这些，就可以开始了。

---

(O(k^3log n)) 的相信大家都会。（什么？不会？赶快去学矩阵快速幂）

首先考虑写出转移矩阵 (A) 和初始行向量 (f)，我们要求的是 (f imes A^n) 的第 (0) 维。

假如我们构造出了一个序列 (c) 使得

[A^n=sumlimits_{i=0}^{k-1}c_iA^i ]

那么有：

[f imes A^n=sumlimits_{i=0}^{k-1}c_i(f imes A^i) ]

[(f imes A^n)_0=sumlimits_{i=0}^{k-1}c_i(f imes A^i)_0 ]

[f_n=sumlimits_{i=0}^{k-1}c_if_i ]

那么就可以 (O(k)) 计算了。

那么 (c) 怎么弄呢？

令 (R(A)=sumlimits_{i=0}^{k-1}c_iA^i)

假如存在 (k) 次多项式 (G(A))，使得 (A^n=F(A)G(A)+R(A))。(标准多项式除法形式)

当 (G(A)=0) 时就有 (A^n=R(A))，所以要求的就是 (A^nmod G(A))，快速幂+多项式除法 (O(klog klog n)) 解决。

那么如何构造 (G(A)=0) 的多项式呢？

看到上面的 Hamilton-Cayley 定理，令 (G(lambda)=det(lambda I-A)) 即可。

手玩一下，发现 (det(lambda I-A)=-sumlimits_{i=0}^{k-1}lambda^ia_{k-i}+lambda^k)。（注意交换行的时候取相反数！！！）

所以 (g_i=-a_{k-i}(i e k),g_k=1)。

时间复杂度 (O(klog klog n))。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=333333,mod=998244353;
#define FOR(i,a,b) for(int i=(a);i<=(b);i++)
#define ROF(i,a,b) for(int i=(a);i>=(b);i--)
#define MEM(x,v) memset(x,v,sizeof(x))
inline int read(){
    int x=0,f=0;char ch=getchar();
    while(ch<'0' || ch>'9') f|=ch=='-',ch=getchar();
    while(ch>='0' && ch<='9') x=x*10+ch-'0',ch=getchar();
    return f?-x:x;
}
int n,k,s,f[maxn],g[maxn],fac[maxn],ans[maxn],prod[maxn],tmp[maxn],lim,l,rev[maxn],invt[maxn],Arev[maxn],Btmp[maxn],Brev[maxn],Brevinv[maxn],C[maxn];
inline int add(int x,int y){return x+y<mod?x+y:x+y-mod;}
inline int sub(int x,int y){return x<y?x-y+mod:x-y;}
inline int mul(int x,int y){return 1ll*x*y%mod;}
inline int qpow(int a,int b){
	int ans=1;
	for(;b;b>>=1,a=mul(a,a)) if(b&1) ans=mul(ans,a);
	return ans;
}
void init(int upr){
	for(lim=1,l=0;lim<upr;lim<<=1,l++);
	FOR(i,0,lim-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
}
void NTT(int *A,int tp){
	FOR(i,0,lim-1) if(i<rev[i]) swap(A[i],A[rev[i]]);
	for(int i=1;i<lim;i<<=1)
		for(int j=0,r=i<<1,Wn=qpow(3,mod-1+tp*(mod-1)/r);j<lim;j+=r)
			for(int k=0,w=1;k<i;k++,w=mul(w,Wn)){
				int x=A[j+k],y=mul(w,A[i+j+k]);
				A[j+k]=add(x,y);
				A[i+j+k]=sub(x,y);
			}
	if(tp==-1){
		int linv=qpow(lim,mod-2);
		FOR(i,0,lim-1) A[i]=mul(A[i],linv);
	}
}
void inv(int *A,int *B,int deg){
	if(deg==1) return void(B[0]=qpow(A[0],mod-2));
	inv(A,B,(deg+1)>>1);
	init(deg<<1);
	FOR(i,0,deg-1) invt[i]=A[i];
	FOR(i,deg,lim-1) invt[i]=0;
	NTT(invt,1);NTT(B,1);
	FOR(i,0,lim-1) B[i]=mul(B[i],sub(2,mul(invt[i],B[i])));
	NTT(B,-1);
	FOR(i,deg,lim-1) B[i]=0;
}
void division(int *A,int *B,int *D,int n,int m){
	if(n<m){
		init(n<<1);
		FOR(i,n+1,lim-1) D[i]=A[i]=0;
		FOR(i,0,n) D[i]=A[i];
		return;
	}
	init(n<<1);
	FOR(i,0,lim-1) Arev[i]=Brev[i]=Brevinv[i]=C[i]=Btmp[i]=D[i]=0;
	FOR(i,n+1,lim-1) A[i]=0;
	FOR(i,m+1,lim-1) B[i]=0;
	FOR(i,0,n) Arev[i]=A[n-i];
	FOR(i,0,n-m) Brev[i]=B[m-i];
	inv(Brev,Brevinv,n-m+1);
	init(n<<1);
	NTT(Arev,1);NTT(Brevinv,1);
	FOR(i,0,lim-1) C[i]=mul(Arev[i],Brevinv[i]);
	NTT(C,-1);
	FOR(i,0,(n-m)/2) swap(C[i],C[n-m-i]);
	FOR(i,n-m+1,lim-1) C[i]=0;
	FOR(i,0,m) Btmp[i]=B[i];
	NTT(Btmp,1);NTT(C,1);
	FOR(i,0,lim-1) D[i]=mul(Btmp[i],C[i]);
	NTT(D,-1);
	FOR(i,0,m-1) D[i]=sub(A[i],D[i]);
	FOR(i,m,lim-1) D[i]=0;
}
int main(){
	n=read();k=read();
	FOR(i,1,k) g[k-i]=(mod-read())%mod;
	FOR(i,0,k-1) f[i]=(read()+mod)%mod;
	g[k]=fac[1]=ans[0]=1;
	while(n){
		if(n&1){
			init(k<<1);
			NTT(fac,1);NTT(ans,1);
			FOR(i,0,lim-1) prod[i]=mul(fac[i],ans[i]);
			NTT(fac,-1);NTT(ans,-1);NTT(prod,-1);
			division(prod,g,ans,(k-1)<<1,k);
		}
		init(k<<1);
		NTT(fac,1);
		FOR(i,0,lim-1) prod[i]=mul(fac[i],fac[i]);
		NTT(fac,-1);NTT(prod,-1);
		division(prod,g,fac,(k-1)<<1,k);
		n>>=1;
	}
	FOR(i,0,k-1) s=add(s,mul(f[i],ans[i]));
	printf("%d
",s);
}