题解 P4723 【【模板】线性递推】

2018-07-05 21:20:12


终于调过了...

本题中的线性递推,特指常系数齐次递推。

常系数齐次线性递推(之后简称线性递推)指的是形如:

$$a(x)=\sum_{i=1}^ka(x-i)f(i)$$

的递推式,其中$a(0)\to a(k-1)$都是已知数。

可以知道,这种递推可以通过构造矩阵来实现快速递推。

对于线性递推,转移矩阵为:

$$A=\begin{pmatrix}a_1&a_2&a_3&\dots&a_{k-2}&a_{k-1}&a_k\\1&0&0&\dots&0&0&0\\0&1&0&\cdots&0&0&0\\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\0&0&0&\cdots&1&0&0\\0&0&0&\cdots&0&1&0\end{pmatrix}$$

这个矩阵乘上列向量$F=\begin{pmatrix}f_{k-1}\\f_{k-2}\\\vdots\\f_0\end{pmatrix}$可以得到$\begin{pmatrix}f_k\\f_{k-1}\\\vdots\\f_1\end{pmatrix}$。

那么,答案变成:

$$\sum_{i=0}^{k-1}(FA^n)_i$$

利用矩阵快速幂,可以得到$O(k^3logn)$的算法。对于本题不可承受。

我们尝试转换思想。矩阵乘法的时间复杂度为$O(k^3)$,而多项式乘法为$O(k^2)$(借助$FFT$可以优化至$O(klogk)$),明显更加优秀,所以考虑将矩阵乘法转换为多项式乘法。

可以发现,$a(x)=\displaystyle\sum_{i=0}^{k-1}(FA^x)_i$这个式子可以通过转移变成$\displaystyle\sum_{i=0}^{k-1}p_ia_i$($p_i$为系数)。那么,我们考虑通过什么操作可以将其转换。

(此处需要特征多项式的相关知识)

定义矩阵$A$的特征多项式为

$$g(\lambda)=det(\lambda E-A)$$

其中$E$为单位矩阵。

那么,根据Cayley-Hamilton 定理 可知,$g(A)=0$。

那么,我们要求的变成了:

$\displaystyle\sum_{i=0}^{k-1}(FA^n)_i=\sum_{i=0}^{k-1}[F[A^n\bmod g(A)]]_i$

此时存在两种做法,一种是直接计算,时间复杂度$O(k^2logn)$,另一种是套用多项式取余,做到$O(klogklogn)$。我们明显是需要后者。

接下来的问题就只剩下如何计算$g(\lambda)$了。

$g(\lambda)=det(\lambda E-A)$

$=det\begin{pmatrix}\lambda-a_1&-a_2&-a_3&\dots&-a_{k-2}&-a_{k-1}&-a_k\\-1&\lambda&0&\dots&0&0&0\\0&-1&\lambda&\cdots&0&0&0\\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\0&0&0&\cdots&-1&\lambda&0\\0&0&0&\cdots&0&-1&\lambda\end{pmatrix}$

考虑对第一行展开分别求和,得:

$g(\lambda)=(\lambda-a_1)A_{1,1}+(-a_2)A_{1,2}+\cdots+(-a_k)A_{1,k}$

$=\lambda^k-a_1\lambda^{k-1}-\cdots-a_k$

这样就可以在$O(k)$的时间复杂度求出特征多项式了。

代码:(大常数警告)

#include<bits/stdc++.h>
#define For(i,a,b) for(i=(a);i<=(b);++i)
#define Forward(i,a,b) for(i=(a);i>=(b);--i)
#define Rep(i,a,b) for(register int i=(a),i##end=(b);i<=i##end;++i)
#define Repe(i,a,b) for(register int i=(a),i##end=(b);i>=i##end;--i)
#define Chkmin(a,b) a=a<b?a:b
using namespace std;
template<typename T>inline void read(T &x){
    T s=0,f=1;char k=getchar();
    while(!isdigit(k)&&k^'-')k=getchar();
    if(!isdigit(k)){f=-1;k=getchar();}
    while(isdigit(k)){s=s*10+(k^48);k=getchar();}
    x=s*f;
}
void file(void){
    freopen("polynomial.in","r",stdin);
    freopen("polynomial.out","w",stdout);
}
const int MAXN=1<<20;

typedef long long ll;

namespace polynomial
{
    static int mod=998244353,gen=3,g[21],rev[MAXN],Len;

    inline int ad(int a,int b){return (a+=b)>=mod?a-mod:a;}

    inline int power(int a,int b)
    {
        static int sum;
        for(sum=1;b;b>>=1,a=(ll)a*a%mod)if(b&1)
            sum=(ll)sum*a%mod;
        return sum;
    }

    inline void predone()
    {
        static int i,j;
        for(i=1,j=2;i<=19;++i,j<<=1)g[i]=power(gen,(mod-1)/j);
    }

    inline void calrev(int Len)
    {
        static int Logl;Logl=(int)floor(log(Len)/log(2)+0.3)-1;
        Rep(i,1,Len-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<Logl);
    }

    inline void NTT(int X[],int typ)
    {
        Rep(i,1,Len-1)if(i<rev[i])swap(X[i],X[rev[i]]);
        static int i,j,k,kk,w,t,wn,r;
        for(k=2,kk=1,r=1;k<=Len;k<<=1,kk<<=1,++r)
        {
            wn=g[r];
            for(i=0;i<Len;i+=k)for(j=0,w=1;j<kk;++j,w=(ll)w*wn%mod)
            {
                t=(ll)w*X[i+j+kk]%mod;
                X[i+j+kk]=ad(X[i+j],mod-t);
                X[i+j]=ad(X[i+j],t);
            }
        }
        if(typ==-1)
        {
            reverse(X+1,X+Len);
            static int invn;invn=power(Len,mod-2);
            Rep(i,0,Len-1)X[i]=(ll)X[i]*invn%mod;
        }
    }

    static int x[MAXN],y[MAXN];
    inline void mul(int a[],int b[])
    {
        memset(x,0,sizeof x);memset(y,0,sizeof y);
        Rep(i,0,(Len>>1)-1)x[i]=a[i],y[i]=b[i];
        NTT(x,1);NTT(y,1);
        Rep(i,0,Len-1)x[i]=(ll)x[i]*y[i]%mod;
        NTT(x,-1);
        Rep(i,0,Len-1)a[i]=x[i];
    }

    static int A[MAXN],B[MAXN];

    void Inv(int *a,int *b,int n)
    {
        if(n==1){b[0]=power(a[0],mod-2);return;}
        Inv(a,b,n>>1);
        Len=n<<1;
        calrev(Len);
        Rep(i,0,(Len>>1)-1)A[i]=a[i],B[i]=b[i];
        NTT(A,1);NTT(B,1);
        Rep(i,0,Len-1)B[i]=(ll)B[i]*B[i]%mod*A[i]%mod;
        NTT(B,-1);
        Rep(i,0,(Len>>1)-1)b[i]=ad(b[i],ad(b[i],mod-B[i]));
        Rep(i,0,Len)A[i]=B[i]=0;
    }

    static int X[MAXN],Y[MAXN],TT[MAXN];

    inline void Div(int *a,int n,int *b,int m)
    {
        if(n<m){Rep(i,0,m-1)b[i]=a[i];return;}
        memcpy(X,a,sizeof X);memcpy(Y,b,sizeof Y);
        reverse(b,b+m+1);reverse(X,X+n+1);
        Rep(i,n-m+1,n)X[i]=0;
        memset(TT,0,sizeof TT);
        for(Len=2;Len<=(n-m+1);Len<<=1);
        Inv(b,TT,Len);
        memcpy(b,TT,sizeof TT);
        while(Len<=(n<<2))Len<<=1;
        calrev(Len);
        mul(X,b);reverse(X,X+n-m+1);
        Rep(i,n-m+1,n)X[i]=0;
        mul(Y,X);
        Rep(i,0,m-1)b[i]=ad(a[i],mod-Y[i]);
        memcpy(a,X,sizeof X);
    }
}
using namespace polynomial;

static int n,F[MAXN],G[MAXN],m,AA[MAXN],KK[MAXN];

static int R[MAXN];

inline void MMM(int *A,int *B,int *Mo,int lp)
{
    Chkmin(lp,m<<1);
    calrev(Len);mul(A,B);
    memcpy(R,Mo,sizeof R);
    Div(A,lp,R,m);
    Rep(i,m,m<<1)A[i]=0;
    Rep(i,0,m-1)A[i]=R[i];
}

inline void KSM(int *A,int x,int *B,int *Mo)
{
    for(int i=1;x;x>>=1,MMM(A,A,Mo,2<<i),++i)
    {
        if(x&1)MMM(B,A,Mo,m<<1);
        if((1<<(i-1))>=(m<<1))--i;
    }
}

int main(void){
//    file();
    predone();
    read(n);read(m);
    for(Len=2;Len<=(m<<1);Len<<=1);
    Rep(i,1,m)read(F[m-i]),F[m-i]=ad(0,mod-F[m-i]);
    F[m]=1;
    Rep(i,0,m-1)read(AA[i]),AA[i]=ad(AA[i],mod);
    G[1]=1;KK[0]=1;
    KSM(G,n,KK,F);
    static int ans=0;
    Rep(i,0,m-1)ans=ad(ans,(ll)AA[i]*KK[i]%mod);
    cout<<ans<<endl;
//  cerr<<1.0*clock()/CLOCKS_PER_SEC<<endl;
    return 0;
}