线性递推式的计算

线性代数

目录

• 线性递推式的计算
  • 特征多项式
    • 特征向量与特征值
    • 算法过程
    • 一个我一直不知道的东西
    • 练习

特征多项式

对于常系数齐次递推式，虽然我们可以用矩阵快速幂求了，但是如果递推式的阶数$k$比较大，那就很麻烦了，我们这里介绍一种神奇的优化。

特征向量与特征值

• 对于一个$n \times n$的矩阵$M$，我们可以理解为一个对于$n$维空间的线性变换；
• 注意到，变换前后，有些向量的张成空间是不变的，即这个变换对于这些向量的作用效果只是拉伸；
• 可能因为这些向量对于这个变换算是一种特征吧，所以我们记这些向量（$\vec{0}$是平凡的，所以通常会忽略）为$M$的特征向量，其拉伸程度$\lambda$为对应的特征值；
• 写成式子就是：$M\vec{v}=\lambda\vec{v}$，意思就是变换$M$对于向量$\vec{v}$的作用和一个标量$\lambda$的作用是相同的，都是拉伸$\lambda$倍。
• 现在考虑求这个特征向量，移项得$M\vec{v}-\lambda\vec{v}=\vec{0}$。
• 可是矩阵和标量之间没有定义减法啊，不过这个$\lambda$从作用上讲也是对于$\vec{v}$的线性变换，所以式子可以写为$(M-\lambda I)\vec{v}=\vec{0}$。
• 这个式子有非零解当且仅当$\det(M-\lambda I)=0$，因为从对空间线性变换的角度考虑就是一个非零向量$\vec{v}$经过$M-\lambda I$变换后变为$\vec{0}$，这个显然只有变换$M-\lambda I$将空间压缩至低维时才有的，此时根据行列式的几何意义有$\det(M-\lambda I)=0$。
• 现在看来，重要的是特征值$\lambda$了，因为它直接决定$\vec{v}$是否存在。
• 将那个行列式展开，我们可以得到一个关于$\lambda$的$n$次多项式$f(\lambda)$，我们称它为矩阵$M$的特征多项式（实际上为了方便，特征多项式通常指$\det(\lambda I -M)$，因为这样就没有一堆$(-1)^x$项了），方程$f(\lambda)=0$称为特征方程。
• 然后有个凯莱-哈密顿定理说$f(M)=0$，就是把多项式$f(\lambda)$看作关于矩阵的多项式把$M$带进去会出来全零矩阵。这个我不太会证，不过感性理解一下就好了吧。

算法过程

• 如果我们写出了递推式的转移矩阵$M$，那么按照矩阵快速幂的方法，求数列第$n$项实际上就是求$M^x$；
• 根据最后取列向量的哪个元素，$x$可以取$[n-k+1,n]$，通常为了减小常数都是取列向量第一个元素，所以就是求$M^{n-k+1}$，但是下面的做法取最后一个元素更方便，于是就是求$M^n$；
• 想办法用一下刚刚说的定理，我们如果求出$x^n=f(x)g(x)+h(x)$，$f(x)$是$M$的特征向量，$h(x)$是$k-1$次多项式，那么将矩阵$M$带入，我们发现只有$h(M)$有用；
• 这就很好了，我们只要求出$M$的特征多项式$f(x)$，快速幂求出$x^n\bmod{f(x)}$就可以了，这个暴力求是$O(k^2\log{n})$的，用FFT可以做到$O(k\log{k}\log{n})$。
• 然后先考虑一下怎么带入矩阵$M$，直接带可是$O(k^4)$，显然不行；
• 不过我们要求的只是一个元素，设$g(x)=x^n\bmod{f(x)}$，那么就是求列向量$\vec{v}=(a_{k-1},a_{k-2},\cdots,a_0)$左乘$g(M)$的第$k$个元素，所以式子可以这样化简：
  $$\sum_{i=0}^{k-1}[x^i]g(x)\times(\vec{v}M^{i})_k=\sum_{i=0}^{k-1}[x^i]g(x)\times a_i$$
• 所以求出$\{a_i\}$的第$0$到$k-1$项就可以了。通常题目是给你的，如果没有，基本上可以用生成函数化为多项式求逆的形式，用FFT可以做到$O(k\log{k})$。
• 现在说一下怎么得到特征多项式，先观察一下：
  $$\begin{vmatrix} \lambda-b_1 & b_2 & \cdots & b_{k-1} & b_k \\ 1 & \lambda & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & \lambda & 0 \\ 0 & 0 & \cdots & 1 & \lambda \\ \end{vmatrix}$$
• 我们发现大部分地方都是$0$，不妨对第一行拉普拉斯展开：
  $$\det(M)=\sum_{i=1}^k(-1)^{i+1}M_{1,i}$$
• 其中每个余子式$M_{1,i}$是很有规律的，举个例子写一下可以发现，$f(\lambda)=\lambda^k-\sum_{i=1}^k a_{k-i}\lambda^i$。
• 于是，求$k$阶常系数线性齐次递推式的第$n$项，我们得到了一个朴素做法$O(k^2\log{n})$，FFT可优化到$O(k\log{k}\log{n})$的优秀做法。

模板题：[BZOJ 4161] Shlw loves matrixI

朴素版：

#include <cstdio>
#include <cstring>
#define mset(a,x,n) memset(a,x,(n)*sizeof(*(a)))
#define mcpy(a,b,n) memcpy(a,b,(n)*sizeof(*(a)))
const int MAXN=2005,P=1e9+7;
int a[MAXN],b[MAXN];
int f[MAXN],g[MAXN],h[MAXN];
inline void inc(int &x,int y){x+=y,x= x>=P ? x-P:x;}
inline int pow(int x,int y){
    int ret=1;
    for(;y;x=(long long)x*x%P,y>>=1){
        if(y&1) ret=(long long)ret*x%P;
    }
    return ret;
}
inline int inv(int x){return pow(x,P-2);}
inline void poly_mul(int a[],int b[],int p[],int m){
    static int c[MAXN<<1];
    mset(c,0,(m<<1)-1);
    for(int i=0;i<m;++i){
        for(int j=0;j<m;++j)
            inc(c[i+j],(long long)a[i]*b[j]%P);
    }
    for(int i=(m<<1)-2;i>=m;--i){
        int k=(long long)c[i]*inv(p[m])%P;
        for(int j=0;j<=m;++j)
            inc(c[i-j],P-(long long)k*p[m-j]%P);
    }
    mcpy(a,c,m);
}
inline void poly_pow(int x[],int y,int p[],int m,int ret[]){
    ret[0]=1;
    for(;y;poly_mul(x,x,p,m),y>>=1){
        if(y&1) poly_mul(ret,x,p,m);
    }
}
int main(){
    int n,k;
    scanf("%d%d",&n,&k);
    for(int i=1;i<=k;++i){
        scanf("%d",&a[i]);
        a[i]=(a[i]+P)%P;
    }
    for(int i=0;i<k;++i){
        scanf("%d",&b[i]);
        b[i]=(b[i]+P)%P;
    }
    if(n<k){
        printf("%d",b[n]);
        return 0;
    }
    for(int i=1;i<=k;++i) f[k-i]=P-a[i];
    f[k]=1;
    g[1]=1;
    poly_pow(g,n,f,k,h);
    int ans=0;
    for(int i=0;i<k;++i)
        inc(ans,(long long)h[i]*b[i]%P);
    printf("%d",ans);
    return 0;
}

一个我一直不知道的东西

• 对于$f_n=\sum_{i=1}^n f_{n-i}g_i$，$g$为$n$项的数列（这里认为其大于$n$的项为$0$），这个显然是个线性递推式，不过我们可以用生成函数来搞事：
• 我们有：
  $$\begin{eqnarray*} F(x) &=& f_0+\sum_{i=1}^{+\infty}x^i\sum_{j=1}^if_{i-j}g_j \\ &=& f_0+\sum_{i=1}^{+\infty}\sum_{j=1}^if_{i-j}x^{i-j}g_jx^j \\ &=& f_0+(F \ast G)(x) \end{eqnarray*}$$
• 所以$F(x)=\frac{f_0}{1-G(x)}$。
• 那么要是求形如$[x^n]\frac{a}{1-f(x)}$，其实就是求一个线性递推式的第$n$项，不要只想到多项式求逆！有时只要求某一项的系数，但是$n$特别大，显然不能求逆啊...可是我以前只知道想求逆...
• 拓展一下可以得知：形如$\frac{a}{f(x)+g(x)}$其中$g(x)$的次数是有限项的式子都可以变为一个线性递推式（就是上面的推导过程中$f_0$变为一个有限项的多项式，也就是数列的前几项不符合卷积形式）。

练习

[NOI 2017] 泳池

• 吉老师的神仙题Orz...
• 首先就是将恰好为$K$的转化为$\leq K$的减去$<K$的。
• 然后考虑DP，设$h(i)$表示高为$1001$、宽为$i$的矩形区域的合法概率，那么$h(n)$即为答案。
• 容易发现每一个安全子矩形都是从海滩向上延伸一定高度后被一个不安全格子卡住，然后分成几部分更窄的矩形继续向上延伸。所以高层的可以看作低层的子问题，我们从高到低DP（感觉从下往上也可以，但是可能不太好写。说实话这题的思路我现在理解的还是不顺畅，不能用一种很连贯的方式讲出来）。
• 那么设$f(i,j)$表示一个宽为$i$、顶部为$1001$，下底边为$j+1$的矩形的合法概率。
• 转移的话其实不难理解，就是枚举一个最低点卡住，两边是更高的矩形：$f(i,j)=\sum_{k=1}^j \sum_{a \geq i}\sum_{b>i}f(a,k-1)q^{(a-i)(j-1)}(1-q)f(b-i,j-k)q^{(b-i)(j-k)}$；
• 有一点要注意：就是由于最低点可能是一段，为了不重复，我们强制分开的两个更高的矩形一边是严格高于最低点，一边是不低于。
• 考虑用前缀和优化，就是设$g_j(i)=\sum_{k>j}f(i,k)q^{i(k-j)}$，每次$g_j(i)$转移到$g_{j+1}(i)$时，加上$f(i,j)$再乘上$q^i$就可以了。那么每层转移是$O(宽度^2)$的。
• 由于下底边为$j+1$的矩形要是对我们的答案有贡献，则其下方必须全为安全区域，所以它的宽不超过$\lfloor\frac{K}{j}\rfloor$，所以总复杂度是：
  $$O\left(\sum_{i=1}^K \left(\lfloor\frac{K}{j}\rfloor\right)^2\right)=O(K^2\zeta(2))=O(\frac{\pi^2}{6}K^2)=O(K^2)$$
• 容易得到$h(i)=g_0(i)+\sum_{j=1}^K h(i-j)(1-q)g_0(j-1)$，就是首先是宽为$i$的海滩都是安全的，然后是枚举一个不安全的断点、右边再拼一个安全的海滩。
• 对于$i>K+1$的$h(i)$，显然没有合法的宽为$i$的安全海滩，即$g_0(i)=0$，那么就是一个常系数线性递推式的形式。
• 对于本题，暴力多项式取模即可，总复杂度$O(K^2\log{n})$。
• 实际上吉老师本来的$K=20000$，复杂度要求$O(K\log{K}\log{n})$，上面的DP部分也要用多项式求逆优化到$O(K\log^2{K})$，就先不研究了Orz

代码：

#include <cstdio>
#include <algorithm>
#include <cstring>
#define mset(a,x,n) memset(a,x,(n)*sizeof(*(a)))
#define mcpy(a,b,n) memcpy(a,b,(n)*sizeof(*(a)))
using std::min;
const int MAXN=1005,P=998244353;
int n,q,K;
int powq[MAXN];
int f[MAXN][MAXN],g[MAXN],h[MAXN<<1];
int a[MAXN],b[MAXN];
inline void inc(int &x,int y){x+=y,x= x>=P ? x-P:x;}
inline int pow(int x,int y){
    int ret=1;
    for(;y;x=(long long)x*x%P,y>>=1){
        if(y&1) ret=(long long)ret*x%P;
    }
    return ret;
}
inline int inv(int x){return pow(x,P-2);}
inline void DP(int K){
    mset(g,0,K+1);
    g[0]=1;
    for(int i=K;i;--i){
        int m=min(n,K/i);
        mset(f[i],0,m+1);
        for(int j=1;j<=m;++j){
            for(int k=1;k<=j;++k){
                inc(f[i][j],(long long)g[k-1]*g[j-k]%P);
                if(k!=K) inc(f[i][j],(long long)g[k-1]*f[i][j-k]%P);
            }
            f[i][j]=(long long)f[i][j]*(1-q+P)%P;
        }
        for(int j=1;j<=m;++j)
            g[j]=(long long)(g[j]+f[i][j])%P*powq[j]%P;
    }
}
inline void poly_mul(int a[],int b[],int ret[],int p[],int m){
    static int c[MAXN<<1];
    mset(c,0,(m<<1)-1);
    for(int i=0;i<m;++i){
        for(int j=0;j<m;++j)
            inc(c[i+j],(long long)a[i]*b[j]%P);
    }
    for(int i=(m<<1)-1;i>=m;--i){
        int k=(long long)c[i]*p[m]%P;
        for(int j=0;j<=m;++j)
            inc(c[i-j],P-(long long)k*p[m-j]%P);
    }
    mcpy(ret,c,m+1);
}
inline void poly_pow(int ret[],int y,int p[],int m){
    static int x[MAXN];
    mset(x,0,m+1),mset(ret,0,m);
    x[1]=1,ret[0]=1;
    for(;y;y>>=1){
        if(y&1) poly_mul(ret,x,ret,p,m);
        poly_mul(x,x,x,p,m);
    }
}
inline int cal(int K){
    if(K==0) return pow((1-q+P)%P,n);
    DP(K);
    mcpy(h,g,K+1);
    for(int i=K+1;i;--i)
        g[i]=(long long)(1-q+P)*g[i-1]%P;
    for(int i=1;i<=K;++i){
        for(int j=1;j<=i;++j)
            inc(h[i],(long long)h[i-j]*g[j]%P);
    }
    if(n<=K) return h[n];
    a[K+1]=1;
    for(int i=K+1;i;--i)
        a[K+1-i]=P-g[i];
    poly_pow(b,n,a,K+1);
    int ret=0;
    for(int i=0;i<=K;++i)
        inc(ret,(long long)b[i]*h[i]%P);
    return ret;
}
inline void pre_cal(){
    powq[0]=1;
    for(int i=1;i<=K;++i)
        powq[i]=(long long)powq[i-1]*q%P;
}
int main(){
    int x,y;
    scanf("%d%d%d%d",&n,&K,&x,&y);
    q=(long long)x*inv(y)%P;
    pre_cal();
    printf("%d",(cal(K)-cal(K-1)+P)%P);
    return 0;
}