【学习小记】常系数齐次线性递推

问题引入

给出数列(g)，满足当(n>m)时

[g_n=sumlimits_{i=1}^{m}g_{n-i} imes a_i ]

当(n<=m)时,(g_n=c_n)

m比较小，n特别大，快速计算(g_n)

Newbie的解法

暴力递推计算

时间复杂度(O(nm))

Pupil的解法

可以将转移和数列都写成(m imes m)的矩阵的形式，矩阵快速幂即可

时间复杂度(O(m^3log n))

Master的解法

我们需要一些数学知识进行铺垫：

Part 1 矩阵的特征值与特征多项式

我们知道一个矩阵乘一个列向量仍然是一个列向量。

若对于m阶矩阵A，有常数(lambda)，非零列向量(vec v)，满足$$lambdavec v=Avec v$$则称(lambda)为矩阵A的特征值，(vec v)为矩阵的特征向量

上式也可以写作$$(lambda I-A)vec v=0$$其中(I)为单位矩阵
此式有解的充要条件是(|lambda I-A|=0)，即矩阵(lambda I-A)的行列式为0

(|lambda I-A|)可以看做是关于(lambda)的一个m次多项式，记作(f(lambda))，(f(lambda))称作矩阵A的特征多项式，对于矩阵A的任意一个特征值(lambda_0)，都有(f(lambda_0)=0)。

Part 2 Hamilton-Cayley theorem

对于矩阵，也一样的定义多项式运算（把多项式中的x换乘矩阵A），加法就是直接对应相加，常数乘法就按位相乘，乘法是矩阵乘法，0次方是单位矩阵，它的结果仍然是一个矩阵。

显然，矩阵多项式满足交换律，即(f(A)g(A)=g(A)f(A))成立。
简单证明：考虑某两项相乘的结果(A^x imes A^y)，由于前后都是A，矩阵乘法满足结合律，因此指数可以任意分配，换成(A^y imes A^x)也是可以的

哈密顿—凯莱定理：对于矩阵A的特征多项式(f(x))，满足(f(A)=0)

证明网上到处都有，此处就不赘述了。

Part 3 求解转移矩阵的特征多项式

回到原题，我们对于Pupil解法的转移矩阵A，求解它的特征多项式
考虑矩阵(lambda I-A)

它长这样：

[lambda I-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} ag{1} } ight) ]

根据行列式的定义，将第一行展开
(|lambda I-A|=(lambda-a_1)A_{1,1}+a_2 imes A_{1,2}+cdots+a_m imes A_{1,m})
其中(A{i,j})表示矩阵A的代数余子式，即挖掉第i行和第j列以后剩下的矩阵的行列式。

我们发现所有的余子矩阵都是下三角矩阵，行列式就是对角线乘积。

化简整理，可得(f(lambda)=|lambda I-A|=lambda^m-sumlimits_{i=0}^{m-1}a_{m-i}lambda ^i)

Part 4 计算答案

我们设要求的数列(g)的初始矩阵为(G)，它是一个m行1列的矩阵（列向量），从第m行到第1行分别为(g_{1dots m})（注意顺序是反的）
实际上我们想知道的(g_n)就是矩阵(A^{n-1}G)的第m行第一列的值。

此时的关键就是(A^{n-1})，因为(n-1)非常大，无法直接计算

然而根据前面的铺垫，我们有(f(A)=0)，(A^{n-1})我们可以看做只有一项的一个关于A的多项式

那么根据多项式除法相关知识，可以得到(A^{n-1}=P(A)f(A)+Q(A))，其中(Q(A))的次数是小于(f(A))的次数也就是小于m的，(Q(A))相当于多项式(A^{n-1})对多项式(f(A))取模

可能会有这样的疑问，(f(A)=0)怎么能作除数呢？
其实并不要紧，我们并不需要知道(f(A))的实际值，我们相当于将(A^{n-1})减去了若干个(f(A))，将次数降低了，而结果不变。

实现上来说，由于(f)的系数已知，我们可以先将式子里的矩阵A换成变量(x)，代入，利用多项式取模算出Q的系数，然后再将x换回A，这样得出来的Q的系数是相同的。并且计算(Q(A) imes G)与(A^{n-1} imes G)的结果是一样的。

为了求出(Q(x))的系数，我们可以采用快速幂的做法，初始(Q_0(x)=x^1)，然后不断的自己与自己相乘，乘完对多项式(f(x))取模
这一部分如果暴力取模，时间复杂度为(O(m^2log n))
如果采用NTT优化多项式取模，时间复杂度为(O(mlog mlog n))
这样求出了(Q(A))的系数，不妨设(Q(A)=sumlimits_{i=0}^{m-1}d_iA^i)
要求矩阵(Q(A) imes G)的第m行第一列的值

也就是(sumlimits_{i=0}^{m-1}d_iA^iG)的第m行第一列
然而(A^iG)的第m行第一列的值就是(g_{i+1})

所以$$g_n=sumlimits_{i=0}^{m-1}d_ig_{i+1}=sumlimits_{i=0}^{m-1}d_ic_{i+1}$$

还有一种情况，前m项并没有直接给出，也是通过递推得出的，暴力递推求前m项的复杂度是(O(m^2))的
考虑优化

考虑数列(g)的一般生成函数(G(x))（与矩阵G不同）
转移序列(a)的一般生成函数(A(x))

由于(G(x))是无限长的一个序列，我们可以得到(G(x)=G(x)A(x)+r)
其中(r)是一个常数，相当于第0项

移项，可以得到$$G(x)={rover 1-A(x)}$$
在模(x^{m+1})意义下多项式求逆即可
时间复杂度是(O(mlog m))的

模板题（[BZOJ4161] Shlw loves matrixI）

Code

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <iostream>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;++i)
#define fod(i,a,b) for(int i=a;i>=b;--i)
#define N 4005
#define mo 1000000007
#define LL long long
using namespace std;
LL f[N],g[N],h[N],s1[N],a[N],u1[N];
int n,m;
void mul(LL *x,LL *y,LL *z)
{
	fo(i,0,2*m-2) u1[i]=0;
	fo(i,0,m-1) fo(j,0,m-1) u1[i+j]=(u1[i+j]+x[i]*y[j])%mo;
	fod(i,2*m-2,m)
	{
		fo(j,0,m) u1[i-m+j]=(u1[i-m+j]-f[j]*u1[i])%mo; 
	}
	fo(i,0,m-1) z[i]=u1[i];
}
int main()
{
	cin>>n>>m;
	fo(i,1,m) scanf("%lld",&a[i]),f[m-i]=-a[i];
	f[m]=1;
	g[1]=1;
	s1[0]=1;
	for(int t=n;t;t>>=1)
	{
		if(t&1) mul(s1,g,s1);
		mul(g,g,g);
	}
	fo(i,0,m-1) scanf("%lld",&h[i]);
	LL ans=0;
	fo(i,0,m-1) ans=(ans+s1[i]*h[i]%mo+mo)%mo;
	printf("%lld
",ans);
}