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

问题引入

给出数列 g g g，满足当 n > m n>m n>m时
g n = ∑ i = 1 m g n − i × a i g_n=\sum\limits_{i=1}^{m}g_{n-i}\times a_i gn=i=1∑mgn−i×ai
当 n < = m n<=m n<=m时, g n = c n g_n=c_n gn=cn

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

Newbie的解法

暴力递推计算

时间复杂度 O ( n m ) O(nm) O(nm)

Pupil的解法

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

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

Master的解法

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

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

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

若对于m阶矩阵A，有常数 λ \lambda λ，非零列向量 v ⃗ \vec v v ，满足 λ v ⃗ = A v ⃗ \lambda\vec v=A\vec v λv =Av 则称 λ \lambda λ为矩阵A的特征值， v ⃗ \vec v v 为矩阵的特征向量

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

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

Part 2 Hamilton-Cayley theorem

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

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

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

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

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

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

它长这样：
(1) λ I − A = ( λ − a 1 − a 2 ⋯ − a m − 1 − a m − 1 λ ⋯ 0 0 0 − 1 ⋯ 0 0 ⋮ ⋮ ⋱ ⋮ ⋮ 0 0 ⋯ − 1 λ ) \lambda I-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} \tag{1} } \right) λI−A=⎝⎜⎜⎜⎜⎜⎛λ−a1−10⋮0−a2λ−1⋮0⋯⋯⋯⋱⋯−am−100⋮−1−am00⋮λ⎠⎟⎟⎟⎟⎟⎞(1)

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

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

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

Part 4 计算答案

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

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

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

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

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

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

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

也就是 ∑ i = 0 m − 1 d i A i G \sum\limits_{i=0}^{m-1}d_iA^iG i=0∑m−1diAiG的第m行第一列
然而 A i G A^iG AiG的第m行第一列的值就是 g i + 1 g_{i+1} gi+1

所以 g n = ∑ i = 0 m − 1 d i g i + 1 = ∑ i = 0 m − 1 d i c i + 1 g_n=\sum\limits_{i=0}^{m-1}d_ig_{i+1}=\sum\limits_{i=0}^{m-1}d_ic_{i+1} gn=i=0∑m−1digi+1=i=0∑m−1dici+1

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

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

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

移项，可以得到 G ( x ) = r 1 − A ( x ) G(x)={r\over 1-A(x)} G(x)=1−A(x)r
在模 x m + 1 x^{m+1} xm+1意义下多项式求逆即可
时间复杂度是 O ( m log ⁡ m ) O(m\log m) O(mlogm)的

模板题（[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);
}