本文是在洛谷博客、github 博客同時發布的P4723 【模板】常系數齊次線性遞推題解。
介紹一個常數小還極其好寫的科技：bostan-mori 算法。
與其他線性遞推算法不同，利用本方法求解線性遞推問題不需要任何矩陣知識，唯一前置知識：多項式乘法
波斯坦-茉莉算法簡介
這個東西是求分式第 \(n\) 項，即 \([x^n]\frac {f(x)}{g(x)}\) 的，而我們知道分式第 n 項和線性遞推式可以很容易地互化的（后文會細說）。於是我們先看如何利用波斯坦-茉莉算法求解分式第 \(n\) 項。
\(\begin{aligned}&[x^n]\frac{f(x)}{g(x)}\\=&[x^n]\frac{f(x)g(-x)}{g(x)g(-x)}\\=&[x^n]\frac{c_{even}(x^2)+x\cdot c_{odd}(x^2)}{v(x^2)}\\=&[x^{n/2}]\frac{c_{even}(x)}{v(x)},N\text{ is even}\\&[x^{n/2}]\frac{c_{odd}(x)}{v(x)},N\text{ is odd} \end{aligned}\)
然后就可以 \(\mathcal{O}(k\log k\log n)\) 求出解。
示例代碼如下：
int divAt(Poly F,Poly G, ll k){
int i;
for(;k;k>>=1){
Poly R=G;
// R=G(-x)
for(i=1;i<R.size();i+=2)R[i]=mod-R[i];
F*=R,G*=R;
for(i=k&1;i<F.size();i+=2)F[i/2]=F[i];
F.resize(i/2);
for(i=0;i<G.size();i+=2)G[i/2]=G[i];
G.resize(i/2);
}
return F.empty()?0:F[0]*qpow(G[0],mod-2)%mod;
}
應用場景：需要注意這一個算法要求模數是質數。
好，回到線性遞推。
從線性遞推到分式第 n 項
如果你會兩者之間的轉換，可以直接跳到第三節。
我們求分式第 \(n\) 項可以化作線性遞推后求解：
設 \(\frac{f(x)}{g(x)}=h(x)\)，其中 \(deg(f)=m,deg(g)=k\)，則 \(f=g\ast h\)。
根據多項式乘法，項數 \(i\) 很大的時候 \(f_i=0\)， \(h\) 就是一個遞推系數為 \(-\frac{g_{1\cdots k}}{g_0}\) 的線性遞推式：
\(0=f_i=g_0h_i+g_1h_{i-1}+\cdots+g_{k}h_{i-k}\)
\(h_i=\sum_{j=1}^{k}\frac{-g_j}{g_0}h_{i-j}\)
而 \(h\) 的前 \(0\cdots m\) 項可以容易地通過 \(f\ast g^{-1}(\bmod x^k)\) 得出，於是可以求解。
\(m>k\) 的時候可以用 \(\frac fg\) 余數進行計算；\(m=k\) 的時候求出的是 \(h_0\cdots h_k\)，移一位即可獲得答案。
參考實現如下（其中 get_an(F,A,n,k) 為與本題相同的含義）：
int div_at(Poly F,Poly G,ll n){
int m=F.size(),k=G.size();
if(m>k){
Poly P=F/G,R=F-P*G;
return div_at(R,G,n)+P.at(n);
}
Poly h=poly.inv(G,k)*F;
if(m==k){
for(int i=0;i<k&&i+1<h.size();i++)h[i]=h[i+1];
n--;
}
h.cut(k-1); // 就是 h %= x^{k-1}
int invg0=qpow(G[0],mod-2);
G*=-invg0;
return get_an(G,h,m,n-1);
}
從分式第 n 項到線性遞推
我們有一個以 \(F_{1\cdots k}\) 為遞推系數的線性遞推序列 \(h\)，想要構造 \(f,g\) 使得 \([x^n]\frac fg=h_n\)
首先，根據多項式乘法，構造項數 \(i\) 特別大的時候（\(f_i=0\)）的遞推關系。套用上面的式子：
\(h_ig_0=\sum_{j=1}^{k}{-g_j}h_{i-j}\)
令 \(g_0=0,g_i=F_i(i\in[1,k])\)即可。
又：\(f=gh(\bmod x^k)\)
所以可以求出 \(f\) 的 \([0,k-1]\) 項。
參考實現：
int getAn(Poly F,Poly A,ll n,int k){
F=Poly{1}-F;
Poly f=A*F;
f.cut(k); // 就是 h %= x^k
return divAt(f,F,n);
}
其中 divAt 套用波斯坦-茉莉算法即可，非常的方便。
代碼
namespace POLY{
int divAt(Poly F,Poly G, ll k){
int i;
for(;k;k>>=1){
Poly R=G;
// R=G(-x)
for(i=1;i<R.size();i+=2)R[i]=mod-R[i];
F*=R,G*=R;
for(i=k&1;i<F.size();i+=2)F[i/2]=F[i];
F.resize(i/2);
for(i=0;i<G.size();i+=2)G[i/2]=G[i];
G.resize(i/2);
}
return F.empty()?0:F[0]*qpow(G[0],mod-2)%mod;
}
int getAn(Poly F,Poly A,ll n,int k){
F=Poly{1}-F;
Poly f=A*F;f.cut(k);
return divAt(f,F,n);
}
}
int main(){
int i,x,n,k;
scanf("%d%d",&n,&k);
Poly F(k+1),A(k);
F[0]=0;
for(i=1;i<=k;i++){
scanf("%d",&x);
F[i]=(x+mod)%mod;
}
for(i=0;i<k;i++){
scanf("%d",&x);
A[i]=(x+mod)%mod;
}
printf("%d\n",POLY::getAn(F,A,n,k));
return 0;
}
完整代碼就是封裝了一些多項式乘除法，就不貼了，不會多項式的大常數選手寫的很垃圾。