常系数齐次线性递推

求一个数列 $f_i$ 的第 $n$ 项，数列满足

f_n=\sum_{i=1}^ka_if_{n-i}

为了这玩应我去把八十中集训那个线性代数课件啃完了qwq

首先矩阵乘法可以做到 $O(k^3logn)$ 。下面考虑对矩阵乘法进行优化。

课件上讲的还是挺好的，~~懒得抄了，直接截图了~~

线性表出的那个结论，归纳法证明即可。

于是我们可以把原来的矩阵表示为一个 $k-1$ 次多项式，表示线性表出这个矩阵的多项式。两个矩阵相乘时暴力展开，得到 $2k-2$ 次多项式，再把多出来的部分用特征多项式‘降次’，我们就可以得到 $O(k^2)$ 的‘矩阵乘法’。

发现这个过程就是多项式乘法和多项式取模的过程，于是可以优化成 $O(klogk)$ 。外面套上快速幂，总体复杂度就是 $O(klogklogn)$ 。

怎么求特征多项式？

那么对于一般的矩阵快速幂，如果它的转移矩阵是固定的，~~可以手玩解方程~~
如果它的转移矩阵是输入的，可以暴力代入点值然后插值，复杂度是 $O(k^4)$ 的（要算k次行列式。。插值就 $k^2$ 插就好了）~~应该不会有这种毒瘤题~~

对于常系数齐次线性递推，可以先写出它的转移矩阵，对第一行进行拉普拉斯展开，可以求得特征多项式为

P(x)=x^k-\sum_{i=1}^ka_ix^{k-i}

最后一个问题，我们求得最后目标矩阵的线性表出多项式之后，怎么还原为值？首先初始列向量为 $f_{k-1}\dots f_1,f_0$ 。左乘 $A^m$ 时的列向量为
$f_{k+m-1}\dots f_{m+1},f_m$ 。发现我们只关注列向量的最后一项，而左乘 $A^0,A^1\dots A^{k-1}$ 时的列向量最后一项恰好为 $f_0,f_1\dots f_{k-1}$ 。故我们根本不用去算矩阵，直接拿多项式各系数去乘对应的初始项就好了。
撒花。~~为什么shadowice的题解那么长啊~~

常数极大近乎要T的代码:（半小时敲完，总共一小时AC，~~说不定考场上真敢写一写？~~）

#include<bits/stdc++.h>
using namespace std;
const int N=2e5+50,mod=998244353;
int n,k,a[N],f[N],r[N],lo[N],rt[N],irt[N],inv[N];unsigned long long P[N];
int Glim(int n){int lim=1;while(lim<n)lim<<=1;return lim;}
void getr(int n){for(int i=1;i<n;i++)r[i]=r[i>>1]>>1^(i&1)<<lo[n];}
int power(int x,int y){
    int z=1;
    for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
    return z;
}
void init(int n){
    int lim=Glim(n);
    for(int i=2,j=0;i<=lim;i<<=1,j++)inv[i]=power(i,mod-2),lo[i]=j;
    int w=power(3,(mod-1)/lim);rt[lim>>1]=1;
    for(int i=(lim>>1)+1;i<lim;i++)rt[i]=1ll*rt[i-1]*w%mod;
    for(int i=(lim>>1)-1;i;i--)rt[i]=rt[i<<1];
    for(int i=1;i<lim;i++)irt[i]=power(rt[i],mod-2);
}
void NTT(int *a,int n,int op){
    for(int i=0;i<n;i++)P[i]=a[r[i]];
    for(int i=1,*w=op>0?rt:irt;i<n;i<<=1)
        for(int t=i<<1,j=0;j<n;j+=t)
            for(int k=0;k<i;k++){
                int y=P[i|j|k]*w[i|k]%mod;
                P[i|j|k]=P[j|k]+mod-y;P[j|k]+=y;
            }
    for(int i=0;i<n;i++)a[i]=P[i]%mod;
    if(op<0)for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv[n]%mod;
}
int c[N],d[N];
void Inv(int *a,int *b,int n){
    if(n==1){b[0]=power(a[0],mod-2);return;}
    int m=n+1>>1,lim=Glim(n+m);Inv(a,b,m);
    for(int i=0;i<lim;i++)c[i]=i<n?a[i]:0,d[i]=i<m?b[i]:0;
    getr(lim);NTT(c,lim,1);NTT(d,lim,1);
    for(int i=0;i<lim;i++)d[i]=1ll*d[i]*(2+mod-1ll*c[i]*d[i]%mod)%mod;
    NTT(d,lim,-1);for(int i=m;i<n;i++)b[i]=d[i];
}
int ra[N],rb[N],irb[N],e[N],g[N];
void Div(int *a,int *b,int *d,int n,int m){
    for(int i=0;i<=n;i++)ra[i]=a[i];reverse(ra,ra+n+1);
    for(int i=0;i<=m;i++)rb[i]=b[i];reverse(rb,rb+m+1);
    for(int i=m+1;i<=n-m;i++)rb[i]=0;Inv(rb,irb,n-m+1);
    int lim=Glim(2*(n-m+1));for(int i=n-m+1;i<lim;i++)ra[i]=irb[i]=0;
    getr(lim);NTT(ra,lim,1);NTT(irb,lim,1);
    for(int i=0;i<lim;i++)c[i]=1ll*ra[i]*irb[i]%mod;
    NTT(c,lim,-1);reverse(c,c+n-m+1);
    for(int i=n-m+1;i<max(lim,m);i++)c[i]=0;
    lim=Glim(2*m);getr(lim);
    for(int i=0;i<lim;i++)e[i]=i<m?b[i]:0,g[i]=i<m?c[i]:0;
    NTT(e,lim,1);NTT(g,lim,1);
    for(int i=0;i<lim;i++)d[i]=1ll*e[i]*g[i]%mod;
    NTT(d,lim,-1);for(int i=0;i<lim;i++)d[i]=i<m?(a[i]+mod-d[i])%mod:0;
}
int Z[N],X[N],ff[N],ans;
int main(){
    scanf("%d%d",&n,&k);a[k]=1;init(k*4);
    for(int i=1;i<=k;i++)scanf("%d",&a[k-i]);
    for(int i=0;i<k;i++)(a[i]=-a[i]%mod)<0?a[i]+=mod:0;
    for(int i=0;i<k;i++)scanf("%d",&f[i]),f[i]%=mod;
    Z[0]=1;X[1]=1;int lim=Glim(2*k);
    for(;n;n>>=1){
        if(n&1){
            for(int i=0;i<lim;i++)ff[i]=X[i];
            getr(lim);NTT(ff,lim,1);NTT(Z,lim,1);
            for(int i=0;i<lim;i++)Z[i]=1ll*Z[i]*ff[i]%mod;
            NTT(Z,lim,-1);Div(Z,a,ff,2*k-2,k);
            for(int i=0;i<lim;i++)Z[i]=i<k?ff[i]:0;
        }
        getr(lim);NTT(X,lim,1);
        for(int i=0;i<lim;i++)X[i]=1ll*X[i]*X[i]%mod;
        NTT(X,lim,-1);Div(X,a,ff,2*k-2,k);
        for(int i=0;i<lim;i++)X[i]=i<k?ff[i]:0;
    }
    for(int i=0;i<k;i++)ans=(ans+1ll*Z[i]*f[i])%mod;
    printf("%d\n",(ans+mod)%mod);
    return 0;
}

附 $O(k^2logn)$ 代码（写MTT是不可能的，这辈子都不可能写MTT的）

#include<bits/stdc++.h>
using namespace std;
const int N=2020,mod=1e9+7;
int n,k,a[N],h[N],ans;
struct node{
    int c[N*2];
    void clear(){memset(c,0,sizeof(c));}
    void init(){memset(c,0,sizeof(c));c[0]=1;}
    node operator *(const node &b){
        node z;z.clear();
        for(int i=0;i<k;i++)
            for(int j=0;j<k;j++)
                z.c[i+j]=(z.c[i+j]+1ll*c[i]*b.c[j])%mod;
        for(int i=2*k-2;i>=k;z.c[i--]=0)
            for(int j=1;j<=k;j++)
                z.c[i-j]=(z.c[i-j]+1ll*a[j]*z.c[i])%mod;
        return z;
    }
}b;
node ksm(node x,int y){
    node z;z.init();
    for(;y;y>>=1,x=x*x)if(y&1)z=z*x;
    return z;
}
int main(){
    scanf("%d%d",&n,&k);
    for(int i=1;i<=k;i++)scanf("%d",&a[i]);
    for(int i=0;i<k;i++)scanf("%d",&h[i]);
    b.c[1]=1;b=ksm(b,n);
    for(int i=0;i<k;i++)ans=(ans+1ll*h[i]*b.c[i])%mod;
    printf("%d\n",(ans+mod)%mod);
    return 0;
}

还有一些东西：由哈密顿卡莱定理，发现如果一个递推关系是连乘一个相同的矩阵，那么这个矩阵的任意一个位置都满足同一个线性递推关系（即矩阵的特征多项式），也可以推出矩阵的任意一些位置的和也满足相同的线性递推。（这样如果一些dp要求某一行或某一列之和之类的就可以求了）
然后还要注意：线性递推关系不一定就需要快速幂之类的：如果一个dp是类似矩阵乘法的样子，那么可以确定它满足线性递推。那么如果要求 $\leq n$ 的所有值，可以求出初值和递推式后用分治FFT解决（线性递推其实就是分治FFT的形式）而初值和递推式怎么求就要看题了。最不济就是暴力dp求初值，用BM算法求递推关系。而如果这个矩阵很有特点，可能可以dp求递推式（考虑行列式求特征多项式）（而往往还能考虑转移路径进一步用组合数优化）然后也可以快速幂求一个很靠后的值。

Prutekoi's blog

常系数齐次线性递推

一些需要记住的结论