多项式全家桶 - 常系数齐次线性递推
常系数齐次线性递推
注意：这篇是一个不用很高线性代数基础的理解方法
（图解党狂喜）
问题描述
\(f_0,f_1...f_{k-1}\) 已知，并已知递推系数 \(a_1,a_2...a_k\)
\(f_n=\sum\limits_{i=1}^{k} f_{n-k}a_k\)
求 \(f_n\)
\(n\le 1145141919810,k\le 114514\)
常规的矩阵快速幂过不去了。我们需要一个带 \(\log\) 的做法。
换一个思路
考虑把转移看成一个带权 DAG 图：点 \(u\) 和点 \(u-i\) 有一个有向边，边权为 \(a_i\)。（\(1\le i\le k\)，后面如果没提到 \(i\) 的取值范围就是这个范围）
然后我们来考虑点 \(n\) 的状态把点 \([0,k-1]\) 算了多少遍。假设我们能求出来“多少遍”，那直接就可以求出 \(f_n\) 了：每个 \(f_{0...k-1}\)，乘以对应的遍数，加起来，就行了。
然后我们考虑用拓扑排序求这个东西。设 \(c\) 表示遍数数组。
我们的边界条件是 \(n\)：一遍，\(c_n=1\)
对于每个当前入度为 \(0\) 的点 \(u\)，遍历它的出边 \(v\)，设边权为 \(w\)。那么 c[v]+=c[u]*w，显然。
然后这么做一遍就可以得到 \(c_{0...k-1}\)，对应的乘起来即可。复杂度 \(O(nk)\)
然后你会问：艹这nm不就是直接递推么
是的
但是我们反过来考虑了个寂寞 ... 吗？
魔法
考虑 \(c\) 的生成函数，设为 \(C(x)=\sum\limits_{i=0}^{\infin} x^ic_i\)。
一开始我们知道这玩意就是 \(x^n\) （在 \(n\) 位置有一个 \(n\)，其余为 \(0\)）
我们令拓扑排序删边的时候，如果一个点的出度为 \(0\) （即出边被删没了），且不是 \([0,k-1]\) 中的点，就把这个点也删除，同时标记它的 \(c\) 为 \(0\) （这样显然没问题，因为这个点已经不用了）
然后这样显然只会剩下 \(c_{0..k-1}\) 有可能有值（其它都是 \(0\)）
考虑一次拓扑排序做了甚么：假设当前入度为 \(0\) 的点是 \(u\)，那相当于：
- 标记 \(c_{u-i}\) 加上 \(c_u\times a_i\)
- 然后标记 \(c_u\) 为 \(0\)，即 \(c_u\) 减去 \(c_u\)
一共操作了 \(k+1\) 个数，并且操作数都是 \(c_u\) 的倍数。
在生成函数上考虑，它相当于减去了 \(c_u\) 倍的：
\(x^u-\sum\limits_{i=1}^{k} x^{u-i}a_i\)
这个东西的系数是 \(\{1,-a_1,-a_2...-a_k\}\)，然后后面都是 \(0\)。把后面的 \(0\) 去掉，只保留这 \(k+1\) 项系数，构成的多项式被称作这个递推式的 特征多项式。设它为 \(F\)
注：如果您有了解过斐波那契的通项是怎么推的，那您应该接触过这个“特征多项式”的定义
斐波那契的特征多项式为：\(x^2-x-1\)
然后每次我们就是把 \(F\) 和当前 \(C\) 的最高位（设为 \(u\)）对齐，然后减去 \(F\) 的 \(c_u\) 倍 —— 这是拓扑排序干的事情。
我们发现它的本质就是在求余数
回想一下小学竖式计算除法的过程，又因为 \(F\) 的最高位为 \(1\)，所以减去 \(F\) 的 \(c_u\) 倍其实就是在做带余除法求余数
所以，\(C=x^n\mod F\)
然后我们只要求出 \(x^n\mod F\) 的值，就可以知道遍数数组 \(c_{0..k-1}\)，对应的乘以初始值 \(f_{0..k-1}\)，就得到了 \(f_n\)
然后这个 \(C\) 怎么算呢？\(n\) 老大了
考虑多项式快速幂 —— 不用 \(\ln/\exp\)，而是普通的倍增快速幂，参考整数的写法，并把“模”的操作改成多项式取模就可以了。
#include <bits/stdc++.h>
using namespace std;
namespace Flandre_Scarlet
{
#define N 2000006
#define GG 3
#define GI 332748118
#define mod 998244353
#define int long long
#define F(i,l,r) for(int i=l;i<=r;++i)
#define D(i,r,l) for(int i=r;i>=l;--i)
#define Fs(i,l,r,c) for(int i=l;i<=r;c)
#define Ds(i,r,l,c) for(int i=r;i>=l;c)
#define MEM(x,a) memset(x,a,sizeof(x))
#define FK(x) MEM(x,0)
#define Tra(i,u) for(int i=G.st(u),v=G.to(i);~i;i=G.nx(i),v=G.to(i))
#define p_b push_back
#define sz(a) ((int)a.size())
#define all(a) a.begin(),a.end()
#define iter(a,p) (a.begin()+p)
#define PUT(a,n) F(i,1,n) printf("%d ",a[i]); puts("");
int I() {char c=getchar(); int x=0; int f=1; while(c<'0' or c>'9') f=(c=='-')?-1:1,c=getchar(); while(c>='0' and c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar(); return ((f==1)?x:-x);}
template <typename T> void Rd(T& arg){arg=I();}
template <typename T,typename...Types> void Rd(T& arg,Types&...args){arg=I(); Rd(args...);}
void RA(int *p,int n) {F(i,1,n) *p=I(),++p;}
int qpow(int a,int b,int m=mod) {int r=1; while(b) {if (b&1) r=r*a%m; a=a*a%m,b>>=1;} return r;}
int pgg[30],pgi[30];
void init()
{
F(i,0,23)
{
pgg[i]=qpow(GG,(mod-1)>>i);
pgi[i]=qpow(GI,(mod-1)>>i);
}
}
#define ad(x,y) ((x)+(y)>mod?(x)+(y)-mod:(x)+(y))
#define dc(x,y) ((x)-(y)<0?(x)-(y)+mod:(x)-(y))
namespace poly
{
int w[N],r[N];
void NTT(int f[N],int lim,int type)
{
F(i,0,lim-1) if (i<r[i]) swap(f[i],f[r[i]]);
for(int mid=1,pp=1;mid<lim;mid<<=1,++pp)
{
int Wn=(type==1?pgg[pp]:pgi[pp]);
w[0]=1; F(i,1,mid-1) w[i]=w[i-1]*Wn%mod;
for(int i=0;i<lim;i+=(mid<<1))
{
for(int j=0;j<mid;++j)
{
register int y=f[i|mid|j]*w[j]%mod;
f[i|mid|j]=dc(f[i|j],y);
f[i|j]=ad(f[i|j],y);
}
}
}
if (type==-1)
{
int ivv=qpow(lim,mod-2);
F(i,0,lim-1) f[i]=f[i]*ivv%mod;
}
}
#define REV F(i,0,lim-1) r[i]=((r[i>>1]>>1)|((i&1)?len:0));
int pool[10][N];
inline void Inv(int f[N],int g[N],int n) // 求逆, pool 0~2
{
int *iv=pool[0],*a=pool[1],*b=pool[2];
F(i,0,4*n) iv[i]=a[i]=b[i]=0;
iv[0]=qpow(f[0],mod-2);
int len=1,lim=1;
for(len=1;len<=(n<<1);len<<=1)
{
lim=len<<1;
F(i,0,4*len) a[i]=b[i]=0;
F(i,0,len-1) a[i]=f[i],b[i]=iv[i];
REV;
NTT(a,lim,1); NTT(b,lim,1);
F(i,0,lim-1) a[i]=(2*b[i]-a[i]*b[i]%mod*b[i]%mod+mod)%mod;
NTT(a,lim,-1);
F(i,0,len-1) iv[i]=a[i];
}
F(i,0,n-1) g[i]=iv[i];
F(i,n,lim-1) g[i]=0;
}
inline void mul(int f[N],int g[N],int n,bool flag=1) // *=, pool 3~4
// flag: 是否对n取膜
{
int len=1,lim=1; while(len<=(n<<1)) len=lim,lim<<=1;
int *a=pool[3],*b=pool[4];
F(i,0,lim-1) a[i]=b[i]=0;
F(i,0,n-1) a[i]=f[i],b[i]=g[i];
REV;
NTT(a,lim,1); NTT(b,lim,1);
F(i,0,lim-1) a[i]=a[i]*b[i]%mod;
NTT(a,lim,-1);
F(i,0,2*n-1) f[i]=a[i]; F(i,2*n,lim) f[i]=0;
if (flag) F(i,n,2*n-1) f[i]=0;
}
void Mod(int f1[N],int f2[N],int n,int m,int R[N]) // 多项式取模, pool 5~6
// 这里只保留了余数, 商开到 pool 里而不是传参数修改了
{
int *a=pool[5],*b=pool[6],*Q=pool[7];
F(i,0,4*n) a[i]=b[i]=0;
F(i,0,n-1) a[i]=f1[n-1-i]; F(i,n-m+1,n-1) a[i]=0; // a=f1_r%(x^(n-m+1))
F(i,0,m-1) b[i]=f2[m-1-i]; F(i,n-m+1,m-1) b[i]=0;
Inv(b,b,n-m+1);
mul(a,b,n-m+1);
F(i,0,n-m) Q[i]=a[n-m-i];
F(i,0,4*n) a[i]=b[i]=0;
F(i,0,n-1) a[i]=f1[i];
F(i,0,m-1) b[i]=f2[i];
int len=1,lim=1; while(len<=n) len=lim,lim<<=1;
REV;
NTT(b,lim,1); NTT(Q,lim,1);
F(i,0,lim-1) b[i]=b[i]*Q[i]%mod;
NTT(b,lim,-1); NTT(Q,lim,-1);
F(i,0,m-2) R[i]=(a[i]-b[i]%mod+mod)%mod; F(i,m-1,lim-1) R[i]=0;
}
void PowMod(int f[N],int p,int g[N],int n,int h[N]) // f^p%g, g有n项, 保存在h
{
int *res=pool[8];
F(i,0,8*n) res[i]=0;
res[0]=1; // res=1
while(p)
{
if (p&1)
{
mul(res,f,n,0); // res*=f
Mod(res,g,2*n,n,res); // res%=g
}
mul(f,f,n,0); // f=f*f
Mod(f,g,2*n,n,f); // f%=g
p>>=1;
}
F(i,0,n-2) h[i]=res[i]; F(i,n,8*n) h[i]=0;
}
}
int n,k;
int a[N],t[N];
void Input()
{
Rd(n,k);
F(i,1,k) t[i]=(I()%mod+mod)%mod;
F(i,0,k-1) a[i]=(I()%mod+mod)%mod;
}
int f[N],g[N],c[N];
void Sakuya()
{
init();
f[1]=1; // 这玩意就是 x
F(i,1,k) g[k-i]=(mod-t[i]); g[k]=1; // 特征多项式
// f^n%g -> c
poly::PowMod(f,n,g,k+1,c);
int ans=0;
F(i,0,k-1) ans+=a[i]*c[i]%mod;
ans%=mod;
printf("%lld\n",ans);
}
void IsMyWife()
{
Input();
Sakuya();
}
}
#undef int //long long
int main()
{
Flandre_Scarlet::IsMyWife();
getchar();
return 0;
}
