多项式，形式幂级数与生成函数学习笔记集合

ix35 · 2020-08-06 11:06:28 · 个人记录

注：从 Typora 复制到洛谷博客后部分公式炸裂，慢慢修吧。

介绍了我所了解的几乎所有有关这个板块的知识内容。

读完这篇文章并不能使你的多项式和计数水平达到一个很高的高度，也许你不能搞清楚为什么 FFT 的本质是循环卷积，幂级数三角函数的组合意义是什么，但至少可以确保你不问出“FFT 做完后为什么只保留实数部分”，“为什么多项式 \ln 之后还是多项式”，“多项式 \exp 这种理论上的东西到底能有什么实际用途”之类的问题。

〇、目录

一、快速傅里叶变换与数论变换

（一）复数与单位根

（二）多项式的表示

（三）快速傅里叶变换

（四）逆离散傅里叶变换

（五）数论变换

（六）任意模数意义下的变换

（七）分治 FFT

（八）高维 FFT

二、一些多项式技巧

（一）秦九韶算法

（二）拉格朗日插值法

（三）多项式除法和取模

（四）多项式快速求值与插值

三、形式幂级数基础

（一）基本运算

（二）泰勒展开

（三）牛顿迭代法

（四）形式幂级数的更多运算

四、生成函数基础

（一）概念

（二）组合意义

一、快速傅里叶变换与数论变换

这一部分将会介绍：复数的基本概念，快速傅里叶变换（FFT），数论变换（NTT）及其部分衍生内容。

（讲复数的基本概念是为了引入单位根，而和原根有关的知识大都和单位根部分的证明相似，且在数论中原根较为熟悉，因此不再单独讲解原根）

（一）复数与单位根

复数，即可以表示成 z=a+bi\ \ (a,b\in \R) 形式的数，其中 a 称为 z 的实部，b 称为 z 的虚部，i 为虚数单位 \sqrt{-1}。

复数与平面直角坐标系中的二维向量一一对应，我们将 z=a+bi 对应向量 (a,b)，于是平面上的每一个点都代表了唯一的复数，此平面称为复平面，横轴称为实轴，纵轴称为虚轴。

对于 z=a+bi，如果 b=0，则 z 为实数，否则 z 为虚数。如果 a=0 且 b\ne 0，则称 z 为纯虚数。

复数的四则运算：

对于 z_1=a+bi,\ \ z_2=c+di，定义它们的四则运算如下：

1. 加法：z_1+z_2=(a+c)+(b+d)i，即实部和虚部分别相加，满足交换律和结合律。
2. 减法：z_1-z_2=(a-c)+(b-d)i，即实部和虚部分别相减。
3. 乘法：z_1\times z_2=(ac-bd)+(ad+bc)i，可用乘法分配率验证，满足交换律、结合律和对加法的分配率。
4. 除法：\dfrac{z_1}{z_2}=\dfrac{ac+bd}{c^2+d^2}+\dfrac{bc-ad}{c^2+d^2}i，可通过分母实数化得到，这里 c,d 不同时为 0。

对于复数 z=a+bi，定义 \bar z=a-bi 与 z 互为共轭复数，满足 z\bar z\in \R。

对于复数 z=a+bi，定义 r=\sqrt {a^2+b^2} 为其模长。

我们可以用类似极坐标的方法来表示一个复数：z=r(\cos \theta+\sin \theta i)，其中 r 为其模长，\theta 为其对应向量与实轴正半轴的夹角，\theta 称为 z 的辐角。

考虑用极式将两复数相乘：

利用三角函数的和角公式，我们可以得到：

于是我们可以总结出：两复数相乘，模长相乘，辐角相加。

接下来考虑方程 x^n=1，它有 n 个不同复数解。

例如，n=3 时，三个解为 x_1=1,\ \ x_{2,3}=\dfrac{1}{2}\pm\dfrac{\sqrt 3}{2}i。

将三个解作在复平面上，发现恰好将单位圆三等分。是否对于所有的 n 都有类似性质呢？

设 x=r(\cos \theta+\sin\theta i) 是方程的解，则 x^n=r^n(\cos(n\theta)+\sin(n\theta)i)=1，这说明：

其中，\theta>0 且最小的一个解是 \theta=\dfrac{2\pi}{n}，我们就将这个根 x=\cos(\dfrac{2\pi}{n})+\sin(\dfrac{2\pi}{n})i 称为 \omega_n，其他的根都可以表示为 x=\cos(\dfrac{2k\pi}{n})+\sin(\dfrac{2k\pi}{n})i=\omega_n^k，方程的这些根 x 都称为 n 次单位根。

我们还可以将单位根表示成 \omega_n^k=e^{\frac{2\pi ik}{n}} 的形式。

（二）多项式的表示

注意：我们讨论的多项式是广义的，即包含了狭义的多项式与单项式。

我们讲了上面这些，最终希望归到离散傅里叶变换，而在 OI 中这部分的研究主要是在多项式这一对象上的。

多项式的常见表示有两种：系数表示和点值表示。

1. 系数表示：即常见的 f(x)=\sum\limits_{i=0}^{n} a_ix^i 的表示方法，n 称为多项式的次数。
2. 点值表示：用 n+1 对点值表示一个 n 次多项式的方法，即给出 (x_0,f(x_0)),\ldots,(x_n,f(x_n))。

其中系数表示便于直观地分析多项式的若干性质，以及对多项式进行一些常规的操作。

但是，一旦设及到乘法（称为卷积），系数表示的多项式乘法使用暴力的复杂度为 O(n^2)，比较慢。

实际上有一种容易实现的多项式乘法优化，叫做 Karatsuba 乘法，基于每次将一个 n 次的多项式乘法转化成 3 个 \dfrac{n}{2} 次的多项式乘法，根据 Master 定理可分析得到其复杂度为 O(n^{\log_2 3})\approx O(n^{1.585})，但效率依然不够高。

如果我们用点值表示进行卷积，(x_0,f(x_0)),\ldots,(x_n,f(x_n)) 与 (x_0,g(x_0)),\ldots,(x_n,g(x_n)) 的卷积就是 (x_0,f(x_0)g(x_0)),\ldots,(x_n,f(x_n)g(x_n))，这个过程是 O(n) 的，且如果 n 等于 f(x) 和 g(x) 的次数之和，那么这样的点值表示可以唯一确定一个 n 次多项式。

这就让我们想到：如果我们能实现一个点值表示和系数表示之间的转换，那么两个系数表示多项式的卷积就可以先转换成点值做卷积以后再转换回来。

（三）快速傅里叶变换

在下面的讨论中，我们假设要操作的多项式为 n-1 次，且 n=2^p\ (p\in \Z^+)，如果次数不足只需补至 2 的整数幂次即可。

我们要把多项式 f(x)=\sum\limits_{i=0}^{n-1}a_i\times x^i 转化成其点值表示。

点值表示中的 x_i 是可以任意选取的，我们不妨拓宽数域，选择 \omega_n^k\ (k\in [0,n)) 作为所有的 x_i。

之所以选择 \omega_n^k 是因为它有许多好的性质，下面列举一些：

性质 1：\omega_{2n}^{2k}=\omega_n^k。

证明：\dfrac{2\times 2k\times \pi}{2n}=\dfrac{2k\pi}{n}，带入之前写过的式子即得证。

性质 2：\omega_n^k=-\omega_n^{k+\frac{n}{2}}。

证明：\omega_{n}^{k+\frac{n}{2}}=\omega_n^k\times \omega_{n}^{\frac{n}{2}}，而 \omega_{n}^{\frac{n}{2}} 的辐角为 \dfrac{2\times \frac{n}{2}\times \pi}{n}=\pi，所以 \omega_{n}^{\frac{n}{2}}=-1，即得证。

性质 3：\omega_{n}^n=1。

证明：由定义可得证。

接下来考虑如何利用单位根快速实现多项式不同表示方法的转换。

首先，对于 0\leq k<\dfrac{n}{2}，我们分别带入 f(\omega_n^k) 和 f(\omega_n^{k+\frac{n}{2}})：

根据性质 2，我们知道：当 i\equiv 0\bmod 2 时，(\omega_n^{k})^i=(\omega_n^{k+\frac{n}{2}})^i，否则 (\omega_n^{k})^i=-(\omega_n^{k+\frac{n}{2}})^i，于是我们可以用 \omega_n^k 来表示 f(\omega_n^{k+\frac{n}{2}})：

把 i\equiv 0\bmod 2 和 i\equiv 1\bmod 2 的两个部分分开，写成两个部分的和：

对于 (\omega_n^k)^{2i}，可以写成 (\omega_n^{2k})^i，从而根据性质 1 写为 (\omega_{\frac{n}{2}}^k)^i；而对于 (\omega_n^k)^{2i+1}，我们提取一个 \omega_n^k 后用类似方法处理，所以可以得到：

为了追求相似的形式，我们对 f(\omega_n^k) 进行类似的变换：

结合以上两式，发现两个括号内的和式是完全一样的，所以我们可以设两个新的多项式：

于是，f(\omega_n^k) 和 f(\omega_n^{k+\frac{n}{2}}) 的关系可以显式地表达如下：

而我们发现，A(\omega_{\frac{n}{2}}^k),\ B(\omega_{\frac{n}{2}}^k) 这两个 \dfrac{n}{2} 次多项式的点值，本质和求 f(\omega_n^k) 没有区别，只是一个规模一半的子问题，所以我们可以递归求解，递归到 n=1 时直接有 f(x)=a_0。

下面来分析这个做法的时间复杂度。

每一个规模为 n 的问题（求所有 f(\omega_n^k)）可以拆分成两个规模为 \dfrac{n}{2} 的子问题（求所有 A(\omega_{\frac{n}{2}}^k) 和 B(\omega_{\frac{n}{2}}^k))，然后用 O(n) 的代价合并（每个 f(\omega_n^k) 可以 O(1) 通过两个 A,B 的点值得到），所以整个过程的复杂度为 O(n\log n)（分析过程非常经典，与归并排序的复杂度证明类似）。

上面这个过程叫做离散傅里叶变换（DFT），而这种基于分治的优化即为快速傅里叶变换（FFT）。

DFT 与一般傅里叶变换的关系：一般傅里叶变换是积分，DFT 是其离散的版本，积分变为了求和，这个可以稍作了解，毕竟 OI 中基本用不到连续版本的傅里叶变换。

当然，这个递归版本写起来有些麻烦，而且常数较大，所以我们通常使用的是一种非递归版本的 FFT。

观察 FFT 过程中的递归树（图片来源于网络，以 n=16 为例）：

其中 \text{Step}\ i 即为第 i 层的递归，写在同一块中的下标即表示同一个多项式的系数。

做完递归以后，我们把下标从左往右写出来（也就是把图中 \text{Step}\ 4 中的数从左到右从上到下写一遍，特别地，最后两个块在图上的位置反了），与其位置对比，可以看出：最终一个位置上的值和这个位置的二进制表示正好是反转的。

简单说明一下：因为每一步只有处于当前多项式偶数次的系数会走到左边，奇数次的系数会走到右边，所以如果 x 的第 i 位为 1，那么在第 \log_2 n-i+1 步中它就会走到右边，否则就会走到左边，而它最终的位置就是由每一步走到左边还是右边来决定的，于是不难得到上面的结论。

于是，我们可以将本来的从上往下递归，改为从下往上递推：

设 r_i 表示 i 的二进制表示反转后的结果，那么首先我们将 a_i 与 a_{r_i} 互换，先得到了最后的状态。

然后，我们第 i 轮就将连着的两个大小为 2^{i-1} 的块合并成一个大小为 2^i 的块（如之前的图中所示，从 \text{Step}\ 3 到 \text{Step}\ 2 就是把左边两个大小为 4 的块合并成一个大小为 8 的块，右边两个大小为 4 的块也合并成一个大小为 8 的块），这一步和递归 FFT 的实现是一样的。

这样常数略小，而且很容易写，复杂度也非常便于分析（总共 \log_2 n 轮，每轮合并总共 O(n)）。

注意在实现时，我们可以直接令当前某一块中的第 i 个下标处存储 \omega_{2^l}^i 的点值（2^l 是当前块长），这样我们将两个函数的运算转化成了两个变量的运算，更加直观。

（四）逆离散傅里叶变换

在下面的讨论中依然认为 n 是 2 的正整数幂次。

考虑如何将多项式的点值表示转换成系数表示。

上面关于单位根的性质讲到了三条，这里再补充一条：

性质 4：\sum\limits_{i=0}^{n-1}\omega_n^{ik}=\begin{cases}n&(k=0) \\0&(k\ne 0) \end{cases}。

证明：当 k=0 时，\sum\limits_{i=0}^{n-1}\omega_n^{ik}=\sum\limits_{i=0}^{n-1}1=n；而当 k\ne 0 时，由等比数列求和公式有 \sum\limits_{i=0}^{n-1}\omega_n^{ik}=\dfrac{1-\omega_n^{kn}}{1-\omega_n^k}，由于 k\ne 0 时分子为 0，所以结果为 0。

我们设 a_i 是 n-1 次多项式 f(x) 的 x^i 项系数，而 A_i 是 \omega_n^i 的点值，那么用矩阵的形式来表达带入点值的运算如下：

记前面这个方阵为 D，我们现在要倒过来通过 A_i 求 a_i，需要找到 D 的逆。

我们先给出结论，再做证明：

我们设左侧的这个方阵为 E，下面尝试证明 ED=nI。

对于 ED 的第 i 行第 j 列元素（设为 M_{ij}），可以得到（方便起见，行列从 0 到 n-1 编号）：

再根据单位根性质 4，我们得到：当 j-i=0 即 i=j 时，M_{ij}=n，否则 M_{ij}=0，这正是单位矩阵 I 的 n 倍，于是得证。

更有趣的是，我们发现：单位根的前三条性质中，将 \omega_n^k 换成 \omega_n^{-k} 后依然全部成立（可以自行验证），所以逆离散傅里叶变换（IDFT）就简单了，只需要在 DFT 过程中将所有的 \omega_n^k 变成 \omega_n^{-k} 即可。

由于 ED=nI，所以最后还需要再将结果除以 n 得到真实的各项系数。

利用同样的优化，即可得到逆快速傅里叶变换（IFFT）。

至此，我们已经能够实现多项式乘法了，只需要先将两个多项式 f(x),g(x) 求 FFT，用点值相乘再求 IFFT，即可得到 h=f*g 的各项系数了，时间复杂度为 O(n\log n)。

（五）数论变换

在下面的讨论中依然认为 n 是 2 的正整数幂次。

利用快速傅里叶变换求多项式乘法，有一些明显的缺点：

1. 精度较低，如果系数较大可能丢失精度产生误差（即使给出的系数都是整数）；
2. 无法取模，只支持较小规模的数据（这一点在讨论幂级数的时候更加突出，当我们遇到泰勒展开等运算时需要大量用到阶乘，FFT 难以正常处理）。
3. 复数运算慢，算法常数大。

基于此，我们需要一种不依赖与复数域进行操作的变换，而它就是数论变换（NTT）。

原本我们一直在复数域中讨论，现在我们换一个角度，从数论中模意义下的数域出发讨论。

进行数论变换时，我们得到的所有值都是在 \bmod 一个大素数的意义下的，这个大素数 p 通常需要满足 p=b\times 2^k+1，其中 2^k\ge n，稍后将解释这为什么是必要的。

然后，我们将原本 FFT 中的单位根替换为模 p 的原根（这里不打算像介绍单位根那样系统介绍一遍原根的相关知识了），设 g 是模 p 的某一个原根，设 \omega_n=g^{\frac{p-1}{n}}，其余和 FFT 完全一致。

我们知道 FFT 以及 IFFT 的核心在于单位根的四个性质，下面来论证一下这四个性质在我们目前的定义下（在模 p 的数域内）依然成立：

性质 1：\omega_{2n}^{2k}=\omega_n^k。

证明：\omega_{2n}^{2k}=g^{\frac{p-1}{2n}\times 2k}=g^{\frac{p-1}{n}\times k}=\omega_n^k，即得证。

性质 2：\omega_n^k=-\omega_n^{k+\frac{n}{2}}。

证明：\omega_{n}^{k+\frac{n}{2}}=\omega_n^k\times \omega_{n}^{\frac{n}{2}}，而 \omega_n^{\frac{n}{2}}=g^{\frac{p-1}{n}\times \frac{n}{2}}=g^{\frac{p-1}{2}}，由于 g^{p-1}=1 且 g 是模 p 的原根，所以不存在 q<p-1 使得 g^{q}=1，所以 g^{\frac{p-1}{2}}=-1，即得证。

性质 3：\omega_{n}^n=1。

证明：\omega_n^n=g^{\frac{p-1}{n}\times n}=g^{p-1}=1，即得证。

性质 4：\sum\limits_{i=0}^{n-1}\omega_n^{ik}=\begin{cases}n&(k=0) \\0&(k\ne 0) \end{cases}。

证明：当 k=0 时，\sum\limits_{i=0}^{n-1}\omega_n^{ik}=\sum\limits_{i=0}^{n-1}1=n；而当 k\ne 0 时，由等比数列求和公式有 \sum\limits_{i=0}^{n-1}\omega_n^{ik}=\dfrac{1-\omega_n^{kn}}{1-\omega_n^k}，由于 k\ne 0 时分子为 0，所以结果为 0（和单位根部分的证明完全相同）。

所以这里我们定义的 \omega_n 符合了之前用到的所有性质，用它去做 FFT 和 IFFT，能支持取模，不用写复杂的复数运算，而且常数小，相比 FFT 确实有许多优势，时间复杂度是一样的 O(n\log n)。

当然也有不方便用 NTT 而可以用 FFT 做的题，例如题目中出现的系数不是整数时，就难以用 NTT 实现。

因为我们 \omega_n 的定义中用到了 g^{\frac{p-1}{n}}，所以 n\mid p-1，而 n 是一个 2 的正整数幂次，所以这就要求 p-1 是一个形如 b\times 2^k 的数了，其中 2^k 要充分大（不小于最大的可能的 n），这样的素数 p 中，最常见的就是 998244353，另一些如 1004535809，469762049 也可以作为 NTT 模数，这三个数的最小原根都是 3，比较方便记忆。

当然，也不是所有题目的模数都是这些形如 b\times 2^k+1 的数，遇到了不同的例子，我们就要用更多方法解决。

（六）任意模数意义下的变换

当多项式乘法的模数不再是常见的 b\times 2^k+1 型数时，只做一次 NTT 便不可行了，这里介绍两种常见的解决方法，实际效率都比较慢，也有一些较快的方法，往往实现较为复杂，而且我也不了解，读者可以查阅一些论文。

1. 三模数 NTT：选取 3 个可做 NTT 的模数（比如上面举出的三个例子），设为 P_1,P_2,P_3，分别进行 NTT，这样我们对于乘积的每一项系数都得到了 3 个模意义下的值 x_1,x_2,x_3，设真实的系数为 x，那么我们有：

   \begin{aligned} x\equiv x_1\bmod P_1 \\ x\equiv x_2\bmod P_2 \\ x\equiv x_3\bmod P_3 \end{aligned}

   我们要求的是 x\bmod p（p 是给定的模数），可以发现 x 的上界是 p^2n（每一项都乘满），当 p 在 int 存储范围内时，这个值大约为 10^{23} 左右（这也是为什么选择三个模数做 NTT 就足够的原因），我们只需用中国剩余定理合并三个同余式即可得到 x\bmod p（可以实现得比较精细，并不需要高精度，但可能需要快速乘）。

2. 分部 FFT：我们发现，只要按照 FFT 做一次，最后再取模就可以了，但是事实上这样不仅数据范围过大，精度也是很差的，所以我们可以通过将系数减小来避免精度损失。

   设 d=\lceil\sqrt p\rceil，对于原多项式的某个系数 a_i，写成 a_i=kb_i+c_i 的形式，其中 0\leq c_i<k，设 F(x)=kF_1(x)+F_2(x),\ G(x)=kG_1(x)+G_2(x)，那么：

   F(x)\times G(x)=k^2F_1(x)G_1(x)+k(F_1(x)G_2(x)+F_2(x)G_1(x))+F_2(x)G_2(x)

   只需要分别计算四个系数较小的多项式乘积，再合并起来即可。

两种方案时间复杂度都是 O(n\log n)，但是常数都较大。

（七）分治 FFT

有时我们需要处理的卷积形式比较特殊，例如：

给定一个序列 g_{1,\ldots n} 和 f_0=1，且满足 f_i=\sum\limits_{j=1}^i f_{i-j}g_j，求序列 f_{1,\ldots,n}。

此时，我们发现 f 自己也会影响自身后续项的取值，暂且不论求逆等其他做法，单纯使用 FFT 是难以解决这个问题的。

此时我们可以考虑将 FFT 与 CDQ 分治的思想结合。

我们将序列 f_{1,\ldots,n} 分治，设 mid=\dfrac{n+1}{2}，先递归计算出 f_{1,\ldots,mid}，然后考虑 f_{1,\ldots,mid} 对 f_{mid+1,\ldots,n} 的贡献，计算完后再递归求解 f_{mid+1,\ldots,n} 内部的贡献即可。

我们新设一个序列 h_i=\begin{cases}f_i & (1\leq i\leq mid) \\ 0 & (mid<i\leq n) \end{cases}，考虑计算序列 h 和序列 g 的卷积，得到的就是前半部分对后半部分的贡献了。

每一层计算贡献的复杂度是 O(n\log n)，因此总的复杂度是 O(n\log^2 n)。

（八）高维 FFT

由于这一部分我自己也没写过，所以只能简单讲讲。

我们之前的讨论都是基于一元多项式的，那么对于二元多项式甚至更多元的多项式，如何进行点值表示和系数表示的转换呢？

以二元多项式为例，设最高次数为 x^{n-1}y^{m-1}，我们假设 n,m 都是 2 的正整数幂次，我们将系数排成一个 n\times m 的矩阵（行列从 0 开始编号），那么我们需要 n\times m 个点值，对于第 i 行第 j 列将要带入的点值，令它就是 x=\omega_n^i,\ y=\omega_m^j。

先考虑每一行，每一行中 x 的次数都是一样的，所以我们主要关心的是 y，如果我们直接把这一行的 m 个数当做一个多项式的系数做一次 FFT，得到的就是一个关于 x^i 的一元多项式的各项系数（i 是这一行的行号）。

接下来，我们再考虑每一列，每一列中 y 的次数都是一样的，而我们将系数从上往下看，它们对应的次数为 x^0,\ldots,x^{n-1}，所以我们直接对这一列做一次 FFT，就得到了这一列最终的点值表示结果。

由于对每行每列都做了 FFT，所以时间复杂度为 O(nm(\log n+\log m))。

对于更高维的情况也是类似的，只要对每一个维度上都做一维的 FFT，最终就能得到完整的 FFT 结果。设每一维的最高次数分别为 n_1,\ldots,n_k，则复杂度就是 O(\Big(\prod\limits_{i=1}^k n_i\Big)\times \Big(\sum\limits_{i=1}^k \log n_i\Big))。

至此，我们初步了解了快速傅里叶变换及其逆变换，数论变换及其逆变换，用它们处理多项式乘法的方法以及一些常见的相关扩展内容，此后的部分中，不再特别区分 FFT 和 NTT，它们都是计算卷积的工具。

卷积的实质是形如 c_x=\sum\limits_{i+j=x}a_ib_j 这样的运算，因此 FFT 并不局限于计算多项式乘法，还可以扩展到其他带有这样卷积形式的式子。

下面试举几例，说明 FFT 在处理卷积问题上的应用（然而与后文关系不大）。

例题 1： [ZJOI 2014] 力

题目大意：

给定 n 个实数 q_1,\ldots,q_n，求出所有 E_i\ (1\leq i\leq n)，其中：

（实际上是当质量为 q_i 的质点位于点 i 时，点 i 的加速度，分子为点 i 受到的万有引力，这应该也是题名由来。）