简介

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

设有数列${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} \times \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$

$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$

证毕。

优化递推

设转移矩阵为$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)$代码

#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) // notice 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);
    #ifdef LOCAL
        system("pause");
    #endif
}