多项式全家桶学习笔记
nikangle
·
·
算法·理论
本文的大部分运算在模意义下进行。若无特殊说明，所有多项式最高次项可以为零。
文中可能有一些不严谨的地方，可以感性理解一下，然后当结论来记。
多项式求逆
对于 n-1 次多项式 f，求 n-1 次多项式 g，使得 f\times g\equiv 1\pmod {x^n}。
注意到 g 本质上是无穷次多项式 f^{-1} 对 x^m 取模的结果，所以定义 g_k 表示 f^{-1} 对 x^k 取模的结果。只要求出任意 g_k\left(k\ge n\right) 即可。
考虑倍增。首先，显然有 g_1=\left(f_0\right)^{-1}。接下来只需要求出对于任意 k，在已知 g_k 的情况下，求出 g_{2k}。
注意到 g_{2k}-g_k\equiv 0\pmod{x^k}，并且若 a\equiv 0\pmod{x^p}，则 a^2\equiv 0\pmod{x^{2p}}，所以考虑将同余式左边平方，并进一步化简：
\begin{aligned}
g_{2k}^2-2g_kg_{2k}+g_k^2&\equiv 0&\pmod{x^{2k}}\\
f\left(g_{2k}^2-2g_kg_{2k}+g_k^2\right)&\equiv 0&\pmod{x^{2k}}\\
g_{2k}-2g_k+fg_k^2&\equiv 0&\pmod{x^{2k}}\\
g_{2k}&\equiv 2g_k-fg_k^2&\pmod{x^{2k}}\\
g_{2k}&\equiv g_k\left(2-fg_k\right)&\pmod{x^{2k}}
\end{aligned}
那么，我们就得到了 g_k 的递推式，单次递推需要两次多项式乘法，时间复杂度为 \mathcal O\left(k\log_2k\right)。单次求逆时间复杂度为 \mathcal O\left(n\log_2n\right)，但是常数比较大。
Code
void getinv(long long*f,long long*g,int n){//n 为 2 的幂次
long long*h=new long long[(n<<1)+5];
for(int i=0;i<(n<<1);i++) g[i]=h[i]=0;
g[0]=pow(f[0],mod-2);
for(int i=2;i<=n;){
for(int j=0;j<i;j++) h[j]=f[j];//此处 i=2k
i<<=1,ntt(g,i,1),ntt(h,i,1);//此处 i=4k
for(int j=0;j<i;j++) h[j]=h[j]*g[j]%mod;
ntt(h,i,-1);
h[0]=mo(2-h[0]+mod);
for(int j=1;j<i;j++) h[j]=mo(mod-h[j]);
ntt(h,i,1);
for(int j=0;j<i;j++) g[j]=g[j]*h[j]%mod;
ntt(g,i,-1);
for(int j=i>>1;j<i;j++) g[j]=0;
}
delete[]h;
}
多项式开根
对于 n-1 次多项式 f，求 n-1 次多项式 g，使得 g^2\equiv f\pmod{x^n}。
类似于多项式求逆，定义 g_k 表示 \sqrt f 对 x^k 取模的结果。只要求出任意 g_k\left(k\ge n\right) 即可。
仍然考虑倍增。g_1=\sqrt {f_0}，使用 Cipolla 二次剩余求解即可。注意二次剩余可能有多个解，以任意一个解作为初始值 g_1 都可以求出一个不同的 g。
接下来，考虑由 g_k 推到 g_{2k}。用类似的方法：
\begin{aligned}
g_{2k}-g_k&\equiv 0&\pmod{x^k}\\
g_{2k}^2-2g_kg_{2k}+g_k^2&\equiv 0&\pmod{x^{2k}}\\
f-2g_kg_{2k}+g_k^2&\equiv 0&\pmod{x^{2k}}\\
2g_kg_{2k}&\equiv f+g_k^2&\pmod{x^{2k}}\\
g_{2k}&\equiv\dfrac{f}{2g_k}+\dfrac{g_k}2&\pmod{x^{2k}}
\end{aligned}
单次递推需要一次多项式求逆，时间复杂度为 \mathcal O\left(k\log_2k\right)。单次开根时间复杂度为 \mathcal O\left(n\log_2n\right)，常数很大。
Code
long long val;
struct node{
long long a,b;
node operator*(const node&x)const{
return node{(a*x.a+b*x.b%mod*val)%mod,(a*x.b+b*x.a)%mod};
}
};
long long cipolla(long long x){
long long tmp;
mt19937 rnd(time(NULL));
do{
tmp=rnd()%mod;
}while(pow((tmp*tmp-x+mod)%mod,mod>>1)==1);
val=(tmp*tmp-x+mod)%mod;
return pow(node{tmp,1},(mod+1)>>1).a;
}
void getsqrt(long long*f,long long*g,int n){//n 为 2 的幂次
long long*h=new long long[(n<<1)+5],h2=new long long[(n<<1)+5],inv=pow(2,mod-2);
for(int i=0;i<n;i++) h[i]=h2[i]=g[i]=0;
g[0]=cipolla(f[0]);//此处 g[0] 取 mod-g[0] 也可以，cipolla 的返回值为任意一个
for(int i=2;i<=n;){
getinv(g,h,i);
for(int j=0;j<i;j++) h2[j]=f[j];
i<<=1,ntt(h,i,1),ntt(h2,i,1);
for(int j=0;j<i;j++) h[j]=h[j]*h2[j]%mod;
ntt(h,i,-1);
for(int j=0;j<i;j++) g[j]=j<(i>>1)?inv*(g[j]+h[j])%mod:0;
}
delete[]h;
}
多项式除法
对于 n 次多项式 f 和 m\left(m\le n\right) 次多项式 g，求 n-m 次多项式 q 和 m-1 次多项式 r，满足 f=g\times q+r。
令 h\left(x\right) 表示将 x 带入多项式 h 得到的新多项式。对于 k 次多项式 h，rev\left(h\right)=x^kh\left(\dfrac 1x\right)，即将 h 的系数向量反过来。
对原式进行进行一个导的推：
\begin{aligned}
f\left(x\right)&=g\left(x\right)\times q\left(x\right)+r\left(x\right)\\
f\left(\dfrac 1x\right)&=g\left(\dfrac 1x\right)\times q\left(\frac 1x\right)+r\left(\frac 1x\right)\\
x^nf\left(\dfrac 1x\right)&=x^mg\left(\dfrac 1x\right)\times x^{n-m}q\left(\dfrac 1x\right)+x^nr\left(\dfrac 1x\right)\\
rev\left(f\right)&=rev\left(g\right)\times rev\left(q\right)+x^{n-m+1}rev\left(r\right)\\
rev\left(f\right)&\equiv rev\left(g\right)\times rev\left(q\right)&\pmod{x^{n-m+1}}\\
rev\left(q\right)&\equiv\dfrac{rev\left(f\right)}{rev\left(g\right)}&\pmod{x^{n-m+1}}
\end{aligned}
发现 q 正好是 n-m 次多项式，所以，我们就可以求出 q 了。然后把 q 代入 f=g\times q+r，就能求出 r。
时间复杂度 \mathcal O\left(n\log_2n\right)，常数很大。
Code
void division(long long*f,int n,long long*g,int m,long long*q,long long*r){//g 的最高次项必为 0
int k=max(4<<__lg(n-m),2<<__lg(n));
long long *h=new long long[k+5],*h2=new long long[k+5];
k=2<<__lg(n-m);
reverse(g,g+1+m),getinv(g,h,k),reverse(g,g+1+m);
for(int i=0;i<k;i++) h2[i]=f[n-i];
for(int i=k;i<(k<<1);i++) h2[i]=0;
k<<=1,ntt(h,k,1),ntt(h2,k,1);
for(int i=0;i<k;i++) h[i]=h[i]*h2[i]%mod;
ntt(h,k,-1);
for(int i=0;i<k;i++) q[i]=i<=n-m?h[i]:0,h[i]=h2[i]=0;
reverse(q,q+1+n-m);
for(int i=0;i<=m;i++) h[i]=g[i];
for(int i=0;i<=n-m;i++) h2[i]=q[i];
k=2<<__lg(n),ntt(h,k,1),ntt(h2,k,1);
for(int i=0;i<k;i++) h[i]=h[i]*h2[i]%mod;
ntt(h,k,-1);
for(int i=0;i<m;i++) r[i]=mo(f[i]-h[i]+mod);
delete[]h,delete[]h2;
}
常系数齐次线性递推
给定 a_0,a_1,\cdots,a_{k-1}，b_1,b_2,\cdots,b_k。a 序列满足递推式 a_n=\sum_{i=1}^ka_{n-i}b_i\left(n>k\right)，求 a_N。
考虑将目标 a_N 展开，直到只包含 a_0 到 a_{k-1}。
具体来说，从大到小依次将 c_ia_i 替换为 c_i\sum_{j=1}^ka_{n-j}b_j，其中，c_i 是展开到此时时 a_i 的系数。
观察上述过程，发现如果把 a_i 表示为 x^i，对 i 展开其实是将当前多项式减去 c_i\left(x^i-\sum_{j=1}^kb_jx^{i-j}\right)=c_ix^{i-k}\left(x^k-\sum_{j=1}^kb_jx^{k-j}\right)。
因此，整个展开过程其实是减去 \left(\sum_{i=k}^nc_ix_{i-k}\right)\left(x^k-\sum_{j=1}^kb_jx^{k-j}\right)，最终得到一个 k-1 次多项式。
不难发现，该过程的本质为将 x^N 对 x^k-\sum_{j-1}^kb_jx^{k-j} 取模，我们一般把这个作为模数的多项式称为该线性递推的特征多项式。
在实现过程中，因为 N 可能非常大，所以直接将 x^N 取模是不现实的。我们可以使用快速幂，每次将两个多项式相乘再取模。
时间复杂度 \mathcal O\left(k\log_2k\log_2N\right)，常数很大。
Code
void times(long long*a,long long*b,long long*f,long long*inv,int n,int k){
long long*h=new long long[(n<<1)+5];
ntt(a,n<<1,1),ntt(b,n<<1,1);
for(int i=0;i<(n<<1);i++) a[i]=a[i]*b[i]%mod;
ntt(a,n<<1,-1),ntt(b,n<<1,-1);//保证 b 不变
for(int i=0;i<=k-2;i++) h[i]=a[(k<<1)-2-i];
for(int i=k-1;i<(n<<1);i++) h[i]=0;
ntt(h,n<<1,1);
for(int i=0;i<(n<<1);i++) h[i]=h[i]*inv[i]%mod;
ntt(h,n<<1,-1);
reverse(h,h+k-1);
for(int i=k-1;i<(k<<1);i++) h[i]=0;
ntt(h,n<<1,1);
for(int i=0;i<(n<<1);i++) h[i]=h[i]*f[i]%mod;
ntt(h,n<<1,-1);
for(int i=0;i<(k<<1);i++) a[i]=i<k?mo(a[i]-h[i]+mod):0;
delete[]h;
}
long long getval(long long*a,long long*b,int k,long long n){//不保证 b 不变
int m=2<<__lg(k);
long long*f=new long long[(m<<1)+5],*inv=new long long[(m<<1)+5];
long long*ans=new long long[(m<<1)+5],*tmp=new long long[(m<<1)+5];
for(int i=0;i<(m<<1);i++) f[i]=ans[i]=0;
ans[n%k]=1,n/=k;//常数优化
reverse(b,b+k);
for(int i=0;i<k;i++) f[i]=mo(mod-b[i]);
f[k]=1,reverse(f,f+1+k),getinv(f,inv,m);
for(int i=k-1;i<m;i++) inv[i]=0;
reverse(f,f+1+k);
ntt(inv,m<<1,1),ntt(f,m<<1,1);//常数优化
while(n){
if(n&1) times(ans,b,f,inv,m,k);
for(int i=0;i<k;i++) tmp[i]=b[i];
times(b,tmp,f,inv,m,k),n>>=1;
}
long long res=0;
for(int i=0;i<k;i++) res=(res+a[i]*ans[i])%mod;
delete[]f,delete[]inv,delete[]ans,delete[]tmp;
return res;
}
泰勒展开
对于函数 f\left(x\right)，求一个与 f\left(x\right) 近似的多项式函数 g\left(x\right)。
首先，我们规定 g 和 f 在 x=0 处相等，于是我们得到了一个零次多项式。
然后，我们觉得 g 跟 f 相差太大，所以我们要求 g' 和 f’ 在 x=0 处相等，于是我们得到了一个一次多项式。
然后，我们觉得 g 跟 f 相差还是太大，所以我们又要求 g'' 和 f'' 在 x=0 处相等，于是我们得到了一个二次多项式。
……
这样一直展开下去，我们就得到了一个与 f 相等的函数。形式化的（令 f^k 表示 f 的 k 阶导数）：
f=\sum_{i=0}^\infty\dfrac{f^i\left(0\right)}{i!}x^i
式子中除以 i! 表示再将 i 次求导得到的函数值积分回去成为 x^i 的系数。
更一般的，我们可以不在点 0 展开，如果我们在 x_0 处展开，则有：
f=\sum_{i=0}^\infty\dfrac{f^i\left(x_0\right)}{i!}\left(x-x_0\right)^i
这就是泰勒展开式，需要说明的是，x_0 可以不是一个数，是一个多项式也可以，我们在后面会用到。
牛顿迭代法
对于函数 f\left(x\right)，快速求出方程 f\left(x\right)=0 相近解，需要注意的是，并不是所有函数的根都可以用牛顿迭代法求解。
首先选取一个 x_0，然后求出 f\left(x\right) 在 \left(x_0,f\left(x_0\right)\right) 处的切线，然后取该切线的根作为新的 x_0。
形式化的，每次迭代执行下式（x_0' 为新的 x_0）：
x_0'=x_0-\dfrac{f\left(x\right)}{f'\left(x\right)}
具体的细节这里就不展开说明了，我们的重点在多项式牛顿迭代上。
多项式牛顿迭代
与普通牛顿迭代类似，可以快速求出方程 f(g)=0 在模 x 的某次方下的解，其中，g 是一个多项式。
设 g_k 表示方程的解 g 对 x^k 取模的结果，将 f 在 g_k 处泰勒展开，并带入 g_{2k}，得：
f\left(g_{2k}\right)=\sum_{i=0}^\infty\dfrac{f^i\left(g_k\right)}{i!}\left(g_{2k}-g_k\right)^i
等式两边对 x^{2k} 取模，得：
0\equiv f\left(g_k\right)+f'\left(g_k\right)\left(g_{2k}-g_k\right)\pmod{x^{2k}}
左侧为 0 是我们的定义，右侧当 g_{2k}-g_k 次数 >1 时，因为 g_{2k}-g_k\equiv 0\pmod{x^k}，所以其对 x^{2k} 取模的值为 0。
\begin{aligned}
f'\left(g_k\right)\left(g_{2k}-g_k\right)&\equiv-f\left(g_k\right)&\pmod{x^{2k}}\\
g_{2k}-g_k&\equiv-\dfrac{f\left(g_k\right)}{f'\left(g_k\right)}&\pmod{x^{2k}}\\
g_{2k}&\equiv g_k-\dfrac{f\left(g_k\right)}{f'\left(g_k\right)}&\pmod{x^{2k}}
\end{aligned}
有了这个式子，只要我们能求出 g_1 和 f'，就可以快速计算出 g 对 x 的某次方取模的值了。
一些例子
刚才的过程可能非常抽象，其实可以直接把最后的式子当成结论来记，接下来我们举一些例子来验证其正确性。
多项式求逆
对多项式 h 求逆，设 f\left(g\right)=\dfrac 1g-h，代入公式，并进行推导：
\begin{aligned}
g_{2k}&\equiv g_k-\dfrac{\frac 1{g_k}-h}{-\frac 1{g_k^2}}&\pmod{x^{2k}}\\
g_{2k}&\equiv g_k+g_k^2\left(\dfrac 1{g_k}-h\right)&\pmod{x^{2k}}\\
g_{2k}&\equiv g_k+g_k-g_k^2h&\pmod{x^{2k}}\\
g_{2k}&\equiv g_k\left(2-g_kh\right)&\pmod{x^{2k}}
\end{aligned}
与我们之前的结论相符。
多项式开根
对多项式 h 开根，设 f\left(g\right)=g^2-h，代入公式，并进行推导：
\begin{aligned}
g_{2k}&\equiv g_k-\dfrac{g_k^2-h}{2g_k}&\pmod{x^{2k}}\\
g_{2k}&\equiv g_k-\dfrac{g_k}2+\dfrac h{2g_k}&\pmod{x^{2k}}\\
g_{2k}&\equiv\dfrac{g_k}2+\dfrac h{2g_k}&\pmod{x^{2k}}
\end{aligned}
也与我们之前的结论相符。
多项式 \ln & \exp
多项式 \ln
对于 n-1 次多项式 f，求多项式 g=\ln f 对 x^n 取模得结果。
对式子两边求导，得 g'=\dfrac{f'}f，直接对 f 求导和求逆，相乘后积分就可以了。
但是有一个小问题，积分无法凭空积出一个常数项，那么我们考虑用定义计算常数项。
因此多项式 $\ln$ 要求目标多项式得常数项为 $1$，并且结果多项式的常数项一定为 $0$。
时间复杂度 $\mathcal O\left(n\log_2n\right)$，常数较大。
### Code
```cpp
void getln(long long*f,long long*g,int n){//n 为 2 的幂次
long long*h=new long long[(n<<1)+5];
getinv(f,g,n);
for(int i=0;i<n;i++) h[i]=f[i+1]*(i+1)%mod;
for(int i=n;i<(n<<1);i++) h[i]=0;
ntt(g,n<<1,1),ntt(h,n<<1,1);
for(int i=0;i<(n<<1);i++) g[i]=g[i]*h[i]%mod;
ntt(g,n<<1,-1);
for(int i=n-1;i;i--) g[i]=g[i-1]*pow(i,mod-2)%mod;
for(int i=n;i<(n<<1);i++) g[i]=0;
g[0]=0,delete[]h;
}
```
## 多项式 $\exp
对于 n-1 次多项式 f，求多项式 g=\exp f 对 x^n 取模得结果。
设函数 h\left(g\right)=\ln g-f，套牛顿迭代公式，得：
\begin{aligned}
g_{2k}&\equiv g_k-\dfrac{\ln g_k-f}{\frac 1{g_k}}&\pmod{x^{2k}}\\
g_{2k}&\equiv g_k-g_k\left(\ln g_k-f\right)&\pmod{x^{2k}}\\
g_{2k}&\equiv g_k\left(1+f-\ln g_k\right)&\pmod{x^{2k}}
\end{aligned}
这样，我们就可以快速计算多项式 \exp 了。但是与多项式 \ln 类似，多项式 \exp 存在初值的问题，并且其运算过程中也涉及多项式 \ln，所以多项式 \exp 要求目标函数常数项为 0，并且结果多项式的常数项一定为 1。
时间复杂度 \mathcal O\left(n\log_2n\right)，常数很大。
Code
void getexp(long long*f,long long*g,int n){//n 为 2 的幂次
long long*h=new long long[(n<<1)+5];//不需要清空，因为在 getln -> getinv 中已清空
for(int i=0;i<(n<<1);i++) g[i]=0;
g[0]=1;
for(int i=2;i<=n;){
getln(g,h,i);
h[0]=1;//f[0]=ln(g)[0]=0
for(int j=1;j<i;j++) h[j]=mo(f[j]-h[j]+mod);
i<<=1,ntt(g,i,1),ntt(h,i,1);
for(int j=0;j<i;j++) g[j]=g[j]*h[j]%mod;
ntt(g,i,-1);
for(int j=i>>1;j<i;j++) g[j]=0;
}
delete[]h;
}
多项式快速幂
对于 n-1 次多项式 f 和非负整数 k，求 f^k\bmod x^n。
由于 NTT 多项式乘法的本质是循环卷积，所以把 f 变换过后把每一位 k 次方再逆变换回去是错误的。
直接 \mathcal O\left(n\log_2n\log_2k\right) 快速幂复杂度不够优秀（通过下面的过程可以发现把 k 对 p\left(p-1\right) 取模也可以），这里也不考虑。
那么我们可以考虑利用关系式 \ln f^k=k\ln f，所以直接把 f 先 \ln 然后乘 k 再 \exp 回去就可以了，k 可以先对模数取模。
还有一个小问题，多项式 \ln 要求 f 的常数项为 1，但实际情况中可能 f 的常数项不为 1，甚至为 0。
其实也不难解决，先将 f 表示为 ax^bg，其中，a,b 是常数，g 是常数项为 1 的多项式。不难证明，只要 f\ne 0，这样的形式总是存在。
然后分别把三部分 k 次方，最后再乘起来即可。a^k 直接算；x^b 的 k 次方就是 x^{kb}；g^k 就用之前说的方法做就行了。
时间复杂度 \mathcal O\left(n\log_2n\right)，常数很大。
Code
void getpow(long long*f,long long*g,int n,long long k){//n 为 2 的幂次
if(!k){
g[0]=1;
for(int i=1;i<n;i++) g[i]=0;
return;
}
int num=0;
while(num<n&&!f[num]) num++;
if((k>=n&&num)||num*k>=n){//k 在 long long 范围内，直接乘 num 可能乘爆
for(int i=0;i<n;i++) g[i]=0;
return;
}
long long val=f[num],inv=pow(val,mod-2);
for(int i=num;i<n;i++) f[i-num]=f[i]*inv%mod;
for(int i=n-num;i<n;i++) f[i]=0;
num*=k,val=pow(val,k),k%=mod;
getln(f,g,n);
for(int i=0;i<n;i++) g[i]=g[i]*k%mod;
long long*h=new long long[n+5];
getexp(g,h,n);
for(int i=0;i<num;i++) g[i]=0;
for(int i=num;i<n;i++) g[i]=h[i-num]*val%mod;
delete[]h;
}
\ln & \exp 优化背包
用 \ln 和 \exp 可以用来优化多重背包计数问题（包括有无穷多相同物品的情况）。
我们从完全背包入手。
对于一个体积为 a_i 物品 i，我们将其视作多项式 f_i=\sum_{j=0}^\infty x^{a_ij}。发现最终要求恰好装满容积为 m 的背包的方案数就是所有物品对应多项式的乘积的 m 次项系数。
但是这样仍然不好计算，直接每次求 \ln 相加后再 \exp 的总时间复杂度均为 \mathcal O\left(nm\log_2m\right)，甚至不如 \mathcal O\left(nm\right) 的 DP。
考虑优化求 \ln 的复杂度。
注意到 f_i=\sum_{j=0}^\infty x^{a_ij}=\dfrac 1{1-x^{a_i}}，套 \ln 的公式，得：
\begin{aligned}
(\ln f_i)’&=\dfrac{f_i'}{f_i}\\
&=\dfrac{\frac{0-\left(-a_ix^{a_i-1}\right)}{\left(1-x^{a_i}\right)^2}}{\frac 1{1-x^{a_i}}}\\
&=\dfrac{a_ix^{a_i-1}}{1-x^{a_i}}\\
&=\sum_{j=1}^\infty a_ix^{a_ij-1}\\
\ln f_i&=\sum_{j=1}^\infty\dfrac{a_i}{a_ij}x^{a_ij}\\
&=\sum_{j=1}^\infty\dfrac 1jx^{a_ij}
\end{aligned}
可以对每个体积，统计有多少 a_i 恰好是这么多，然后对每个权值暴力加，时间复杂度为 \mathcal O\left(\sum_{i=1}^n\left\lfloor\dfrac mi\right\rfloor\right)=\mathcal O\left(m\ln m\right)。
这样，我们就会做完全背包了。考虑拓展到多重背包。
如果体积为 a_i 的物品有 b_i 个，那么它的多项式应为 f_i=\sum_{j=0}^{b_i}x^{a_ij}，没有完全背包那么好的性质。
考虑对原式进行一点变形，并进行推导：
\begin{aligned}
f_i&=\sum_{j=0}^\infty x^{a_ij}-\sum_{j=b_i+1}^\infty x^{a_ij}\\
&=\dfrac 1{1-x^{a_i}}-x^{a_i\left(b_i+1\right)}\dfrac 1{1-x^{a_i}}\\
&=\left(1-x^{a_i\left(b_i+1\right)}\right)\dfrac 1{1-x^{a_i}}
\end{aligned}
因为要求 $\exp$，总时间复杂度 $\mathcal O(n\log_2n)$。
### Code
```cpp
void knapsack(int n,int m,int*a,int*b,long long*ans){//b[i]=0 表示第 i 个物品有无穷多个
int*num=new int[m+5],k=2<<__lg(m);
long long*inv=new long long[m+5],*f=new long long[k+5];
for(int i=1;i<=m;i++) num[i]=0,inv[i]=pow(i,mod-2);
for(int i=0;i<k;i++) f[i]=0;
for(int i=1;i<=n;i++){
if(a[i]<=n){
num[a[i]]++;
if(b&&a[i]*(b[i]+1)<=m) num[a*(b+1)]--;
}
}
for(int i=1;i<=m;i++){
long long tmp=mo(num[i]+mod);
for(int j=i;j<=m;j+=i) f[j]=(f[j]+tmp*inv[j/i])%mod;
}
getexp(f,ans,k);
delete[]inv,delete[]f;
}
```
## $\ln$ & $\exp$ 的组合意义
直接分析不是很好分析，我们结合例子来分析。
设 $f_i$ 表示 $i$ 个点的有标号连通简单图数量，$F_i$ 表示 $i$ 个点的有标号简单图数量，定义 $f_0=F_0=0$。
首先，我们可以写出由 $f$ 到 $F$ 的递推式：
$$
F_k=\sum_{i=1}^k\dfrac1{i!}\sum_{\sum_{j=1}^ia_j=k}\dfrac{k!}{\prod_{j=1}^ia_j!}\prod_{j=1}^if_{a_j}
$$
让我们来分析一下这个式子的意义：
- $\sum_{i=1}^k\sum_{\sum_{j=1}^ia_j=k}$ 表示枚举所有将 $k$ 顶点划分成连通块的方案，$i$ 表示划分成几个连通块，$a$ 序列表示每个连通块有多少点。注意此处限制 $a_j\ge1$。
- $\prod_{j=1}^if_{a_j}$ 表示计算直接连起来的方案数。
- $\dfrac{k!}{\prod_{j=1}^ia_j!}$ 表示将所有顶点重新按原顺序标号的方案数（可重集排列）。
- $\dfrac1{i!}$ 表示最后将所有连通块无序化。注意此处因为已经对顶点重标号了，所以就算 $a$ 中有重复元素也能正确计算。
接下来我们稍微化简一下这个式子：
$$
\begin{aligned}
\dfrac{F_k}{k!}&=\sum_{i=1}^k\dfrac1{i!}\sum_{\sum_{j=1}^ia_j=k}\prod_{j=1}^i\dfrac{f_{a_j}}{a_j!}\\
\dfrac{F_k}{k!}&=\sum_{i=0}^\infty\dfrac1{i!}\sum_{\sum_{j=1}^ia_j=k}\prod_{j=1}^i\dfrac{f_{a_j}}{a_j!}
\end{aligned}
$$
此时我们可以仅要求 $a_j\ge0$，因为如果有 $a_j=0$，那么 $f_{a_j}=0$，所以后面的乘积一定为 $0$。将 $\sum_{i=1}^k$ 改为 $\sum_{i=0}^\infty$ 也是不影响答案的，因为如果 $i=0$，那么不存在满足条件的 $a$ 数组；如果 $i>k$，必然有一个 $a_j=0$，所以内层 $\sum$ 的结果一定为 $0$。
先把这个式子放在这里，接下来，我们从多项式 $\exp$ 的角度来推一下式子。
记 $\hat f$ 表示 $f$ 的指数生成函数（EGF），根据定义，有：
$$
\begin{aligned}
\hat f=\sum_{i=0}^\infty\dfrac{f_i}{i!}x^i\\
\hat F=\sum_{i=0}^\infty\dfrac{F_i}{i!}x^i
\end{aligned}
$$
将 $\exp$ 函数在 $0$ 处泰勒展开，并带入 $\hat f$，得（令 $g_i$ 表示多项式 $g$ 的 $i$ 次项系数）：
$$
\begin{aligned}
\exp\hat f&=\sum_{i=0}^\infty\dfrac 1{i!}\hat f^i\\
\left(\exp\hat f\right)_k&=\sum_{i=0}^\infty\dfrac 1{i!}\sum_{\sum_{j=1}^ia_j=k}\prod_{j=1}^i\dfrac{f_{a_j}}{a_j!}
\end{aligned}
$$
根据之前的式子，我们不难发现，$\exp\hat f=\hat F$。
所以，我们可以得出多项式 $\exp$ 的组合意义（$\ln$ 反之）：由有标号单组划分方案数的 EGF 进行 $\exp$ 之后就可以得到有标号多组划分方案数的 EGF。
# 多项式多点求值
给定一个 $n$ 次多项式 $f$ 和 $m$ 个整数 $x_0,x_1,\cdots,x_{m-1}$，求 $f\left(x_0\right),f\left(x_1\right),\cdots,f\left(x_{m-1}\right)$。
考虑分治，假设当前在区间 $\left[l,r\right]$，考虑对 $f$ 降次。
注意到把 $f$ 对 $g=\prod_{i=l}^rx-x_i$ 取模是不影响答案的，因为 $g$ 在 $x_l,x_{l+1},\cdots,x_r$ 处的取值都是 $0$。
这样，当我们在一个区间 $[l,r]$ 时，我们可以把 $f$ 控制在 $r-l+1$ 次，保证了时间复杂度。
由于多项式取模常数很大，所以当区间长度比较小（大概在 $\le50$ 左右）的时候，直接暴力反而会跑得更快。
实现时跑两遍分治，第一遍处理区间的 $g$，第二遍处理答案。
总时间复杂度为 $\mathcal O\left(m\log_2^2m+n\log_2n\right)$，常数很大。
### Code
```cpp
#define ls(v) v<<1
#define rs(v) v<<1|1
void solve(long long*f,int n,long long*x,long long*y,int m){//循环实现分治，需要 m 为 2 的幂次，答案在 y 数组
long long **g=new long long*[(m<<1)+5],*h=new long long[m+5],*h2=new long long[m+5];
bool *vis=new bool[(m<<1)+5];
for(int i=(m<<1)-1;i;i--){
vis[i]=0;
int l=(i^(1<<__lg(i)))<<(__lg(m)-__lg(i)),r=l+(1<<(__lg(m)-__lg(i)))-1,len=r-l+1;
if(i==1){//保证第一次取模不会 RE
g[i]=new long long[max(len,2<<__lg(n))+5];
for(int j=2<<__lg(n);j>len;j--) g[i][j]=0;
}
else g[i]=new long long[len+2];
if(l==r){
g[i][0]=mo(mod-x[i-m]),g[i][1]=1;
continue;
}
for(int j=0;j<=(len>>1);j++) h[j]=g[ls(i)][j],h2[j]=g[rs(i)][j];
for(int j=(len>>1)+1;j<len;j++) h[j]=h2[j]=0;
ntt(h,len,1),ntt(h2,len,1);
for(int j=0;j<len;j++) h[j]=h[j]*h2[j]%mod;
ntt(h,len,-1);
g[i][0]=g[ls(i)][0]*g[rs(i)][0]%mod;
g[i][len]=g[ls(i)][len>>1]*g[rs(i)][len>>1]%mod;
for(int j=1;j<len;j++) g[i][j]=h[j];
}
division(f,n,g[1],m,h);//只要余数，不要商，需要稍微改一下 division 函数
for(int i=0;i<=m;i++) g[1][i]=h[i];
for(int i=2;i<(m<<1);i++){
if(vis[i]){
if(i<m) vis[ls(i)]=vis[rs(i)]=1;
continue;
}
int l=(i^(1<<__lg(i)))<<(__lg(m)-__lg(i)),r=l+(1<<(__lg(m)-__lg(i)))-1,len=r-l+1;
division(g[i>>1],(len<<1)-1,g[i],len,h);//同上
if(len<=50){
if(i<m) vis[ls(i)]=vis[rs(i)]=1;
for(int j=l;j<=r;j++){
long long now=1;
y[j]=0;
for(int k=0;k<len;k++,now=now*x[j]%mod) y[j]=(y[j]+h[k]*now)%mod;
}
}
else for(int j=0;j<len;j++) g[i][j]=h[j];
}
for(int i=1;i<(m<<1);i++) delete[]g[i];
delete[]g,delete[]h,delete[]h2,delete[]vis;
}
#undef ls
#undef rs
```
## 一种特殊形式
给定一个 $n$ 次多项式 $f$ 和 $m$ 以及常数 $c$，求 $f\left(c^0\right),f\left(c^1\right),\cdots,f\left(c^{m-1}\right)$。
注意到答案 $f\left(c^i\right)=\sum_{j=0}^n f_jc^{ij}$，其中，$f_j$ 表示 $f$ 的 $j$ 次项系数。式子中，$i,j$ 有乘法关系，不好处理，考虑拆成加减法进行卷积。
考虑式子 $ij=\binom{i+j}2-\binom i2-\binom j2$，其中 $\binom ab$ 表示组合数。这样，原式就变成了 $f\left(c^i\right)=c^{-\binom i2}\sum_{i=0}^nf_jc^{-\binom j2}c^\binom{i+j}2$，这个式子里面只有关于 $j$ 和 $i+j$ 的式子，就可以用卷积处理了，时间复杂度 $\mathcal O\left(n\log_2n\right)$。
当模数比较小的时候，也可以找到一个原根，然后用这种方法求出所有的函数值（需要 MTT）。
### Code
```cpp
void solve(long long*f,int n,long long c,int m,long long*y){
int k=2<<__lg(n+n+m);
long long sum=1,now=1,inv=pow(c,mod-2),*g=new long long[k+5],*h=new long long[k+5];
for(int i=0;i<k;i++) g[i]=h[i]=0;
for(int i=0;i<=n;i++) g[n-i]=f[i]*sum%mod,sum=sum*now%mod,now=now*inv%mod;//递推幂次，减小常数
sum=now=1;
for(int i=0;i<n+m;i++) h[i]=sum,sum=sum*now%mod,now=now*c%mod;
ntt(g,k,1),ntt(h,k,1);
for(int i=0;i<k;i++) g[i]=g[i]*h[i]%mod;
ntt(g,k,-1);
sum=now=1;
for(int i=0;i<m;i++) y[i]=g[i+n]*sum%mod,sum=sum*now%p,now=now*inv%p;
delete[]g,delete[]h;
}
```
# 多项式快速插值
给出 $n$ 个点 $\left(x_0,y_0\right),\left(x_1,y_1\right),\cdots,\left(x_{n-1},y_{n-1}\right)$，求一个 $n-1$ 次多项式 $f$，使得对于任意 $0\le i<n$，有 $f\left(x_i\right)=y_i$。
考虑分治，假设当前在区间 $\left[l,r\right]$，合并左右区间的 $f_{\left[l,mid\right]},f_{\left(mid,r\right]}$（以下简写为 $f_l,f_r$），得到当前 $f$。
考虑构造 $g_l,g_r$，使得 $g_l$ 是仅在左区间的 $x$ 处取值均为 $0$ 的 $mid-l+1$ 次多项式，$g_r$ 是仅在右区间的 $x$ 处取值均为 $0$ 的 $r-mid$ 次多项式，然后让 $f=f_lg_r+f_rg_l$。
令 $g_l=\prod_{i=l}^{mid}x-x_i,g_r=\prod_{i=mid+1}^rx-x_i$。但是我们发现这样定义的 $g_l,g_r$ 会使原本 $f_l,f_r$ 正确的函数值变化，并且无法避免。
考虑在递归时就将 $g$ 带来的变化量放到要求的函数值上，这样乘上 $g$ 变化回来时会将值修正而不是变错。
具体的，在当前区间，将 $f\left(x_l\right),f\left(x_{l+1}\right),\cdots,f\left(x_{mid}\right)$ 分别除以 $g_r\left(x_l\right),g_r\left(x_{l+1}\right),\cdots,g_r\left(x_{mid}\right)$ 再下传，右侧同理。
发现求所有 $g\left(x\right)$ 需要多点求值，所以直接这么做复杂度是 $\mathcal O\left(n\log_2^3n\right)$ 的，不够优秀，考虑优化。
观察这个过程，发现对于每一个 $x_i$，其实最终它是除以了 $\prod_{j\ne i}x_i-x_j$。
如果没有 $j\ne i$ 的限制，这个式子是好处理的，但如果去掉这个限制，函数值一定为 $0$，并且需要除以 $0$ 来还原，不可行。
构造多项式 $t=\sum_{i=1}^n\prod_{j\ne i}x-x_j$，对其多点求值。只有当 $\sum$ 中的 $i$ 是我们需要的 $i$ 时，$\prod$ 的值不为 $0$，并且就是我们需要的值。
$t$ 可以在分治中维护出来，但是有更简单的求法，可以发现 $t=g'$，这样，求出 $g$ 后直接求导就可以得到 $t$。
如果递归实现，可以把多点求值的分治和合并答案的分治合并起来，只需要分治两次。第一次处理 $g$（这里的 $g$ 恰好与求值需要的 $g$ 相同），第二次先求值再合并（注意求值可以在小数据范围下暴力，插值必须从最底层开始）。
时间复杂度 $\mathcal O\left(n\log_2^2n\right)$，常数很大。
### Code
```cpp
#define ls(v) v<<1
#define rs(v) v<<1|1
void solve(int n,long long*x,long long*y,long long*f){
int m=2<<__lg(n);
long long **g=new long long*[(m<<1)+5],**t=new long long*[(m<<1)+5];
long long *h=new long long*[(m<<1)+5],*h2=new long long*[(m<<1)+5];
bool *vis=new bool[(m<<1)+5];
for(int i=(m<<1)-1;i;i--){
int l=(i^(1<<__lg(i)))<<(__lg(m)-__lg(i)),r=l+(1<<(__lg(m)-__lg(i)))-1,len=r-l+1;
g[i]=new long long[len+2],t[i]=new long long[len+2],vis[i]=0;
if(l>=n){
g[i][0]=1;
for(int j=1;j<=len;j++) g[i][j]=0;
continue;
}
if(l==r){
g[i][0]=mo(mod-x[i-m]),g[i][1]=1;
continue;
}
for(int j=0;j<=(len>>1);j++) h[j]=g[ls(i)][j],h2[j]=g[rs(i)][j];
for(int j=(len>>1)+1;j<len;j++) h[j]=h2[j]=0;
ntt(h,len,1),ntt(h2,len,1);
for(int j=0;j<len;j++) h[j]=h[j]*h2[j]%mod;
ntt(h,len,-1);
g[i][0]=g[ls(i)][0]*g[rs(i)][0]%mod;
g[i][len]=g[ls(i)][len>>1]*g[rs(i)][len>>1]%mod;
for(int j=1;j<len;j++) g[i][j]=h[j];
}
for(int i=0;i<m;i++) t[1][i]=g[1][i+1]*(i+1)%mod;
for(int i=2;i<(m<<1);i++){
if(vis[i]){
if(i<m) vis[ls(i)]=vis[rs(i)]=1;
continue;
}
int l=(i^(1<<__lg(i)))<<(__lg(m)-__lg(i)),r=l+(1<<(__lg(m)-__lg(i)))-1,len=r-l+1;
if(l>=n) continue;
division(t[i>>1],g[i],min(len<<1,n-(l&(INT_MAX^len)))-1,min(len,n-l),h);
if(min(len,n-l+1)<=50){
r=min(r,n-1),len=r-l+1;
vis[ls(i)]=vis[rs(i)]=1;
for(int j=l;j<=r;j++){
long long val=0,now=1;
for(int k=0;k<len;k++,now=now*x[j]%mod){
val=(val+now*h[k])%mod;
}
y[j]=y[j]*pow(val,mod-2)%mod;
}
}
else for(int j=0;j<len;j++) t[i][j]=h[j];
}
for(int i=(m<<1)-1;i;i--){
int l=(i^(1<<__lg(i)))<<(__lg(m)-__lg(i)),r=l+(1<<(__lg(m)-__lg(i)))-1,len=r-l+1;
if(l>=n){
for(int j=0;j<len;j++) t[i][j]=0;
continue;
}
if(l==r){
t[i][0]=y[i-m];
continue;
}
for(int j=0;j<(len>>1);j++) h[j]=t[ls(i)][j],h2[j]=g[rs(i)][j];
for(int j=len>>1;j<len;j++) h[j]=h2[j]=0;
h2[len>>1]=g[rs(i)][len>>1];
ntt(h,len,1),ntt(h2,len,1);
for(int j=0;j<len;j++) h[j]=h[j]*h2[j]%mod;
ntt(h,len,-1);
for(int j=0;j<len;j++) t[i][j]=h[j];
for(int j=0;j<(len>>1);j++) h[j]=t[rs(i)][j],h2[j]=g[ls(i)][j];
for(int j=len>>1;j<len;j++) h[j]=h2[j]=0;
h2[len>>1]=g[ls(i)][len>>1];
ntt(h,len,1),ntt(h2,len,1);
for(int j=0;j<len;j++) h[j]=h[j]*h2[j]%mod;
ntt(h,len,-1);
for(int j=0;j<len;j++) t[i][j]=mo(t[i][j]+h[j]);
}
for(int i=0;i<n;i++) f[i]=t[1][i];
for(int i=0;i<(m<<1);i++) delete[]g[i],delete[]t[i];
delete[]g,delete[]t,delete[]h,delete[]h2,delete[]vis;
}
#undef ls
#undef rs
```
# 上升幂多项式 & 下降幂多项式
OI 中唯一的作用是用来快速计算第一类 Stirling 数。
记 $f_i$ 表示第一类 Stirling 数的第 $i$ 行的的普通生成函数（OGF），则：
$$
f_i=\sum_{j=0}^ix^js_{i,j}
$$
写出第一类 Stirling 数的递推公式，并进行推导（$f_{i,j}$ 表示 $f_i$ 的 $j$ 次项系数）：
$$
\begin{aligned}
s_{n,m}&=s_{n-1,m-1}+\left(n-1\right)s_{n-1,m}\\
f_{n,m}&=f_{n-1,m-1}+\left(n-1\right)f_{n-1,m}\\
f_{n,m}x^m&=xf_{n-1,m-1}x^{m-1}+\left(n-1\right)f_{n-1,m}x^m\\
f_n&=\left(x+n-1\right)f_{n-1}
\end{aligned}
$$
由此，我们可以写出 $f$ 的通项公式：
$$
f_n=\prod_{i=0}^{n-1}x+i=x^{\overline{n}}
$$
接下来，我们考虑如何快速计算 $x^{\overline{n}}$。
考虑倍增。首先，$+1$ 操作是容易的，直接递推即可，单次时间复杂度线性。
考虑 $\times 2$ 操作，假设我们当前知道 $x^{\overline{m}}$，如何求出 $x^{\overline{2m}}$。
注意到 $x^{\overline{2m}}=x^{\overline m}\left(x+m\right)^{\overline{m}}$，那么只要我们求出了 $\left(x+m\right)^{\overline{m}}$，就可以在 $\mathcal O\left(m\log_2m\right)$ 的时间内求出 $x^{\overline{2m}}$。
将 $x=x+m$ 带入 $x^{\overline{m}}$，并进行推导（$a_i$ 表示 $x^{\overline{m}}$ 的 $i$ 次项系数）：
$$
\begin{aligned}
\left(x+m\right)^{\overline{m}}&=\sum_{i=0}^ma_i\left(x+m\right)^i\\
&=\sum_{i=0}^ma_i\sum_{j=0}^ix^jm^{i-j}\binom ij\\
&=\sum_{i=0}^mx^i\sum_{j=i}^ma_jm^{j-i}\binom ji\\
&=\sum_{i=0}^m\dfrac{x^i}{i!}\sum_{j=i}^ma_jj!\dfrac{m^{j-i}}{\left(j-i\right)!}
\end{aligned}
$$
内层 $\sum$ 是一个卷积的形式，可以直接用多项式乘法处理，这样就可以处理 $\times 2$ 操作了。
注意到两种递推推出 $x^{\overline{2m}}$ 的 $x^k$ 次项系数都只用到了 $x^{\overline{m}}$ 的 $x^0,x^1,\cdots,x^k$ 系数，所以这个做法也可以处理 $n$ 很大时的前缀第一类 Stirling 数。
单次时间复杂度 $\mathcal O\left(m\log_2m\right)$，总时间复杂度 $\mathcal O(n\log_2n)$，处理前 $k$ 个时间复杂度 $\mathcal O(k\log_2 k\log_2n)$（因为有多项式乘法）。
### Code
```cpp
void Stirling(int n,long long*f){
if(n==0) return f[0]=1,void();
int k=2<<__lg(n);
long long *g=new long long[k+5],*h=new long long[k+5];
long long *fac=new long long[n+5],*inv=new long long[n+5];
fac[0]=1;
for(int i=1;i<=n;i++) fac[i]=fac[i-1]*i%mod;
inv[n]=pow(fac[n],mod-2);
for(int i=n-1;~i;i--) inv[i]=inv[i+1]*(i+1)%mod;
f[0]=0,f[1]=1;
for(int i=2;i<k;i++) f[i]=0;
for(int i=__lg(n)-1;~i;i--){
int m=n>>(i+1),k=4<<__lg(m);
long long mi=1;
for(int j=0;j<=m;j++,mi=mi*m%mod) g[j]=f[m-j]*fac[m-j]%mod,h[j]=mi*inv[j]%mod;
for(int j=m+1;j<k;j++) g[j]=h[j]=0;
ntt(g,k,1),ntt(h,k,1);
for(int j=0;j<k;j++) g[j]=g[j]*h[j]%mod,h[j]=0;
ntt(g,k,-1);
for(int j=0;j<=m;j++) h[j]=g[m-j]*inv[j]%mod;
ntt(f,k,1),ntt(h,k,1);
for(int j=0;j<k;j++) f[j]=f[j]*h[j]%mod;
ntt(f,k,-1);
if(n&(1<<i)){
m=n>>i;
for(int j=m;~j;j--) f[j]=(f[j-1]+f[j]*(m-1))%mod;
}
}
delete[]g,delete[]h,delete[]fac,delete[]inv;
}
```