题解 P4723 【【模板】常系数齐次线性递推】
BJpers2
2020-06-30 21:06:15
本问题的大致做法基本上是全网统一的，都是快速幂求$x^n$对特征多项式$p$取模的结果。
翻阅网上的题解，我发现它们对这种做法的解释都用到了线性代数内容，包括矩阵特征值，行列式，C-H定理等。
这些内容对于有基础的选手来说并不至于非常艰深，然而不可否认的是，一个事先并没有了解过这些线性代数内容的选手，面对一篇充满全新概念的题解，很可能会产生理解的障碍。
事实上，我们并不需要使用任何线性代数的知识来导出这个做法。
---
### 1. 返璞归真
让我们忘掉矩阵，行列式和特征值，采用小学层面的方法求斐波那契数列（一个经典的例子）的第$n$项。
我们换一个角度，倒过来求。假设我们现在欲求$fib_5$。
我们可以列出如下过程：$fib_5=fib_4+fib_3=2fib_3+fib_2=3fib_2+2fib_1=5fib_1+3fib_0$
最后带入$f_0,f_1$计算即可。
这当中的推导过程不是什么秘密，不过是每次取最前面那一项拿递推式展开罢了。
### 2. 本质
考虑一下上述过程本质上究竟干了什么，你可以先停一分钟自己好好思考一下。如果你想到了，本题解对你来说就已经没有阅读的价值。
> 上述过程，就是通过每次消去最高次项，求$x^5$对斐波那契数列的特征多项式——$x^2-x-1$取模的结果的过程。
如果你还没有反应过来，我可以展示一下：
$$
\begin{aligned}
& 0x^0+0x^1+0x^2+0x^3+0x^4+1x^5 \\
\xrightarrow{-1x^3(x^2-x-1)} & 0x^0+0x^1+0x^2+1x^3+1x^4 \\
\xrightarrow{-1x^2(x^2-x-1)} & 0x^0+0x^1+1x^2+2x^3 \\
\xrightarrow{-2x^1(x^2-x-1)} & 0x^0+2x^1+3x^2 \\
\xrightarrow{-3x^0(x^2-x-1)} & 3x^0+5x^1
\end{aligned}
$$
两者之间的逻辑关联很简单，因为前者是每次取和式里下标最大那一项，替换成两个下标较小的项之后加回去；而后者是每次从多项式里减去若干倍的$x^2$，并替换成相应倍数的$x+1$之后加回去。
### 3. 推广
理解了上述逻辑之后，我们能轻易将其推广到一般线性其次递推式的情形：
$$
f_i=p_1f_{i-1}+p_2f_{i-2}+\dots+p_kf_{i-k}
$$
要求$f_n$，只要求多项式$x^n$对
$$
p(x)=x^k-p_1x^{k-1}-p_2x^{k-2}-\dots-p_kx^0
$$
取模的结果。因为「多项式取模」的过程恰好与「每次按递推式展开下标最大的那项」的过程一致，都可以用「替换」来解读。
### 4. 多项式取模
接下来基本上就是多项式取模的板子，会多项式取模的可以直接跳过。
我们写一个快速幂，求多项式$x$的$n$次方，只不过每次乘法都对上面提到的特征多项式$p(x)$取模。
那么需要解决的问题变为快速求一个$l=2k-2$次多项式对一个$k$次多项式取模的结果。
也即需把一个$l$次多项式$a$展开成$a=pq+r$的形式，其中$q$是商，$r$是余数（也就是我们的目标）。
我们设$a^R$为多项式$a$的反转（第$l$项和第$0$项互换，第$l-1$项和第$1$项互换，以此类推），易得$a^R\equiv p^Rq^R \bmod x^{l-k+1}$。
因此我们只需求出$p^R$在模$x^{l-k+1}$意义下的逆元，乘上$a^R$即可得到$q^R$。最后利用$r=a-pq$即可求出我们需要的$r$。
### 5. 后记
本片题解介绍的思路应该有很多前人想到，但是令我比较吃惊的是没有人把它发表出来（或者是发表了我没有看到）。希望这篇题解能让初次接触线性递推的选手不要被纷繁复杂的数学公式劝退，少走弯路，更深入地理解这个算法。
### 6. 代码
稍微有点压行。
```cpp
#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';
}
```