常系数齐次线性递推初探

2018-09-18 15:03:40


同步于 这里

简介

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

设有数列$\{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 (要求$\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);
}