同步于 这里

简介

什么是常系数齐次线性递推？

设有数列\{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$个解，可以写作$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 只要直接上多项式取模就行了。 最后的总复杂度$\mathcal{O}(k \log_2 k \log_2 n)

我们手动模拟一下多项式取模的过程，其实可以发现我们手动向前压的过程就是在暴力取模。

例题

BZOJ4161 (要求\mathcal{O}(k^2\log_2 n))

洛谷P4732 (要求\mathcal{O}(k\log_2 k\log_2 n))

洛谷P3824 (要求\mathcal{O}(k^2\log_2 n))

洛谷P4732的\mathcal{O}(k\log_2 k\log_2 n)的代码太长了，就不粘在这里了。感兴趣的同学可以来这里看。

附上BZOJ4161的\mathcal{O}(k^2\log_2 n)代码

#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);
}