P3803
【模板】多项式乘法（FFT）
零、前言
因为内容实在太多，不同部分字母表示的含义可能有所差别。
虽然已经极力精简了语言。但这篇文章仍然达到了 20KB 的可怕大小。
其中不免有疏漏之处，恳请来赏眼的巨佬能指出错误。
本来我是不想给文章分章节标号的。但后来发现真的写得太长了。。。
不知不觉写了 3 个多小时。
估计这篇文章会被一大堆人 D。初中还不会多项式的我真的是屑。
不要 D 我啊。
一、概论
FFT（快速傅里叶离散变换）是一种可以快速实现 DFT 和 IDFT 的方法。FFT 可以快速计算加法卷积（$O(n\log n)$）。
这个板子就是最简单的例子，我们所作的多项式乘法，就是加法卷积的一种。
这题我们分为两个过程，DFT 和 IDFT。
DFT（离散傅里叶变换）的作用就是将多项式的系数表达变为点值表达。
IDFT（离散傅里叶变换的逆变换）的作用就是将多项式的点值表达变为系数表达。
二、基本思路
至于我们为什么要转化为点值表达再转化为系数表达，是因为在点值表达下，我们可以更加快速地进行多项式乘法。
一般的多项式乘法是将系数卷积，这个时间复杂度显然是 $O(nm)$ 的。
一个显然的事实是 $n+1$ 个点可以唯一确定一个次数不超过 $n$ 的多项式。所以我们尝试用点值表达多项式。
比如 $F(x)$ 过 $(x_0,F(x_0))$，$(x_1,F(x_1))$，$\cdots$，$(x_{n+m},F(x_{n+m}))$。我们用这 $n+m+1$ 个点确定了 $F(x)$ 这个多项式。
同理可以用 $(x_0,G(x_0))$，$(x_1,G(x_1))$，$\cdots$，$(x_{n+m},G(x_{n+m}))$ 这 $n+m+1$ 个点来表示 $G(x)$ 这个多项式。
我们分别找到了 $x_0,x_1,\cdots,x_{n+m}$ 上 $F(x)$ 和 $G(x)$ 的值，显然我们可以很轻松的求出 $F(x)G(x)$ 在这 $n+m+1$ 个点上的值。
第一次听这个事实可能觉得有点奇怪，很有可能是 whk 做多了导致的。
这么说：$x_0$ 上 $F(x)G(x)$ 的值就是 $F(x_0)G(x_0)$，显然你知道 $F(x_0)$ 和 $G(x_0)$，然后你一乘就完了。
就是这样。所以很显然，在知道 $F(x)$ 和 $G(x)$ 的点值表达之后，可以在 $O(n+m)$ 的时间求出 $F(x)G(x)$ 的点值表达。
这样以来，我们的时间复杂度就取决于点值表达和系数表达来回转换的复杂度了。
首先我们知道，想要系数表达 $\rightarrow$ 点值表达，这个点的 $x$ 在哪里取完全取决于我们。
但是很容易知道，在系数表达 $\rightarrow$ 点值表达的过程中如果直接往这个 $n$ 次多项式里代入一般数值，并采取一般方法计算的话，复杂度仍然高达 $O(n^2)$。
并且在点值表达 $\rightarrow$ 系数表达这个过程中，需要搞出 $n+m+1$ 个线性方程组成线性方程组，然后高斯消元的复杂度是 $O(n^3)$ 的，用不起。
接下来我们需要解决点的就是如何加速这两个过程。
三、单位根
首先，我们要知道，直接代入一般数值是不行的。我们要代入一些更优（du）美（liu）的东西。
比如说，我们可以考虑复数。当然一般的复数还不行，我们代入的是单位根。
满足 $x^n-1=0$ 的 $x$ 被称作 $n$ 次单位根。需要注意的是，这里的 $x$ 可以是复数，所以情况会比一般的复杂一些。
这个图是我直接窃来的。
图中是一个复平面上的单位元。
我们设 $\omega_n$ 表示 $n$ 次单位根。在复平面上，它正好代表了 $\overrightarrow{OA}$ 这个向量。它的数值就是 $\exp(i\dfrac{2\pi}{n})=\cos{\dfrac{2\pi}{n}}+i\sin{\dfrac{2\pi}{n}}$。
因为复数相乘在复平面上表现为向量的模长相乘，辅角相加。所以这里如果我们给 $\omega_n$ 取 $k$ 次幂的话，$\omega_n^k$ 就对应了复平面上一个辅角占 $\dfrac{k}{n}$ 圆、长度为 1 的向量。
设 $n=2m$，我们需要用到的性质是：
$(\omega_n^k)^2=\omega_m^k$。考虑占多少份单位圆就是显然的了。
$(\omega_n^{m+k})^2=\omega_m^{m+k}=\omega_m^m\omega_m^k=\omega_m^k$。
$\omega_n^{m+k}=\omega_n^m\omega_n^k=-1\cdot\omega_n^k=-\omega_n^k$。
四、DFT 与 IDFT
了解了单位根，现在我们来看一下 DFT 和 IDFT 是什么。
接下来我们给出 DFT 的定义。
设 $\mathbf{f}$ 是长度为 $n$ 的数列，$\forall k\in[0,n)$，令：
$$\hat{f_k}=\sum_{i=0}^{n-1}f_i(\omega_n^k)^i$$
则称 $\mathbf{\hat{f}}$ 为 $\mathbf{f}$ 的 DFT（离散傅里叶变换）。
这里我们设 $\mathbf{f}$ 为多项式 $F(x)$ 的系数。那么 $\hat{f}_k$ 显然就是 $F(x)$ 在 $x=\omega_n^k$ 处的点值表示。
我们再给出 IDFT 的定义：
设 $\mathbf{\hat{f}}$ 是长度为 $n$ 的数列，$\forall k\in[0,n)$，令：
$$f_k=\dfrac{1}{n}\sum_{i=0}^{n-1}\hat{f_i}(\omega_n^{-k})^i$$
则称 $\mathbf{f}$ 是 $\mathbf{\hat{f}}$ 的 IDFT（离散傅里叶逆变换）。
我毫无根据地给出了这个 IDFT 的定义，但实际上定义的这两个式子是等价的。
下面我们给出关于这两个式子等价的推导。
我们将 DFT 的定义式代入 IDFT 的定义式的右边：
$$\dfrac{1}{n}\sum_{i=0}^{n-1}\hat{f_i}(\omega_n^{-k})^i=\dfrac{1}{n}\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}f_j(\omega_n^i)^j(\omega_n^{-k})^i$$
后面的是关于 $\omega_n$ 的次幂，我们可以把它们乘到一起：
$$\dfrac{1}{n}\sum_{i=0}^{n-1}\hat{f_i}(\omega_n^{-k})^i=\dfrac{1}{n}\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}f_j\omega_n^{(j-k)i}$$
交换一下求和顺序：
$$\dfrac{1}{n}\sum_{i=0}^{n-1}\hat{f_i}(\omega_n^{-k})^i=\dfrac{1}{n}\sum_{j=0}^{n-1}f_j\sum_{i=0}^{n-1}\omega_{n}^{(j-k)i}$$
$k$ 是一个已知的数，$j$ 是我们枚举的数。
分类讨论 $j$ 与 $k$ 的关系：
当 $j=k$ 的时候，有：
$$\sum_{i=0}^{n-1}\omega_{n}^{(j-k)i}=\sum_{i=0}^{n-1}1=n$$
当 $j\not=k$ 的时候，有：
$$\sum_{i=0}^{n-1}\omega_{n}^{(j-k)i}=\sum_{i=0}^{n-1}(\omega_n^{j-k})^i$$
这是个等比数列求和，也就是：
$$\sum_{i=0}^{n-1}\omega_{n}^{(j-k)i}=\omega_n^{j-k}\dfrac{1-(\omega_n^{j-k})^n}{1-\omega_n^{j-k}}$$
很显然 $(\omega_{n}^{j-k})^n=(\omega_n^n)^{j-k}=1$。所以：
$$\sum_{i=0}^{n-1}\omega_{n}^{(j-k)i}=\omega_n^{j-k}\dfrac{1-1}{1-\omega_n^{j-k}}=0$$
讨论完了，现在我们回顾上式：
$$\dfrac{1}{n}\sum_{i=0}^{n-1}\hat{f_i}(\omega_n^{-k})^i=\dfrac{1}{n}\sum_{j=0}^{n-1}f_j\sum_{i=0}^{n-1}\omega_{n}^{(j-k)i}$$
$\sum_{j=0}^{n-1}f_j\sum_{i=0}^{n-1}\omega_{n}^{(j-k)i}$ 这个式子的右边在 $j$ 枚举到 $k$ 的时候是 $nf_j=nf_k$，其余时候都是 $0$，那么它的值就是 $nk$。所以：
$$\dfrac{1}{n}\sum_{i=0}^{n-1}\hat{f_i}(\omega_n^{-k})^i=\dfrac{1}{n}nf_k=f_k$$
就是说，DFT 和 IDFT 的两个定义式是等价的。
五、Cooley-Tukey 算法
知道了 DFT 和 IDFT，我们仍不能高效地实现 DFT 和 IDFT（$O(n^2)$）。
接下来我们就要使用 Cooley-Tukey 算法优化。（具体关于这个名字可以看 这里）
Cooley-Tukey 算法的基本思想是分治。
设 $n=2m$。考虑将 $F(x)$ 的项按次数的奇偶分类。
$$F(x)=\sum_{i=0}^{n-1}f_ix^i=\sum_{i=0}^{m-1}f_{2i}x^{2i}+\sum_{i=0}^{m-1}f_{2i+1}x^{2i+1}=\sum_{i=0}^{m-1}f_{2i}x^{2i}+x\sum_{i=0}^{m-1}f_{2i+1}x^{2i}$$
这个结构很有意思。
我们定义 $F_l(x)$（可以理解为 $F_{even}(x)$，表示系数的下标为偶数）和 $F_r(x)$（$F_{odd}$，奇数）：
$$F_l(x)=\sum_{i=0}^{m-1}f_{2i}x^i$$
$$F_r(x)=\sum_{i=0}^{m-1}f_{2i+1}x^i$$
那么就会有：
$$F(x)=F_l(x^2)+xF_r(x^2)$$
这就是我们分治的基础。
这张图也是窃来的。
可以把每一个小方块理解为系数。连在一起的小方块就是这些系数组成的多项式，表示我们要对这个多项式进行 DFT。
而接下来，就是算法的核心：蝴蝶操作。
$\forall k\in[0,n)$，将 $\omega_n^k$ 和 $\omega_n^{m+k}$ 分别代入 $F(x)$。就有：
$$F(\omega_n^k)=F_l((\omega_n^k)^2)+\omega_n^kF_r((\omega_n^k)^2)=F_l(\omega_m^k)+\omega_n^kF_r(\omega_m^k)$$
$$F(\omega_n^{m+k})=F_l((\omega_n^{m+k})^2)+\omega_n^{m+k}F_r((\omega_n^{m+k})^2)=F_l(\omega_m^{k})-\omega_n^{k}F_r(\omega_m^{k})$$
这两个推导分别用到了单位根的性质 1 和 2、3。
现在，也就是说，如果我们知道了 $\forall k\in [0,m)$，$\omega_m^k$（共 $m$ 个）上 $F_l(x)$ 和 $F_r(x)$ 的取值，我们可以在 $O(n)$ 的时间得到 $F(x)$ 在 $\forall k\in[0,n)$，$\omega_n^k$（共 $n$ 个）上的取值。
所以根据主定理：
$$T(n)=2T(\dfrac{n}{2})+\Theta(n)=\Theta(n\log n)$$
也就是我们再 $O(n\log n)$ 下求出了 $F(x)$ 在 $\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}$ 共 $n$ 个点上的取值。
换句话说，我们成功找到了 $F(x)$ 的点值表达，完成了 DFT。
肯定有人会问 IDFT 怎么办。
实际上和 DFT 几乎一模一样。
我们根据 IDFT 的定义式：
$$f_k=\dfrac{1}{n}\sum_{i=0}^{n-1}\hat{f}_i(\omega_n^{-k})^i$$
把刚才 DFT 求得的点值（也就是式子中的 $\hat{f}_i$）当作我们刚才过程中的系数，并且把代入的 $\omega_n^k$ 换为 $\omega_n^{-k}$。
用刚才的 Cooley-Tukey 算法，我们就得到了所有的 $nf_k$。所以系数就是这些值乘上 $n^{-1}$ 即可。
六、递归 FFT 代码实现
我们手动把这个次数调成 2 的次幂，方便求解。因此这个数组的范围也要比题目要求的大一些。
复数自己写结构体，就那几行。
因为 DFT 和 IDFT 在递归过程中的差别就是代入单位根的次幂不同，所以我们直接写入一个函数。
最后输出系数要 f[i].x/n，因为 C++ 默认这种式子换成整型都是下取整的，所以我们要手动实现四舍五入，加个 0.5 即可。
代码（递归实现）：
1 |
|
七、非递归 FFT 实现与优化
显然现在洛谷的机子跑得比较快，这份代码已经可以通过了。
但是，递归的写法有较大的常数，并且码量稍长，我们来考虑一种非递归的写法优化。
这种写法需要用到位逆序置换。
我们设 $\operatorname{rev}(x)$ 表示 $x$ 翻转其二进制位所得到的的数。
再令：
$$f_i’=f_{\operatorname{rev}(i)}$$
则每次对 $\mathbf{f}$ 进行蝴蝶操作（分奇偶）变成了对 $\mathbf{f’}$ 上相邻两个序列的操作。
把 $\mathbf{f}$ 转化为 $\mathbf{f’}$ 的过程称为位逆序置换。
可以看下面这个图来理解这个东西。
还是窃来的。
上面的二进制数表示都是正常的。
你发现分奇偶实际上就是按二进制数从低到高，是 0 的分为一类，是 1 的分为另一类。
那如果把这个二进制数反过来，就是按二进制数从高到低分了。
从高到低分，等于按大小来分。
那现在问题是，如何利用这个来非递归实现 FFT。
枚举处理问题的层数。先将最底层的子问题全部解决，然后再枚举上一层，再解决。
还有个问题是 $\operatorname{rev}(i)$ 的求法。
其实可以用一个简单的递推解决。
代码（非递归实现）：
1 |
|
八、NTT（快速数论变换 Number Theorem Transform）
因为 FFT 需要用到单位根，不可避免的就要设计浮点数的运算，导致精度出现问题，常数大。
并且，一个事实已经被证明，那就是在复数域 $\mathbb{C}$ 中，只有单位根能实现 DFT 和 IDFT。
然而很多计数题都只需要在模意义下完成，需要一个更快的替代品。
某人说，要有 NTT，于是便有了 NTT。
考虑一个模意义下雨单位根类似的东西：原根。
定义 $g$ 为模 $p$ 的原根，当且仅当 $g^0,g^1,\cdots,g^{p-2}$ 在模 $p$ 意义下互不相同。
可以证明的是，模素数的原根总是存在的。（事实上，有原根的数只有 $2,4,2p^c,p^c$，这里 $p$ 是奇素数）
原根和本原单位根实际上非常类似。换句话说，在模 $p$ 意义下，$g$ 可以被看作一个 $p-1$ 次本原单位根。
设 $n$ 满足 $n\mid (p-1)$，令 $\omega_n=g^{\frac{p-1}{n}}$，那么就会有：
$$\omega_n^n=(g^{\frac{p-1}{n}})^n=g^{p-1}\equiv 1\pmod{p}$$
最后一步根据费马小定理。
并且 $\omega_n^0,\omega_n^1,\omega_n^2,\cdots,\omega_n^{n-1}$ 互不相同。
当然，求原根并不是一件容易的事。
所以大部分的原根靠背，因为一般 NTT 对模数要求很高。因为分治算法中 $n$ 一般是 2 的幂，为了满足上述 $n\mid (p-1)$，$p-1$ 要有 2 的次幂，且次数尽量要高。
比较常见的模数有：
$$998244353=7\cdot17\cdot2^{23}+1,g=3$$
$$1004535809=479*2^{21}+1,g=3$$
我从 Fuyuki 的主页里直接窃来一张模数表（表中 $p=r\cdot2^k+1$，原根为 $g$）：
1 | p r k g |
实际上没有必要记太多，够用就行了。
优化的话可以预处理单位根，再加上判断取模。
这题我们直接取一个较大的模数（比如 $998244353$），然后根据它这个题的数据范围，显然最后答案取模与不取模是一样的。
代码：
1 |
|