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

问题引入

给出数列$g$，满足当$n>m$时
$g_n=\sum\limits_{i=1}^{m}g_{n-i}\times a_i$
当$n<=m$时,$g_n=c_n$

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

Newbie的解法

暴力递推计算

时间复杂度$O(nm)$

Pupil的解法

可以将转移和数列都写成$m\times m$的矩阵的形式，矩阵快速幂即可

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

Master的解法

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

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

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

若对于m阶矩阵A，有常数$\lambda$，非零列向量$\vec v$，满足$\lambda\vec v=A\vec 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\times A^y$，由于前后都是A，矩阵乘法满足结合律，因此指数可以任意分配，换成$A^y \times A^x$也是可以的

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

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

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

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

它长这样：
$\lambda I-A= \left( { \begin{matrix} \lambda-a_1 &amp; -a_2 &amp; \cdots &amp;-a_{m-1} &amp; -a_m \\ -1 &amp; \lambda &amp; \cdots &amp; 0 &amp;0 \\ 0 &amp; -1 &amp;\cdots &amp; 0 &amp; 0\\ \vdots &amp; \vdots &amp; \ddots &amp; \vdots &amp;\vdots \\ 0 &amp; 0 &amp; \cdots &amp; -1 &amp; \lambda \end{matrix} \tag{1} } \right)$

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

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

化简整理，可得$f(\lambda)=|\lambda I-A|=\lambda^m-\sum\limits_{i=0}^{m-1}a_{m-i}\lambda ^i$
负号都被行列式里面逆序对个数的负号消掉了。

Part 4 计算答案

我们设要求的数列$g$的初始矩阵为$G$，它是一个m行1列的矩阵（列向量），从第m行到第1行分别为$g_{1\dots 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)\times G$与$A^{n-1}\times G$的结果是一样的。

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

也就是$\sum\limits_{i=0}^{m-1}d_iA^iG$的第m行第一列
然而$A^iG$的第m行第一列的值就是$g_{i+1}$

所以$g_n=\sum\limits_{i=0}^{m-1}d_ig_{i+1}=\sum\limits_{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)={r\over 1-A(x)}$
在模$x^{m+1}$意义下多项式求逆即可
时间复杂度是$O(m\log 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\n",ans);
}