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

· · 题解

本问题的大致做法基本上是全网统一的，都是快速幂求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. 代码

稍微有点压行。

#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';
}