从分式第n项到线性递推
Potassium
2021-08-11 21:25:29
介绍一个常数小还好写的科技：[bostan-mori 算法](https://arxiv.org/pdf/2008.08822.pdf)。
[博客园](https://www.cnblogs.com/Potassium/p/15130342.html) [github 博客](https://potassiumwings.github.io/2021/08/11/from_n_th_element_in_fraction_to_linear_recurrence_array/)
不需要任何矩阵知识，前置知识：多项式乘法
应用场景：**需要注意这一个算法要求模数是质数。**
# 波斯坦-茉莉算法简介
这个东西是求分式第 $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}$
$n=0$ 时，显有 $[x^0]\frac fg=\frac{f_0}{g_0}$。
然后就可以 $\mathcal{O}(k\log k\log n)$ 求出解。
示例代码如下：
```c++
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)` 为与本题相同的含义）：
```c++
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]$ 项。
参考实现：
```c++
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` 套用波斯坦-茉莉算法即可，非常的方便。
# 代码
```c++
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;
}
```
[完整代码](https://www.luogu.com.cn/record/55601927)就是封装了一些多项式乘除法，就不贴了，不会多项式的大常数选手写的很垃圾