LG4723 【模板】常系数线性递推

P4723 【模板】常系数齐次线性递推

题目描述

求一个满足$k$阶齐次线性递推数列${a_i}$的第$n$项。

即：$a_n=\sum\limits_{i=1}^{k}f_i \times a_{n-i}$

输入输出格式

输入格式：

第一行两个数$n$,$k$，如题面所述。

第二行$k$个数，表示$f_1 \ f_2 \ \cdots \ f_k$

第三行$k$个数，表示$a_0 \ a_1 \ \cdots \ a_{k-1}$

输出格式：

一个数，表示 $a_n \% 998244353$ 的值

输入输出样例

输入样例#1： 复制

6 4
3 -1 0 4
-2 3 1 5

输出样例#1： 复制

73

说明

$N = 10^{9} , K = 32000$

题解

先修《数学选修4-2：矩阵与变换》。

常系数线性递推

给出递推式\(f_n=\sum_{i=1}^ka_if_{n-i}\)，和初值条件\(f_1,f_2,\dots ,f_k\)，求\(f_n\)。

可以用生成函数解一下，然后多项式求逆，\(O(n\log n)\)。当然这对于\(n=10^9\)的数据范围是不行的。

矩阵快速幂解法

由递推式，构造矩阵和列向量

\[A=\begin{bmatrix} 0 & 1 & 0 & \dots & 0\\ 0 & 0 & 1 & \dots & 0\\ 0 & 0 & 0 & \dots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \dots & 1\\ a_k & a_{k-1} & a_{k-2} & \dots & a_1 \end{bmatrix} ,F=\begin{bmatrix} f_1\\ f_2\\ f_3\\ \vdots\\ f_k \end{bmatrix} \]

计算\(A^{n-k}F\)即可，\(O(k^3\log n)\)。然而这对于\(k=32000\)的数据范围还是不行。

矩阵的特征值和特征向量

\(A^nF\)是线性变换的形式，自然也可以用特征值与特征向量来求解。

特征多项式是\(f_A(\lambda)=|A-\lambda I|\)，于是特征方程为

\[\det\begin{bmatrix} -\lambda & 1 & 0 & \dots & 0\\ 0 & -\lambda & 1 & \dots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \dots & 1\\ a_k & a_{k-1} & a_{k-2} & \dots & a_1-\lambda \end{bmatrix}=0 \]

手动高斯消元，消成上三角可得

\[(a_1-\lambda+\frac{a_2}{\lambda}+\frac{a_3}{\lambda^2}+\dots+\frac{a_{k}}{\lambda^{k-1}})(-\lambda)^{k-1}=0\\ a_k+a_{k-1}\lambda+a_{k-2}\lambda^2+\dots+a_1\lambda^{k-1}=\lambda^k \]

然后呢，把特征值解出来然后用特征向量那套理论吗？虽然这种方法可行，但是只能做低次的（4次及以下），所以我们需要新科技。

矩阵的多项式

对于\(n\)次多项式\(f(x)\)，将矩阵\(A\)看做自变量带入，得

\[f(A)=a_0E+\sum_{i=1}^na_iA^i \]

记\(f(A)\)为\(A\)的\(n\)次多项式。与另一个\(A\)的\(m\)次多项式\(g(A)\)，其乘法运算满足交换律，即

\[f(A)g(A)=g(A)f(A) \]

Cayley-Hamilton 定理

特征多项式\(f_A(x)=\sum_{i=0}^{k-1}a_{k-i}x^i-x^k\)，这里注意下标。

定理：\(f_A(A)=0\)，即矩阵被自己的特征多项式化零。记忆方法：\(f_A(A)=|A-AI|=0\)。

推论：\(A^n=q(A)f_A(A)+r(A)=r(A)\)，其中\(r(A)=A^n\mod f_A(A)\)。

所以\(A^nF=r(A)F=\sum_{i=0}^{k-1}r_iA^iF\)，\(O(k^4)\)。？？？

多项式优化

将Cayley-Hamilton定理的推论变成普通多项式：\(x^n=q(x)f_A(x)+r(x)\)，其中\(r(x)=x^n\mod f_A(x)\)。

于是我们可以对多项式\(x\)做快速幂并取模\(f_A(x)\)，即可得出\(r(x)\)的系数，即\(r(A)\)的系数，\(O(k\log k\log n)\)。

而\(A^iF=\begin{bmatrix}f_{1+i} & f_{2+i} & \dots & f_{k+i}\end{bmatrix}^T\)，所以\(f_{n+k}=\sum_{i=0}^{k-1}r_if_{k+i}\)。

问题转化成了如何求\(f\)的前\(2k\)项。这时使用生成函数和多项式求逆，\(O(k\log k)\)。

于是我们便得到了\(O(k\log k\log n+k\log k+k)\)的优秀解法。

小优化

刚才说的快速幂是指\(A^{n-k}\)，如果我们做到\(A^{n-1}\)，就只需要前\(k\)项了，那么就用不着什么生成函数了。\(O(k\log k\log n+k)\)。

UPD：然后我发现\(a_i\)系数要翻转，还要取相反数，是我把特征方程写反了？的确数学书上的应该是\(|\lambda I-A|\)，但这没有影响。我看差异是我和其他人转移矩阵列的不一样，我把他们的做了一个中心对称（这是讲义的风格），然后特征方程都不一样了？我觉得这不是我以现有的水平能搞明白的问题。但是大众的写法的优势在于他们的最高次项系数为1，在做暴力多项式取模的时候很好办。

UPD：最近我发现最开始的那个矩阵含有 \(a\) 的最后一行写反了……

void num_trans(polynomial&a,int inverse){
    int limit=a.size(),len=log2(limit);
    static vector<int> bit_rev;
    if(bit_rev.size()!=limit){
        bit_rev.resize(limit);
        for(int i=0;i<limit;++i) bit_rev[i]=bit_rev[i>>1]>>1|(i&1)<<(len-1);
    }
    for(int i=0;i<limit;++i)if(i<bit_rev[i]) swap(a[i],a[bit_rev[i]]);
    for(int step=1;step<limit;step<<=1){
        int gn=fpow(inverse==1?g:g_inv,(mod-1)/(step<<1));
        for(int even=0;even<limit;even+=step<<1){
            int odd=even+step,gk=1;
            for(int k=0;k<step;++k,gk=mul(gk,gn)){
                int t=mul(gk,a[odd+k]);
                a[odd+k]=add(a[even+k],mod-t),a[even+k]=add(a[even+k],t);
            }
        }
    }
    if(inverse==-1){
        int lim_inv=fpow(limit,mod-2);
        for(int i=0;i<limit;++i) a[i]=mul(a[i],lim_inv);
    }
}
polynomial poly_inv(polynomial a,int n){ // mod x^n
    polynomial b(1,fpow(a[0],mod-2));
    if(n==1) return b;
    int limit;
    for(limit=2;limit<n;limit<<=1){
        polynomial a1(a.begin(),a.begin()+limit);
        a1.resize(limit<<1),num_trans(a1,1);
        b.resize(limit<<1),num_trans(b,1);
        for(int i=0;i<limit<<1;++i) b[i]=mul(add(2,mod-mul(a1[i],b[i])),b[i]);
        num_trans(b,-1),b.resize(limit);
    }
    a.resize(limit<<1),num_trans(a,1);
    b.resize(limit<<1),num_trans(b,1);
    for(int i=0;i<limit<<1;++i) b[i]=mul(add(2,mod-mul(a[i],b[i])),b[i]);
    num_trans(b,-1),b.resize(n);
    return b;
}
polynomial poly_div(polynomial f,polynomial g){ // return the quotient
    int n=f.size()-1,m=g.size()-1;
    reverse(g.begin(),g.end()),g.resize(n-m+1),g=poly_inv(g,n-m+1);
    reverse(f.begin(),f.end()),f.resize(n-m+1);
    int limit=1<<int(ceil(log2(2*(n-m)+1)));
    f.resize(limit),g.resize(limit);
    num_trans(f,1),num_trans(g,1);
    for(int i=0;i<limit;++i) f[i]=mul(f[i],g[i]);
    num_trans(f,-1),f.resize(n-m+1);
    return reverse(f.begin(),f.end()),f;
}
polynomial poly_mod(polynomial f,polynomial g){ // return the reminder
    int n=f.size()-1,m=g.size()-1;
    polynomial q=poly_div(f,g);
    int limit=1<<int(ceil(log2(n+1)));
    g.resize(limit),q.resize(limit);
    num_trans(g,1),num_trans(q,1);
    for(int i=0;i<limit;++i) g[i]=mul(g[i],q[i]);
    num_trans(g,-1);
    for(int i=0;i<m;++i) f[i]=add(f[i],mod-g[i]);
    return f.resize(m),f;
}

int n,k;
void mul(polynomial&a,polynomial b,co polynomial&p){
    static co int limit=1<<int(ceil(log2(2*k-1)));
    a.resize(limit),b.resize(limit);
    num_trans(a,1),num_trans(b,1);
    for(int i=0;i<limit;++i) a[i]=mul(a[i],b[i]);
    num_trans(a,-1),a.resize(2*k-1);
    a=poly_mod(a,p);
}
int main(){
    read(n),read(k);
    polynomial a(k),f(k);
    for(int i=1;i<=k;++i) a[k-i]=(mod-read<int>()%mod)%mod; // [-1e9,1e9]
    for(int i=0;i<k;++i) f[i]=(read<int>()%mod+mod)%mod;
    if(n<k) return printf("%d\n",f[n]),0;
    a.push_back(1);
    polynomial rmd(1,1),tmp(2);tmp[1]=1;
    for(;n;n>>=1,mul(tmp,tmp,a))
        if(n&1) mul(rmd,tmp,a);
    int ans=0;
    for(int i=0;i<k;++i) ans=add(ans,mul(rmd[i],f[i]));
    printf("%d\n",ans);
    return 0;
}