【洛谷P4723】【模板】—线性递推（多项式取模）

传送门

---

考虑说要求一个
$a_n=sum_{i=1}^{k}a_{n-i}f_i$

写成矩阵的形式就是

$A_n=A_{n-1}F=A_0F^n$

实际上$F^n$最后就是一个长度为$k$的向量
满足$a_n=sum_{i=1}^{k}c_ia_i$
这样的形式

由$Cayley-Hamilton$定理可以得到$F$的特征多项式

$g(λ)=det(λI-F)=λ^k-f_{1}λ^{k-1}-f_{2}λ^{k-2}……-f_k=0$

而且$g$最高次只有$k-1$次

所以我们可以倍增多项式取模$O(klog_klog_n)$求出$F^n\% g$
由于$g(λ)=0$，所以$F^n\% g=F^n$

然后就完了

注意读入的数最小有$-1e9$

#include<bits/stdc++.h>
using namespace std;
#define gc getchar
inline int read(){
	char ch=gc();
	int res=0,f=1;
	while(!isdigit(ch))f^=ch=='-',ch=gc();
	while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
	return f?res:-res;
}
#define re register
#define pb push_back
#define cs const
#define pii pair<int,int>
#define fi first
#define se second
#define ll long long
#define poly vector<int>
#define bg begin
#define int long long
cs int mod=998244353,G=3;
inline int add(int a,int b){return (a+=b)>=mod?a-mod:a;}
inline void Add(int &a,int b){(a+=b)>=mod?(a-=mod):0;}
inline int dec(int a,int b){return (a-=b)<0?a+mod:a;}
inline void Dec(int &a,int b){(a-=b)<0?(a+=mod):0;}
inline int mul(int a,int b){return 1ll*a*b>=mod?1ll*a*b%mod:a*b;}
inline void Mul(int &a,int b){a=mul(a,b);}
inline int ksm(int a,int b,int res=1){
	for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));return res;
}
inline void chemx(ll &a,ll b){a<b?a=b:0;}
inline void chemn(int &a,int b){a>b?a=b:0;}
cs int N=128005,C=19;
poly w[C+1];
inline void init_w(){
	for(int i=1;i<=C;i++)w[i].resize(1<<(i-1));
	int wn=ksm(G,(mod-1)/(1<<C));
	w[C][0]=1;
	for(int i=1;i<(1<<(C-1));i++)w[C][i]=mul(w[C][i-1],wn);
	for(int i=C-1;i;i--)
	for(int j=0;j<(1<<(i-1));j++)
	w[i][j]=w[i+1][j<<1];
}
int rev[N<<2];
inline void init_rev(int lim){
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)*(lim>>1));
}
inline void ntt(poly &f,int lim,int kd){
	for(int i=0;i<lim;i++)if(i>rev[i])swap(f[i],f[rev[i]]);
	for(int mid=1,l=1;mid<lim;mid<<=1,l++)
	for(int i=0;i<lim;i+=(mid<<1))
		for(int j=0,a0,a1;j<mid;j++)
			a0=f[i+j],a1=mul(f[i+j+mid],w[l][j]),
			f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
	if(kd==-1&&(reverse(f.bg()+1,f.bg()+lim),1))
		for(int inv=ksm(lim,mod-2),i=0;i<lim;i++)Mul(f[i],inv);
}
inline poly operator +(poly a,poly b){
	poly c;int lim=max(a.size(),b.size());c.resize(lim);
	a.resize(lim),b.resize(lim);
	for(int i=0;i<lim;i++)c[i]=add(a[i],b[i]);return c;
}
inline poly operator -(poly a,poly b){
	poly c;int lim=max(a.size(),b.size());c.resize(lim);
	a.resize(lim),b.resize(lim);
	for(int i=0;i<lim;i++)c[i]=dec(a[i],b[i]);return c;
}
inline poly operator *(poly a,int b){
	for(int i=0;i<a.size();i++)Mul(a[i],b);return a;
}
inline poly operator /(poly a,int b){
	for(int i=0,inv=ksm(b,mod-2);i<a.size();i++)Mul(a[i],inv);
	return a;
}
inline poly operator *(poly a,poly b){
	int deg=a.size()+b.size()-1,lim=1;
	if(deg<=128){
		poly c(deg,0);
		for(int i=0;i<a.size();i++)
		for(int j=0;j<b.size();j++)
			Add(c[i+j],mul(a[i],b[j]));
		return c;
	}
	while(lim<deg)lim<<=1;init_rev(lim);
	a.resize(lim),ntt(a,lim,1);
	b.resize(lim),ntt(b,lim,1);
	for(int i=0;i<lim;i++)Mul(a[i],b[i]);
	ntt(a,lim,-1),a.resize(deg);
	return a;
}
inline poly Inv(poly a,int deg){
	poly c,b(1,ksm(a[0],mod-2));
	for(int lim=4;lim<(deg<<2);lim<<=1){
		init_rev(lim);
		c=a,c.resize(lim>>1);
		c.resize(lim),ntt(c,lim,1);
		b.resize(lim),ntt(b,lim,1);
		for(int i=0;i<lim;i++)Mul(b[i],dec(2,mul(b[i],c[i])));
		ntt(b,lim,-1),b.resize(lim>>1);
	}b.resize(deg);return b;
}
inline poly operator /(poly a,poly b){
	int lim=1,deg=(int)a.size()-(int)b.size()+1;
	reverse(a.bg(),a.end());
	reverse(b.bg(),b.end());
	while(lim<=deg)lim<<=1;
	b=Inv(b,lim),b.resize(deg);
	a=a*b,a.resize(deg);
	reverse(a.bg(),a.end());
	return a;
}
inline poly operator %(poly a,poly b){
	int deg=(int)a.size()-(int)b.size()+1;
	if(deg<0)return a;
	poly c=a-(a/b)*b;
	c.resize(b.size()-1);
	return c;
}
int n,m,a[N];
poly f,g;
signed main(){
	init_w();
	n=read(),m=read();
	f.resize(m+1);
	for(int i=1;i<=m;i++)f[m-i]=read(),f[m-i]=mod-(f[m-i]%mod+mod)%mod;
	f[m]=1;
	for(int i=0;i<m;i++)a[i]=read(),a[i]=(a[i]%mod+mod)%mod;
	poly res,g;
	res.resize(m+1),res[0]=1;
	g.resize(m+1),g[1]=1;
	for(;n;n>>=1,g=g*g%f){if(n&1)res=res*g%f;}
	int anc=0;
	for(int i=0;i<res.size();i++)Add(anc,mul(res[i],a[i]));
	cout<<anc;
}