大砖数学（三）生成函数和多项式全家桶
Amidst
·
·
算法·理论
前言
行列式被推到第四篇去了，这里是插播。
后缀算法原定于本篇之前更新，但是好像赶不出进度，就先搁置了。
我爱积积本来要更新的，但是本来本来，没有本来的啦！
后记：后缀算法更完了，这篇还没开始写。
很多东西可能来自这篇，给 wangyizhi 一个二作，但是他退役了，那我就放心贺了。
生成函数
生成函数部分为幼儿园小班难度，如果您已经上过幼儿园小班，建议跳过。
我不知道怎么开头。
生成函数（Generating Function，又称母函数），是利用多项式的各项指数和系数来表示状态和结果的函数。
骗你的，我也不知道这东西标准定义是啥。
我们可以利用指数来表示状态，用函数之间的运算来表示状态的转移，最终对应项的系数就是我们想要的结果，一般来说，这个状态是一个数。
普通生成函数
普通生成函数（OGF）可以把数列转化成多项式。具体的，序列 \{a_n\} 的 OGF 为
F(x)=\sum_{i=0}^n a_ix_i
如果数列无限长，生成函数就变成了幂级数，例如
F(x)=1+x+x^2+\cdots
幂级数有封闭形式，即形式幂级数意义下的恒等，以 F(x)=1+x+x^2+\cdots 为例：
F(x)&=1+x+x^2+\cdots \\
xF(x)&=x+x^2+x^3+\cdots \\
(1-x)F(x)&=1 \\
F(x)&=\frac{1}{1-x}
\end{aligned}
生成函数一个比较重要的性质是，它的乘法相当于序列的卷积，即
\operatorname{OGF}(a*b)=\operatorname{OGF}(a)\times \operatorname{OGF}(b)
其中
(a*b)_i=\sum_{j=0}^ia_jb_{i-j}
生成函数可以用作有递推式的情况下求通项。例如，斐波那契数列：
我们有递推式 f_i=f_{i-1}+f_{i-2}，不妨设 F(x)=\operatorname{OGF}(f)，有
F(x)=xF(x)+x^2F(x)+x
解得
F(x)=\frac{x}{1-x-x^2}
设
F(x)=\frac A{1-ax} + \frac B{1-bx}=\frac{(A+B)-(Ab+Ba)x}{1-(a+b)x+abx^2}
则
\begin{cases}
A+B=0 \\
Ab+Ba=-1 \\
a+b=1 \\
ab=-1
\end{cases}
解得
\begin{cases}
A=\frac 1 {\sqrt 5} \\
B=-\frac 1 {\sqrt 5} \\
a=\frac{1+\sqrt 5}{2} \\
b=\frac{1-\sqrt 5}{2}
\end{cases}
又由于
\frac{A}{1-ax}=\sum_{i\ge 0}Aa^ix^i
所以
F(x)=\sum_{i=0} \frac 1{\sqrt 5} \left( \left( \frac{1+\sqrt 5}{2}\right)^i - \left( \frac{1-\sqrt 5}{2}\right)^i \right) x^i
故
f_n=\frac 1{\sqrt 5} \left( \left( \frac{1+\sqrt 5}{2}\right)^n - \left( \frac{1-\sqrt 5}{2}\right)^n \right)
指数型生成函数
定义数列 \{a_n\} 的指数型生成函数（EGF）为
\hat F(x)=\sum_{i=0}^n a_i \frac{x^i}{i!}
考虑一个经典问题：
把 n 个有标号的求染成红色或蓝色，求方案数。
设其中有 i 个红球的方案数为 f_i，有 i 个蓝球的方案为 g_i，则答案为
h_n&=n! \sum_{i+j=n} \binom{n}{i} f_ig_j \\
&=n! \sum_{i+j=n} \frac{f_i}{i!} \cdot \frac{g_j}{j!}
\end{aligned}
则
\frac{h_n}{n!}=\sum_{i+j=n} \frac{f_i}{i!} \cdot \frac{g_j}{j!}
记 \hat F(x)=\operatorname{EGF}(f),\hat G(x)=\operatorname{EGF}(g),\hat H(x)=\operatorname{EGF}(h)，则
\hat H(x)=\hat F(x) \hat G(x)
结论可以推广到多个数列的情况，我们有
\hat G(x)=\prod_i \hat F_i(x)
和
g_n=\sum_{\sum a_i=n} \binom{n}{a_1,a_2,\dots,a_n}\prod f_{i,a_i}
等价。
另外，对 e^x 做泰勒展开，有
e^x=\sum_{i=0}^\infty \frac{x^i}{i!}=\operatorname{EGF}(\{1,1,1,\dots\})
EGF 的微分（导数）也有其组合意义。
令
\hat F(x)=\sum_{i\ge 0} a_i \frac{x^i}{i!}
有
\hat F'(x)=\sum_{i\ge 1}a_i \frac{ix^{i-1}}{i!}=\sum_{i \ge 0} a_{i+1} \frac{x^i}{i!}
相当于给每项都乘上了 i+1，也就是说，我们可以将其理解为有 i+1 个有标号元素的结构，从 i+1 个元素中选取一个，也就是说，这是在结构中标记了一个元素的方案数，于是我们可以将其用作含有某个特定元素的计数。
二元生成函数
一个变量不够用，就再引入一个，到后面消元。
多项式全家桶
多项式全家桶部分为幼儿园小中班衔接难度，如果您已经上过幼儿园中班或上过小中班衔接，建议跳过。
任意模数多项式乘法
我不会啊。
但我们有三模 NTT。
泰勒展开
设函数 f 在点 a 的某邻域内具有足够阶导数，则：
f(x)=\sum_{k=0}^n \frac{f^{(k)}(a)}{k!}(x-a)^k+o((x-a)^n)
邻域，即 (a-\epsilon,a+\epsilon)。
这里的 o((x-a)^n) 表示比 (x-a)^n 小一个量级。
多项式牛顿迭代
已知二元函数 G(x,y)，求一多项式 G(x,f(x))\equiv 0 \pmod {x^n}。
考虑倍增，先对 n=1 求出答案。
现在假设我们求出了模 x^k 意义下的解 f'，将 G(x,f(x)) 对 f(x) 在 f(x)=f'(x) 处做一阶泰勒展开，令 f(x)=f'(x)+\Delta(x)，显然 \Delta(x) 是 x^k 量级，有
G(x,f'+\Delta)=G(x,f')+\frac{\partial G}{\partial y}(x,f'(x))\cdot \Delta+O(\Delta^2)
忽略 $O(\Delta^2)$，这样精度在 $x^{2k}$，即精度平方，并令 $G(x,f(x))$ 为 $0$：
$$\begin{aligned}
G(x,f')+G_y(x,f')\cdot \Delta &= 0 \\
\Delta &= -\frac{G(x,f')}{G_y(x,f')}
\end{aligned}$$
于是
$$f(x) \equiv f'(x)-\frac{G(x,f')}{G_y(x,f')} \pmod {x^{2k}}$$
这里用 $G_y$ 表示了上述偏导。
如此迭代，$O(\log n)$ 次即可求出 $f(x)$。
### 多项式乘法逆
使用多项式牛顿迭代。
设给定函数为 $h(x)$，则 $G(x,y)=\frac{1}{y}-h(x)$。设 $f'(x)$ 为求出的模 $x^{\lceil \frac{n}{2} \rceil}$ 意义下的答案，我们有
$$\begin{aligned}
f(x)&\equiv f'(x)-\frac{\frac{1}{f'(x)}-h(x)}{-\frac{1}{f'^2(x)}} \\
&\equiv 2f'(x) - f'^2(x)h(x) \pmod {x^n}
\end{aligned}$$
考虑时间复杂度。由 master 定理，每次都要做 FFT，得到
$$T(n)=T(\frac{n}{2})+O(n\log n)=O(n\log n)$$
总时间复杂度 $O(n\log n)$。
:::info[参考代码]
```cpp
void invNTT(int f[],int ans[],int tmpa[],int rev[],int n,int mod,int G)
{
if(n==1)
{
ans[0]=quick_pow(f[0],mod-2,mod);
return;
}
invNTT(f,ans,tmpa,rev,n+1>>1,mod,G);
int l=-1,lim=1;
while(lim<=(n<<1))
{
lim<<=1;
l++;
}
for(int i=0;i<=lim;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<l);
for(int i=0;i<n;i++)
tmpa[i]=f[i];
for(int i=n;i<=lim;i++)
tmpa[i]=0;
NTT(tmpa,rev,lim,1,mod,G);
NTT(ans,rev,lim,1,mod,G);
for(int i=0;i<=lim;i++)
ans[i]=(2-tmpa[i]*ans[i]%mod+mod)%mod*ans[i]%mod;
NTT(ans,rev,lim,-1,mod,G);
for(int i=n;i<=lim;i++)
ans[i]=0;
}
```
:::
### 多项式 ln
设给定的函数为 $f(x)$，考虑先微分再积分，有
$$\begin{aligned}
\frac{\text{d} \ln f(x)}{\text{d} x}&\equiv \frac{f'(x)}{f(x)} \\
\text{d} \ln f(x) &\equiv \frac{f'(x)}{f(x)} \text{d} x \\
\ln f(x)&\equiv \int \frac{f'(x)}{f(x)} \text{d} x \pmod {x^n}
\end{aligned}$$
求导求逆乘法积分，做完了，时间复杂度 $O(n\log n)$。
:::info[参考代码]
```cpp
void devNTT(int f[],int g[],int n,int mod)
{
for(int i=1;i<n;i++)
g[i-1]=i*f[i]%mod;
g[n-1]=0;
}
void lnNTT(int f[],int g[],int n,int mod,int G)
{
memset(a,0,sizeof a);
memset(b,0,sizeof b);
devNTT(f,a,n,mod);
invNTT(f,b,tmpa,rev,n,mod,G);
int l=-1,lim=1;
while(lim<=(n<<1))
{
lim<<=1;
l++;
}
for(int i=0;i<=lim;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<l);
NTT(a,rev,lim,1,mod,G);
NTT(b,rev,lim,1,mod,G);
for(int i=0;i<=lim;i++)
a[i]=a[i]*b[i]%mod;
NTT(a,rev,lim,-1,mod,G);
for(int i=1;i<n;i++)
g[i]=a[i-1]*quick_pow(i,mod-2)%mod;
g[0]=0;
}
```
:::
### 多项式 exp
设给定函数为 $h(x)$，则 $G(x,y)=\ln y-h(x)$。设 $f'(x)$ 为求出的模 $x^{\lceil \frac{n}{2} \rceil}$ 意义下的答案，我们有
$$\begin{aligned}
f(x)&\equiv f'(x)-\frac{\ln f'(x)-h(x)}{\frac{1}{f'(x)}} \\
&\equiv f'(x)(1-\ln x+h(x)) \pmod {x^n}
\end{aligned}$$
求 $\ln$ 加法乘法迭代，复杂度分析同求逆，时间复杂度 $O(n\log n)$。
:::info[参考代码]
```cpp
void expNTT(int f[],int g[],int n,int mod,int G)
{
if(n==1)
{
g[0]=1;
return;
}
expNTT(f,g,n+1>>1,mod,G);
lnNTT(g,tmpb,n,mod,G); // 1-ln(F_0(x))
tmpb[0]=(f[0]+1-tmpb[0]+mod)%mod;
for(int i=1;i<n;i++)
tmpb[i]=(f[i]-tmpb[i]+mod)%mod;
int l=-1,lim=1;
while(lim<=(n<<1))
{
lim<<=1;
l++;
}
for(int i=0;i<=lim;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<l);
NTT(tmpb,rev,lim,1,mod,G);
NTT(g,rev,lim,1,mod,G);
for(int i=0;i<=lim;i++)
g[i]=g[i]*tmpb[i]%mod;
NTT(g,rev,lim,-1,mod,G);
for(int i=n;i<=lim;i++)
g[i]=tmpb[i]=0;
}
```
:::
### 多项式快速幂
我们有
$$f(x)^k = \exp (k \ln f(x))$$
做完了，时间复杂度 $O(n\log n)$。
:::info[参考代码]
```cpp
INLINE void powNTT(int *f,int *g,int n,const string &K)
{
static int h[N];
int tmpk=0,tk=0,k=0;
for(auto ch:K)
{
k=((k<<3)+(k<<1)+(ch^48))%mod;
tk=((tk<<3)+(tk<<1)+(ch^48))%(mod-1);
tmpk=min(1000000,(tmpk<<3)+(tmpk<<1)+(ch^48));
}
int c=0;
while(f[c]==0&&c<n)
c++;
if((f[0]==0&&tmpk>=n)||c*k>=n)
{
for(int i=0;i<n;i++)
g[i]=0;
return;
}
int tmpn=n;
for(int i=c;i<n;i++)
f[i-c]=f[i];
for(int i=n-c+1;i<n;i++)
f[i]=0;
n-=c;
int iv=quick_pow(f[0],mod-2),v=quick_pow(f[0],tk);
for(int i=0;i<n;i++)
f[i]=1ll*f[i]*iv%mod;
lnNTT(f,g,n);
for(int i=0;i<n;i++)
g[i]=1ll*g[i]*k%mod;
expNTT(g,h,n);
for(int i=0;i<c*k;i++)
g[i]=0;
for(int i=c*k;i<tmpn;i++)
g[i]=1ll*h[i-c*k]*v%mod;
}
```
:::
### 多项式开方
一种思路是 $\sqrt {f(x)} = f(x) ^{1/2}$，然后按照快速幂的方法求解。
当然，也可以多项式牛顿迭代，设给定函数为 $h(x)$，则 $G(x,y)=y^2-h(x)$。设 $f'(x)$ 为求出的模 $x^{\lceil \frac{n}{2} \rceil}$ 意义下的答案，我们有
$$\begin{aligned}
f(x)&\equiv f'(x)-\frac{f'^2(x)-h(x)}{2 f'(x)} \\
&\equiv \frac{f'(x)+h(x)}{2f'(x)} \pmod {x^n}
\end{aligned}$$
时间复杂度 $O(n\log n)$。
:::info[参考代码]
```cpp
int quadrasique(int x)
{
mt19937 mt(time(0));
int s=mt()%(mod-1)+1;
while(quick_pow(w=(s*s%mod-x+mod)%mod,mod-1>>1)<=1)
s=mt()%(mod-1)+1;
fcomplex<int> tmp=fcomplex<int>(s,1);
tmp=quick_pow(tmp,mod+1>>1);
return min(tmp.x,mod-tmp.x);
}
void sqrtNTT(int f[],int g[],int n)
{
static int h[N];
if(n==1)
{
g[0]=quadrasique(f[0]);
assert(g[0]*g[0]%mod==f[0]);
return;
}
sqrtNTT(f,g,n+1>>1);
int l=-1,lim=1;
while(lim<=(n<<1))
{
lim<<=1;
l++;
}
for(int i=0;i<lim;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<l);
invNTT(g,h,rev,n);
NTT(g,rev,lim,1);
NTT(h,rev,lim,1);
for(int i=0;i<lim;i++)
g[i]=g[i]*g[i]%mod;
NTT(g,rev,lim,-1);
for(int i=n;i<lim;i++)
g[i]=0;
for(int i=0;i<n;i++)
g[i]=(g[i]+f[i])%mod;
NTT(g,rev,lim,1);
for(int i=0;i<lim;i++)
g[i]=g[i]*inv2%mod*h[i]%mod;
NTT(g,rev,lim,-1);
for(int i=n;i<lim;i++)
g[i]=0;
for(int i=0;i<lim;i++)
h[i]=0;
}
```
:::