常系数齐次线性递推学习笔记

定义

对于数列$f$，如果有递推式

$$f_n=\sum _{i=1}^k{a}_i\times f_{n-i}(n\ge k)$$

并且已知前$k$项，可称其为$k$ 阶齐次线性递推数列。

这个算法可以快速求出$f$的第$n$项，复杂度为$O(k^2\log n)$,可优化至$O(k\log k\log n)$

前置知识

矩阵快速幂

不难看出矩阵快速幂可以解决这个问题

我们构造$k$转移矩阵

$$A=\left(\begin{array}{cccccc}{a}_1& a_2& a_3& \cdots & a_{k-1}& a_k\\ 1& 0& 0& \cdots & 0& 0\\ 0& 1& 0& \cdots & 0& 0\\ 0& 0& 1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & 1& 0\end{array}\right)$$

对于一个状态

$$S=\left(\begin{array}{c}{f}_{i+k-1}\\ f_{i+k-2}\\ ⋮\\ ⋮\\ f_{i+1}\\ f_i\end{array}\right)$$

左乘上转移矩阵即可得到下一个状态

$$AS=\left(\begin{array}{c}{f}_{i+k}\\ f_{i+k-1}\\ ⋮\\ ⋮\\ f_{i+2}\\ f_{i+1}\end{array}\right)$$

我们只需要求出$A^{n}S$就可以算出答案

直接求的复杂度是$O(k^3\log n)$，不可接受

接下来将对矩阵快速幂进行优化

特征向量、值、多项式、方程

对于一个$k$阶方阵$A$，如果有列向量$x$和实数$\lambda$，满足

$$Ax=\lambda x$$

我们称$x$为$A$的特征向量，$\lambda$为$A$对于$x$的特征值。

对这个式子变形

$$(\lambda I-A)x=0$$

其中$I$为$k$阶单位矩阵

鬼知道我为什么要写这个

我们记关于$\lambda$的多项式

$$f(\lambda )=|\lambda I-A|$$

为$A$的特征多项式

$f(\lambda )=0$称为特征方程，不过我们用不上

Caylay-Hamilton 定理

定理内容：对于方阵$A$,设$f(\lambda )$为其特征多项式，有$f(A)=0$

你没看错，$f(A)$是矩阵的多项式，就是把矩阵$A$代进去，计算$a_{0}I+a_{1}A+a_2{A}^2+\cdots +a_k{A}^k$得到一个零矩阵

证明有太多前置知识，我还不会，希望有生之年能补上。

代数余子式

余子式：一个矩阵$A$去掉$A_{i,j}$所在的第$i$行第$j$列后的行列式称为余子式（对就是矩阵树那样），记为$a_{i,j}$,是个具体的数。

代数余子式：$(-1{)}^{i+j}{a}_{i,j}$

引理 矩阵的行列式等于它第一行的代数余子式之和。

对于第$1$行第$i$个元素，因为是第一行，所以贡献的逆序对个数为$i-1$,刚好和代数余子式的$(-1{)}^{i+j}$消掉，去掉第$1$行第$i$列之后顺序没有任何影响，所以乘一个$a_{i,j}$即可

算法流程

考虑之前的矩阵快速幂，最后推到了$A^{n}S$

我们考虑怎么快速求出$A^n$

我们发现$A^n$是一个矩阵的多项式，考虑前面的CH定理

我们设$A$的特征多项式为$f(\lambda )$

那么我们可以把$A^n$写成

$$A^n=f(A)g(A)+h(A)$$

因为$f(A)=0$,所以

$$A^n=h(A)$$

也就是说，$A^n$和$h(A)$在矩阵意义上等价

那么如果我们知道$A$的特征多项式，我们就可以用一般意义下的多项式取模（具体做法后面讲）算出$h(A)$的表达式,并且最高项次数不超过$k-1$

我们设求出的

$$h(x)=\sum _{i=0}^{k-1}{b}_i{x}^i$$

考虑我们要求的

$$A^{n}S$$

因为$A^n=h(A)$,所以等于

$$\sum _{i=0}^{k-1}{b}_i{A}^{i}S$$

$A^{i}S$相当于是第$i+1$个状态，我们写成矩阵形式

$$\sum _{i=0}^{k-1}{b}_i\left(\begin{array}{c}{f}_{i+k-1}\\ f_{i+k-2}\\ ⋮\\ ⋮\\ f_{i+1}\\ f_i\end{array}\right)=A^{n}S=\left(\begin{array}{c}{f}_{n+k-1}\\ f_{n+k-2}\\ ⋮\\ ⋮\\ f_{n+1}\\ f_n\end{array}\right)$$

仔细观察，发现我们只需要最下面的元素

所以

$$\sum _{i=0}^{k-1}{b}_i{f}_i=f_n$$

也就是把$b$和题目给的前$k$项对应乘起来就好了

---

现在的问题是如何求特征多项式

$$\lambda I-A=\left(\begin{array}{ccccccc}\lambda -a_1& -a_2& -a_3& -a_4& \cdots & -a_{k-1}& -a_k\\ -1& \lambda & 0& 0& \cdots & 0& 0\\ 0& -1& \lambda & 0& \cdots & 0& 0\\ 0& 0& -1& \lambda & \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \cdots & \lambda & 0\\ 0& 0& 0& 0& \cdots & -1& \lambda \end{array}\right)$$

根据上面的引理，

$$|\lambda I-A|=(\lambda -a_1)a_{1,1}+\sum _{i=2}^k(-1{)}^{i+1}(-a_i)a_{1,i}$$

其中$a_{i,j}$为余子式

因为去掉第$1$行和任意一列后剩下的是一个下三角矩阵，并且对角线上删除的列的左边是$-1$,右边是$\lambda$,所以

$$a_{1,i}=(-1{)}^{i-1}{\lambda }^{k-i}$$

所以

$$\begin{gathered} |\lambda I-A|=(\lambda -a_1){\lambda }^{k-1}+\sum _{i=2}^k{a}_i{\lambda }^{k-i} \\ ={\lambda }^k-\sum _{i=1}^k{a}_i{\lambda }^{k-i} \end{gathered}$$

即

$$f(\lambda )={\lambda }^k-\sum _{i=1}^k{a}_i{\lambda }^{k-i}$$

就是特征多项式

实现

前面的多项式取模部分是挖了坑的

$$A^n=f(A)g(A)+h(A)$$

我们已知$f(A)$为特征多项式，要求出余数$h(A)$

因为$n$很大，所以我们倍增计算并顺便取模

常用的暴力取模好写的多，复杂度$O(k^2\log n)$，具体实现见代码注释

BZOJ4161 Shlw loves matrixI

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#define MAXN 4005//数组开两倍
using namespace std;
int n,k,F[MAXN],R[MAXN],H[MAXN];
const int MOD=1e9+7;
typedef long long ll;
inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;}
inline void mul(int* f,int* g,int* to)
{
	static int t[MAXN];
	memset(t,0,sizeof(t));
	for (int i=0;i<k;i++)
		for (int j=0;j<k;j++)
			t[i+j]=add(t[i+j],(ll)f[i]*g[j]%MOD);//求卷积
	for (int i=(k<<1);i>=k;t[i--]=0)
		for (int j=1;j<=k;j++) 
			t[i-j]=add(t[i-j],(ll)t[i]*F[j]%MOD);//因为最高项次数为1，后面都是读入的数的相反数，所以读入不取负，这里的减改成加即可
	for (int i=0;i<k;i++) to[i]=t[i];
}
inline void qpow(int* f,int p)
{
	static int ans[MAXN];
	memset(ans,0,sizeof(ans));
	ans[0]=1;
	while (p)
	{
		if (p&1) mul(f,ans,ans);
		mul(f,f,f);p>>=1;
	}
	for (int i=0;i<k;i++) f[i]=ans[i];
}
int main()
{
	scanf("%d%d",&n,&k);
	for (int i=1;i<=k;i++) scanf("%d",&F[i]),F[i]<0&&(F[i]+=MOD);
	for (int i=0;i<k;i++) scanf("%d",&H[i]),H[i]<0&&(H[i]+=MOD);
	R[1]=1;
	qpow(R,n);
	int ans=0;
	for (int i=0;i<k;i++) ans=add(ans,(ll)R[i]*H[i]%MOD);
	printf("%d\n",ans);
	return 0;
}

---

如果不幸遇到了毒瘤出题人，则需要用多项式全家桶中的多项式取模

复杂度为$O(k\log k\log n)$,常数较大

Luogu P4723

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#include <algorithm>
#define MAXN 131072
using namespace std;
inline int read()
{
	int ans=0,f=1;
	char c=getchar();
	while (!isdigit(c)) (c=='-')&&(f=-1),c=getchar();
	while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();
	return f*ans;
}
const int MOD=998244353;
typedef long long ll;
inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;}
inline int dec(const int& x,const int& y){return x<y? x-y+MOD:x-y;}
inline int qpow(int a,int p)
{
	int ans=1;
	while (p)
	{
		if (p&1) ans=(ll)ans*a%MOD;
		a=(ll)a*a%MOD;p>>=1;
	}
	return ans;
}
#define inv(x) qpow(x,MOD-2)
int rt[2][24];
int l,lim,r[MAXN];
inline void init(){lim=1<<l;for (int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));}
inline void NTT(int* a,int type)
{
	for (int i=0;i<lim;i++) if (i<r[i]) swap(a[i],a[r[i]]);
	for (int L=0;L<l;L++)
	{
		int mid=1<<L,len=mid<<1,Wn=rt[type][L+1];
		for (int s=0;s<lim;s+=len)
			for (int k=0,w=1;k<mid;k++,w=(ll)w*Wn%MOD)
			{
				int x=a[s+k],y=(ll)w*a[s+mid+k]%MOD;
				a[s+k]=add(x,y);a[s+mid+k]=dec(x,y);
			}
	}
	if (type)
	{
		int t=inv(lim);
		for (int i=0;i<lim;i++) a[i]=(ll)a[i]*t%MOD;
	}
}
void getinv(int* A,int* B,int n)
{
	static int f[MAXN],t[MAXN];
	if (n==1) return (void)(*B=inv(*A));
	getinv(A,t,(n+1)>>1);
	for (l=0;(1<<l)<(n<<1);l++);
	init();
	for (int i=0;i<n;i++) f[i]=A[i];
	for (int i=n;i<lim;i++) f[i]=t[i]=0;
	NTT(f,0);NTT(t,0);
	for (int i=0;i<lim;i++) B[i]=(ll)t[i]*dec(2,(ll)f[i]*t[i]%MOD)%MOD;
	NTT(B,1);
	for (int i=n;i<lim;i++) B[i]=0;
}
void getmod(int* A,int* B,int n,int m)
{
	static int f[MAXN],g[MAXN],h[MAXN],t[MAXN];
	for (int i=0;i<=n;i++) f[i]=A[n-i];
	for (int i=0;i<=m;i++) g[i]=B[m-i];
	getinv(g,t,n-m+1);
	for (int i=n-m+1;i<lim;i++) f[i]=0;
	NTT(f,0);NTT(t,0);
	for (int i=0;i<lim;i++) h[i]=(ll)f[i]*t[i]%MOD;
	NTT(h,1);
	for (int i=n-m+1;i<lim;i++) h[i]=0;
	reverse(h,h+n-m+1),reverse(g,g+m+1);
	for (l=0;(1<<l)<=n;l++);
	init();
	for (int i=m+1;i<lim;i++) g[i]=0;
	for (int i=n-m+1;i<lim;i++) h[i]=0;
	NTT(g,0);NTT(h,0);
	for (int i=0;i<lim;i++) g[i]=(ll)g[i]*h[i]%MOD;
	NTT(g,1);
	for (int i=0;i<m;i++) A[i]=dec(A[i],g[i]);
	for (int i=m;i<lim;i++) A[i]=0;
}
int n,k;
inline void mul(int* A,int* B,int* C,int* M)
{
	static int f[MAXN],g[MAXN];
	for (int i=0;i<k;i++) f[i]=A[i],g[i]=B[i];
	for (l=0;(1<<l)<(k<<1);l++);
	init();
	for (int i=k;i<lim;i++) f[i]=g[i]=0;
	NTT(f,0);NTT(g,0);
	for (int i=0;i<lim;i++) C[i]=(ll)f[i]*g[i]%MOD;
	NTT(C,1);
	getmod(C,M,(k-1)<<1,k);
}
inline void qpow(int* A,int* R,int p)
{
	static int ans[MAXN];
	memset(ans,0,sizeof(ans));
	ans[0]=1;
	while (p)
	{
		if (p&1) mul(ans,A,ans,R);
		mul(A,A,A,R);
		p>>=1;
	}
	for (int i=0;i<k;i++) A[i]=ans[i];
}
int f[MAXN],a[MAXN],h[MAXN];
int main()
{
	rt[0][23]=qpow(3,119);rt[1][23]=inv(rt[0][23]);
	for (int i=23;i>=1;i--) rt[0][i-1]=(ll)rt[0][i]*rt[0][i]%MOD,rt[1][i-1]=(ll)rt[1][i]*rt[1][i]%MOD;
	n=read(),k=read();
	for (int i=1;i<=k;i++) (f[k-i]=-read())<0&&(f[k-i]+=MOD);
	for (int i=0;i<k;i++) (a[i]=read())<0&&(a[i]+=MOD);
	f[k]=1;
	h[1]=1;
	qpow(h,f,n);
	int ans=0;
	for (int i=0;i<k;i++) ans=add(ans,(ll)h[i]*a[i]%MOD);
	printf("%d\n",ans);
	return 0;
}