专栏文章

常系数齐次线性递推初探

ZZhang_RQ
算法·理论参与者 30已保存评论 31

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
31 条
当前快照
1 份
快照标识符
@mhz5rnoy
此快照首次捕获于
2025/11/15 01:55
4 个月前
此快照最后确认于
2025/11/29 05:24
3 个月前
查看原文
同步于 这里

简介

什么是常系数齐次线性递推？
设有数列{a0,a1,a2⋯ }\{a_0,a_1,a_2 \cdots \}满足递推关系an=∑i=1kan−i∗fia_n=\sum\limits_{i=1}^{k}a_{n-i}*f_i，则称该数列满足k阶齐次线性递推关系。

矩阵乘法

这个很普及的东西还是提一下吧。
假设我们有一个满足k阶齐次线性递推关系的数列a0,a1,a2⋯ {a_0,a_1,a_2 \cdots}，它满足的齐次线性递推关系为an=∑i=1kan−i∗fia_n=\sum\limits_{i=1}^{k}a_{n-i}*f_i，现在我们要求ana_n。
假设我们现在维护着一个列矩阵，它的行数为k。
A=[anan−1⋯an−k+2an−k+1]A=\begin{bmatrix}a_{n} \\a_{n-1} \\ \cdots \\ a_{n-k+2} \\ a_{n-k+1} \end{bmatrix}
我们考虑如何让这个矩阵的所有元素都向前递推一格，不难设计出k行k列的转移矩阵。 M=[f1f2f3f4⋯fk−2fk−11000⋯000100⋯00⋯⋯⋯⋯⋯⋯⋯0000⋯10]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} 我们可以得到 [f1f2f3f4⋯fk−2fk−11000⋯000100⋯00⋯⋯⋯⋯⋯⋯⋯0000⋯10]×[an−1an−2⋯an−k+1an−k]=[anan−1⋯an−k+2an−k+1]\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} 现在我们要求ana_n，根据矩阵乘法的结合律，易得到： ([f1f2f3f4⋯fk−2fk−11000⋯000100⋯00⋯⋯⋯⋯⋯⋯⋯0000⋯10])n×[ak−1ak−2⋯a1a0]=[an+k−1an+k−2⋯an+1an]\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} 所以我们只需要算出Mn×AM^n \times A，然后取最后一个数就可以了。
可以使用矩阵快速幂，复杂度O(k3log2n)\mathcal{O}(k^3log_2n)

特征多项式

特征值，特征向量:

若有常数λ\lambda,向量v⃗\vec v，满足λv⃗=Av⃗\lambda\vec v=A\vec v ，则称向量v⃗\vec v为矩阵AA的一组特征向量，λ\lambda为矩阵AA的一组特征值。
秩为k的矩阵有k组线性不相关的特征向量。

特征多项式

对关系式进行变换：(λI−A)v⃗=0(\lambda I-A)\vec v=0。
这个等式有解的充要条件是det(λI−A)=0det(\lambda I-A)=0,大致可以看做向量被拍扁之类的。
可以得到一个kk次多项式，这个多项式的值等于00时把这个方程称为矩阵AA的特征方程。这个多项式称为矩阵AA的特征多项式。
特征多项式记为f(x)=det(λI−A)f(x)=det(\lambda I-A)。
其中，det()det()为行列式函数，II为单位矩阵。
kk个解自然就是kk个特征值。（并不需要解出来，怎么处理后面会说）
跟据算数基本定理，最终的多项式有kk个解，可以写作f(x)=∏i(λi−x)f(x)=\prod_i(\lambda_i-x)。

Cayley-Hamilton定理

根据Cayley-Hamilton定理，可知f(A)=Of(A)=O，OO为0矩阵。
下面给出一个简单证明：
f(A)=∏i(λiI−A)f(A)=\prod\limits_{i}(\lambda_{i} I - A)
考虑将f(A)f(A)得到的矩阵分别乘特征向量，如果最后都得到了00矩阵，因为这几个特征向量线性不相关，那么可证明f(A)f(A)乘以任意向量都会得到00矩阵，从而f(A)f(A)为00矩阵。
现在问题转换为证明f(A)f(A)乘任意特征向量都会得到00矩阵。
先证明：(λiI−A)×(λjI−A)=(λjI−A)×(λiI−A)(\lambda_i I-A) \times (\lambda_j I -A) = (\lambda_j I -A) \times (\lambda_i I -A)
展开即可，不再赘述。
现在考虑任意一个特征向量 v⃗i\vec v_i ,f(A)×v⃗i=∏j!=i×(λjI−A)×(λiI−A)×v⃗i由特征向量和特征值的定义可知，(λiI−A)×v⃗i=0∴f(A)×v⃗i=O\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} 证毕。

优化递推

设转移矩阵为MM。
根据矩阵快速幂那套理论，我们只要求出来MnM^n就可以了。
我们考虑MM的特征多项式f(x)f(x)，这是一个kk次多项式。
我们将MM带入，就会得到一个关于MM的kk次多项式。
我们可以将MnM^n分解为Mn=f(M)×g(M)+R(M)M^n=f(M) \times g(M) +R(M)。
由于f(M)=Of(M)=O，所以Mn=R(M)M^n=R(M)，R(M)R(M)是一个次数不超过k−1k-1的多项式。
也就是说，我们只要求出来Mn%f(M)M^n \% f(M) 就万事大吉了。
但是我们怎么求呢？
我们考虑快速幂的过程（就是倍增）。
假设我们现在已知 g(M)=M2i%f(M)g(M)=M^{2^i} \% f(M)，现在要求h(M)=M2i+1%f(M)h(M)= M^{2^{i+1}} \% f(M)。
一个直接的想法就是令H(M)=g(M)×g(M)H(M)=g(M) \times g(M)。
但是这样做，H(M)H(M) 的次数是2k−22k-2的。
那我们考虑原本的递推关系，an=∑i=1kan−i∗fia_n=\sum\limits_{i=1}^{k}a_{n-i}*f_i。
不难得到Mn=∑i=1kMn−i×fiM^n=\sum\limits_{i=1} ^{k} M^{n-i} \times f_{i}。
所以我们可以用这个式子将多余出来的系数都向前压一位。
这样我们就得到了一个O(k2log⁡2n)\mathcal{O}(k^2\log_2 n)的做法。
那有没有优化的余地呢？ 我们从倍增的过程入手，可以发现H(M)=g(M)×g(M)H(M)=g(M) \times g(M)的过程可以由FFT加速至O(klog2k)\mathcal{O}(klog_2k)。
现在只要解决压系数的过程就行了。
我们只要让h(M)=H(M)%f(M)h(M) =H(M) \% f(M)就行了。
等等，f(M)f(M)怎么求？
根据定义，f(x)=det(xI−M)f(x)=det(xI - M)，得到 f(x)=∣xI−M∣=[x−a1−a2−a3⋯−ak−2−ak−1−ak−1x0⋯0000−1x⋯00000−1⋯000⋮⋮⋮⋱⋮⋮⋮000⋯−1x0000⋯0−1x]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−a1)M11+(−a2)M12+⋯+(−ak)M1n=xk−a1xk−1−a2xk−2−⋯−akf(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 Mi,jM_{i,j}代表MM的算术余子式。
只要直接上多项式取模就行了。
最后的总复杂度O(klog⁡2klog⁡2n)\mathcal{O}(k \log_2 k \log_2 n)
我们手动模拟一下多项式取模的过程，其实可以发现我们手动向前压的过程就是在暴力取模。

例题

BZOJ4161 (要求O(k2log⁡2n)\mathcal{O}(k^2\log_2 n))
洛谷P4732 (要求O(klog⁡2klog⁡2n)\mathcal{O}(k\log_2 k\log_2 n))
洛谷P3824 (要求O(k2log⁡2n)\mathcal{O}(k^2\log_2 n))
洛谷P4732的O(klog⁡2klog⁡2n)\mathcal{O}(k\log_2 k\log_2 n)的代码太长了，就不粘在这里了。感兴趣的同学可以来这里看。
附上BZOJ4161的O(k2log⁡2n)\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);
}

评论

共 31 条评论，欢迎与作者交流。

正在加载评论...