常系数齐次线性递推初探
Zhang_RQ
2018-09-18 15:03:40
同步于 [这里](https://zhang-rq.github.io/2018/06/28/%E5%B8%B8%E7%B3%BB%E6%95%B0%E9%BD%90%E6%AC%A1%E7%BA%BF%E6%80%A7%E9%80%92%E6%8E%A8%E5%AD%A6%E4%B9%A0%E7%AC%94%E8%AE%B0/)
# 简介
什么是常系数齐次线性递推？
设有数列$\{a_0,a_1,a_2 \cdots \}$满足递推关系$a_n=\sum\limits_{i=1}^{k}a_{n-i}*f_i$，则称该数列满足k阶齐次线性递推关系。
# 矩阵乘法
这个很普及的东西还是提一下吧。
假设我们有一个满足k阶齐次线性递推关系的数列$ {a_0,a_1,a_2 \cdots}$，它满足的齐次线性递推关系为$a_n=\sum\limits_{i=1}^{k}a_{n-i}*f_i$，现在我们要求$a_n$。
假设我们现在维护着一个列矩阵，它的行数为k。
$$A=\begin{bmatrix}a_{n} \\a_{n-1} \\ \cdots \\ a_{n-k+2} \\ a_{n-k+1} \end{bmatrix}$$
我们考虑如何让这个矩阵的所有元素都向前递推一格，不难设计出k行k列的转移矩阵。
$$M=\begin{bmatrix}f_1 &f_2 &f_3 &f_4 & \cdots &f_{k-2} &f_{k-1} \\ 1 &0 &0 &0 & \cdots &0 &0 \\ 0 &1 &0 &0 & \cdots &0 &0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots\\ 0 &0 &0 &0 & \cdots &1 &0\end{bmatrix}$$
我们可以得到
$$\begin{bmatrix} f_1 &f_2 &f_3 &f_4 & \cdots &f_{k-2} &f_{k-1} \\ 1 &0 &0 &0 & \cdots &0 &0 \\ 0 &1 &0 &0 & \cdots &0 &0\\ \cdots & \cdots & \cdots & \cdots &\cdots & \cdots & \cdots\\ 0 &0 &0 &0 & \cdots &1 &0 \end{bmatrix} \times \begin{bmatrix} a_{n-1} \\ a_{n-2} \\ \cdots \\ a_{n-k+1} \\ a_{n-k}\end{bmatrix} =\begin{bmatrix} a_{n} \\ a_{n-1} \\ \cdots \\ a_{n-k+2} \\ a_{n-k+1} \end{bmatrix}$$
现在我们要求$a_n$，根据矩阵乘法的结合律，易得到：
$$\left( \begin{bmatrix} f_1 &f_2 &f_3 &f_4 & \cdots &f_{k-2} &f_{k-1} \\ 1 &0 &0 &0 & \cdots &0 &0 \\ 0 &1 &0 &0 & \cdots &0 &0\\ \cdots & \cdots& \cdots & \cdots & \cdots & \cdots & \cdots\\ 0 &0 &0 &0 & \cdots &1 &0 \end{bmatrix} \right) ^n \times \begin{bmatrix} a_{k-1} \\ a_{k-2} \\ \cdots \\ a_{1} \\ a_{0}\end{bmatrix} =\begin{bmatrix} a_{n+k-1} \\ a_{n+k-2} \\ \cdots \\ a_{n+1} \\ a_{n}\end{bmatrix}$$
所以我们只需要算出$M^n \times A$，然后取最后一个数就可以了。
可以使用矩阵快速幂，复杂度$\mathcal{O}(k^3log_2n)$
# 特征多项式
## 特征值，特征向量:
若有常数$\lambda$,向量$\vec v$，满足$\lambda\vec v=A\vec v$ ，则称向量$\vec v$为矩阵$A$的一组特征向量，$\lambda$为矩阵$A$的一组特征值。
**秩为k**的矩阵有**k组线性不相关**的特征向量。
## 特征多项式
对关系式进行变换：$(\lambda I-A)\vec v=0$。
这个等式有解的充要条件是$det(\lambda I-A)=0$,大致可以看做向量被拍扁之类的。
可以得到一个$k$次多项式，这个多项式的值等于$0$时把这个方程称为矩阵$A$的**特征方程**。这个多项式称为矩阵$A$的**特征多项式**。
特征多项式记为$f(x)=det(\lambda I-A)$。
其中，$det()$为行列式函数，$I$为单位矩阵。
$k$个解自然就是$k$个特征值。（并不需要解出来，怎么处理后面会说）
跟据算数基本定理，最终的多项式有$k$个解，可以写作$f(x)=\prod_i(\lambda_i-x)$。
## Cayley-Hamilton定理
根据Cayley-Hamilton定理，可知$f(A)=O$，$O$为0矩阵。
下面给出一个简单证明：
$f(A)=\prod\limits_{i}(\lambda_{i} I - A)$
考虑将$f(A)$得到的矩阵分别乘特征向量，如果最后都得到了$0$矩阵，因为这几个特征向量线性不相关，那么可证明$f(A)$乘以任意向量都会得到$0$矩阵，从而$f(A)$为$0$矩阵。
现在问题转换为证明$f(A)$乘任意特征向量都会得到$0$矩阵。
先证明：$(\lambda_i I-A) \times (\lambda_j I -A) = (\lambda_j I -A) \times (\lambda_i I -A)$
展开即可，不再赘述。
现在考虑任意一个特征向量 $\vec v_i$ ,$$\begin{aligned} f(A) \times \vec v_{i} = \prod _{j!=i} \times (\lambda_{j} I - A) \times (\lambda_i I -A) \times \vec v_{i} \\ \text{由特征向量和特征值的定义可知，} \\ (\lambda_i I -A) \times \vec v_{i}=0 \\ \therefore f(A) \times \vec v_{i} =O \end{aligned} $$
证毕。
## 优化递推
设转移矩阵为$M$。
根据矩阵快速幂那套理论，我们只要求出来$M^n$就可以了。
我们考虑$M$的特征多项式$f(x)$，这是一个$k$次多项式。
我们将$M$带入，就会得到一个关于$M$的$k$次多项式。
我们可以将$M^n$分解为$M^n=f(M) \times g(M) +R(M)$。
由于$f(M)=O$，所以$M^n=R(M)$，$R(M)$是一个次数不超过$k-1$的多项式。
也就是说，我们只要求出来$M^n \% f(M)$ 就万事大吉了。
但是我们怎么求呢？
我们考虑快速幂的过程（就是倍增）。
假设我们现在已知 $g(M)=M^{2^i} \% f(M)$，现在要求$h(M)= M^{2^{i+1}} \% f(M)$。
一个直接的想法就是令$H(M)=g(M) \times g(M)$。
但是这样做，$H(M)$ 的次数是$2k-2$的。
那我们考虑原本的递推关系，$a_n=\sum\limits_{i=1}^{k}a_{n-i}*f_i$。
不难得到$M^n=\sum\limits_{i=1} ^{k} M^{n-i} \times f_{i}$。
所以我们可以用这个式子将多余出来的系数都向前压一位。
这样我们就得到了一个$\mathcal{O}(k^2\log_2 n)$的做法。
那有没有优化的余地呢？
我们从倍增的过程入手，可以发现$H(M)=g(M) \times g(M)$的过程可以由FFT加速至$\mathcal{O}(klog_2k)$。
现在只要解决压系数的过程就行了。
我们只要让$h(M) =H(M) \% f(M)$就行了。
等等，$f(M)$怎么求？
根据定义，$f(x)=det(xI - M)$，得到
$$f(x) = |x I - M| = \begin{bmatrix} x- a_1 & -a_2 & -a_3 & \cdots & -a_{k - 2} & -a_{k - 1} & -a_k \\ -1 & x & 0 & \cdots & 0 & 0 & 0 \\ 0 & -1 & x & \cdots & 0 & 0 & 0 \\ 0 & 0 & -1 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -1 & x & 0 \\ 0 & 0& 0 & \cdots & 0 & -1 & x\end{bmatrix}$$
我们对其进行第一列展开，得到
$$f(x) = (x - a_1)M_{11} + (-a_2)M_{12} + \cdots + (-a_k)M_{1n} = x ^ k - a_1 x ^ {k - 1} - a_2x ^ {k - 2} - \cdots - a_k$$
$M_{i,j}$代表$M$的算术余子式。
只要直接上多项式取模就行了。
最后的总复杂度$\mathcal{O}(k \log_2 k \log_2 n)$
我们手动模拟一下多项式取模的过程，其实可以发现我们手动向前压的过程就是在暴力取模。
# 例题
[BZOJ4161](https://www.lydsy.com/JudgeOnline/problem.php?id=4161) (要求$\mathcal{O}(k^2\log_2 n)$)
[洛谷P4732](https://www.luogu.org/problemnew/show/P4723) (要求$\mathcal{O}(k\log_2 k\log_2 n)$)
[洛谷P3824](https://www.luogu.org/problemnew/show/P3824) (要求$\mathcal{O}(k^2\log_2 n)$)
洛谷P4732的$\mathcal{O}(k\log_2 k\log_2 n)$的代码太长了，就不粘在这里了。感兴趣的同学可以来[这里](https://github.com/Zhang-RQ/OI_Code/blob/master/%E6%B4%9B%E8%B0%B7/4723.cpp)看。
附上[BZOJ4161](https://www.lydsy.com/JudgeOnline/problem.php?id=4161)的$\mathcal{O}(k^2\log_2 n)$代码
```cpp
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cctype>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<vector>
#include<set>
#include<map>
#include<queue>
#include<stack>
#include<cassert>
typedef long long ll;
typedef unsigned long long ull;
using namespace std;
const int P=1000000007;
const int MAXN=4010;
int n,k,ans;
int f[MAXN],h[MAXN];
struct Matrix{ //其实是多项式
int a[MAXN];
Matrix (){memset(a,0,sizeof a);}
int& operator [] (const int &i) {return a[i];}
int operator [] (const int &i) const {return a[i];}
inline Matrix operator * (const Matrix &rhs) const
{
Matrix ret;
for(int i=0;i<k;i++)
for(int j=0;j<k;j++)
(ret[i+j]+=1ll*a[i]*rhs[j]%P)%=P;
for(int i=2*k-2;i>=k;ret[i--]=0)
for(int j=1;j<=k;j++) //这里就是多项式取模优化的地方
(ret[i-j]+=1ll*ret[i]*f[j]%P)%=P; //可以认为是暴力向前压系数
return ret;
}
}res;
Matrix ksm(Matrix a,int b)
{
Matrix ret;
ret[0]=1;
for(;b;a=a*a,b>>=1) if(b&1) ret=ret*a;
return ret;
}
int main()
{
scanf("%d%d",&n,&k);
for(int i=1;i<=k;i++) scanf("%d",&f[i]),f[i]=f[i]>0?f[i]:f[i]+P;
for(int i=0;i<k;i++) scanf("%d",&h[i]),h[i]=h[i]>0?h[i]:h[i]+P;
if(n<k) printf("%d\n",h[n]);
res[1]=1;ans=0;
res=ksm(res,n);
for(int i=0;i<k;i++) ans=(ans+1ll*res[i]*h[i]%P)%P;
printf("%d\n",ans);
}
```