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

1. 返璞归真

让我们忘掉矩阵，行列式和特征值，采用小学层面的方法求斐波那契数列（一个经典的例子）的第nn项。

我们换一个角度，倒过来求。假设我们现在欲求fib5fib5​。

我们可以列出如下过程：fib5=fib4+fib3=2fib3+fib2=3fib2+2fib1=5fib1+3fib0fib5​=fib4​+fib3​=2fib3​+fib2​=3fib2​+2fib1​=5fib1​+3fib0​

最后带入f0,f1f0​,f1​计算即可。

这当中的推导过程不是什么秘密，不过是每次取最前面那一项拿递推式展开罢了。

2. 本质

考虑一下上述过程本质上究竟干了什么，你可以先停一分钟自己好好思考一下。如果你想到了，本题解对你来说就已经没有阅读的价值。

上述过程，就是通过每次消去最高次项，求x5x5对斐波那契数列的特征多项式——x2−x−1x2−x−1取模的结果的过程。

如果你还没有反应过来，我可以展示一下：

0x0+0x1+0x2+0x3+0x4+1x5→−1x3(x2−x−1)0x0+0x1+0x2+1x3+1x4→−1x2(x2−x−1)0x0+0x1+1x2+2x3→−2x1(x2−x−1)0x0+2x1+3x2→−3x0(x2−x−1)3x0+5x1−1x3(x2−x−1)​−1x2(x2−x−1)​−2x1(x2−x−1)​−3x0(x2−x−1)​​0x0+0x1+0x2+0x3+0x4+1x50x0+0x1+0x2+1x3+1x40x0+0x1+1x2+2x30x0+2x1+3x23x0+5x1​

两者之间的逻辑关联很简单，因为前者是每次取和式里下标最大那一项，替换成两个下标较小的项之后加回去；而后者是每次从多项式里减去若干倍的x2x2，并替换成相应倍数的x+1x+1之后加回去。

3. 推广

理解了上述逻辑之后，我们能轻易将其推广到一般线性其次递推式的情形：

fi=p1fi−1+p2fi−2+⋯+pkfi−kfi​=p1​fi−1​+p2​fi−2​+⋯+pk​fi−k​

要求fnfn​，只要求多项式xnxn对

p(x)=xk−p1xk−1−p2xk−2−⋯−pkx0p(x)=xk−p1​xk−1−p2​xk−2−⋯−pk​x0

取模的结果。因为「多项式取模」的过程恰好与「每次按递推式展开下标最大的那项」的过程一致，都可以用「替换」来解读。

4. 多项式取模

接下来基本上就是多项式取模的板子，会多项式取模的可以直接跳过。

我们写一个快速幂，求多项式xx的nn次方，只不过每次乘法都对上面提到的特征多项式p(x)p(x)取模。

那么需要解决的问题变为快速求一个l=2k−2l=2k−2次多项式对一个kk次多项式取模的结果。

也即需把一个ll次多项式aa展开成a=pq+ra=pq+r的形式，其中qq是商，rr是余数（也就是我们的目标）。

我们设aRaR为多项式aa的反转（第ll项和第00项互换，第l−1l−1项和第11项互换，以此类推），易得aR≡pRqR mod xl−k+1aR≡pRqRmodxl−k+1。

因此我们只需求出pRpR在模xl−k+1xl−k+1意义下的逆元，乘上aRaR即可得到qRqR。最后利用r=a−pqr=a−pq即可求出我们需要的rr。

5. 后记

本片题解介绍的思路应该有很多前人想到，但是令我比较吃惊的是没有人把它发表出来（或者是发表了我没有看到）。希望这篇题解能让初次接触线性递推的选手不要被纷繁复杂的数学公式劝退，少走弯路，更深入地理解这个算法。

6. 代码

稍微有点压行。

#include<iostream>
#include<cstdio>
#include<algorithm>
#define FF FOF(i,0,L)
#define FOR(i,a,b) for(int i=a;i<=b;i++)
#define FOF(i,a,b) for(int i=a;i< b;i++)
#define ROF(i,a,b) for(int i=a;i>=b;i--)
#define CON(a,b) DFT(a,L);FF a[i]=1ll*a[i]*b[i]%P;IFT(a,L);
using namespace std;
const int N=150150,P=998244353;
int n,m,K,L,o,A,Q;
int rv[N],W[N],p[N],f[N],t[N],a[N],b[N];
void ipt(int&x){scanf("%d",&x);x%=P;x+=x>>31&P;}
int qpw(int x,int y){int z=1;for(;y;y>>=1,x=1ll*x*x%P)if(y&1) z=1ll*z*x%P;return z;}
void ini(int n){
	for(L=1;L<=n;L<<=1,o++);
	FF rv[i]=rv[i>>1]>>1|(i&1)<<(o-1);
	int V=qpw(3,P>>o);W[L>>1]=1;
	FOF(i,(L>>1)+1,L) W[i]=1ll*W[i-1]*V%P;
	ROF(i,(L>>1)-1,1) W[i]=W[i<<1];
}
void DFT(int*A,int L){
	static unsigned long long B[N];
	int u=o-__builtin_ctz(L),t;
	FF B[i]=A[rv[i]>>u];
	for(int i=1;i<L;i<<=1)for(int j=0,s=i<<1;j<L;j+=s)FOF(k,0,i)
		t=B[i+j+k]*W[i+k]%P,B[i+j+k]=B[j+k]+P-t,B[j+k]+=t;
	FF A[i]=B[i]%P;
}
void IFT(int*A,int L){
	reverse(A+1,A+L);DFT(A,L);
	int V=P-(P-1)/L;
	FF A[i]=1ll*A[i]*V%P;
}
int upl(int n){return 1<<(32-__builtin_clz(n));}
void INV(int*A,int*B,int n){
	static int C[N],L;
	if(!n) return B[0]=qpw(A[0],P-2),void();
	INV(A,B,n>>1);L=upl(n<<1);
	FF C[i]=i>n?0:A[i];
	DFT(C,L);DFT(B,L);
	FF B[i]=(2-1ll*C[i]*B[i]%P+P)*B[i]%P;IFT(B,L);
	FF B[i]=i>n?0:B[i],C[i]=0;
}
void MUL(int*A,int*B){
	static int C[N];
	FF C[i]=A[i];DFT(C,L);
	CON(B,C);
	FF C[i]=i>K?0:B[n-i];
	CON(C,t);
	FF C[i]=i>K?0:C[i];
	reverse(C,C+K+1);
	CON(C,p);
	FF (B[i]+=P-C[i])%=P;
}
int main(){
	scanf("%d%d",&Q,&m);
	ini(m<<1);n=m-1<<1;K=n-m;
	FOR(i,1,m) ipt(p[i]);
	FOF(i,0,m) ipt(f[i]);
	p[0]=P-1;
	INV(p,t,K);DFT(t,L);
	reverse(p,p+m+1);DFT(p,L);
	a[1]=b[0]=1;
	for(;Q;Q>>=1,MUL(a,a))if(Q&1) MUL(a,b);
	FOF(i,0,m) (A+=1ll*b[i]*f[i]%P)%=P;
	cout<<A<<'\n';
}