题解 P4723 【【模板】线性递推】
Great_Influence
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 定理](https://en.wikipedia.org/wiki/Cayley%E2%80%93Hamilton_theorem) 可知，$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)$,另一种是套用[多项式取余](https://www.luogu.org/problemnew/show/P4512)，做到$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)$的时间复杂度求出特征多项式了。
代码:(大常数警告)
```cpp
#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;
}
```