多项式四则运算

前言

在数学中，由若干个单项式相加组成的代数式叫做多项式（若有减法：减一个数等于加上它的相反数）。多项式中的每个单项式叫做多项式的项，这些单项式中的最高项次数，就是这个多项式的次数。其中多项式中不含字母的项叫做常数项。
在数学中，多项式（polynomial）是指由变量、系数以及它们之间的加、减、乘、幂运算（非负整数次方）得到的表达式。对于比较广义的定义，一个或零个单项式的和也算多项式。按这个定义，多项式就是整式。实际上，还没有一个只对狭义多项式起作用，对单项式不起作用的定理。$0$ 作为多项式时，次数定义为负无穷大（或 $0$）。单项式和多项式统称为整式。
（摘自：百度百科）

多项式的定义初中就有，没必要多说吧。本文将会积累所有学了的有关多项式的知识，不定期更新内容。

本博客全文所有内容均不属于提高组考试范围，你要是硬说多项式加法可能有……那我也不能反驳。

多项式求和，单项式乘多项式

合并同类项，没必要多说吧。多项式求差同理。

单项式乘多项式就拆个括号。

多项式求积

假设对于两个多项式 $A,B$，有：

$A=\sum_{i=0}^{n-1} a_ix^i \\ B=\sum_{i=0}^{m-1} b_ix^i$

考虑如何求出两个多项式的积 $C=A\ast B$。

系数表示法，点值表示法

考虑上面的这两个式子，如果我们暴力地求解，不难得到最后的结果为：

$C=\sum_{i=0}^{n+m-2} \sum_{j=0}^i a_j b_{i-j}x^i$

暴力求解是 $\operatorname{O}(n^2)$ 的，复杂度不太能接受的亚子。我们考虑优化。

我们知道这种表示多项式的方法是最常见的，叫做系数表示法，众所不周知多项式的表示还可以用一种方法叫做点值表示法。根据初中数学知识，两点确定一个一次函数，三点确定一个二次函数……于是伟大的数学家们证明了一个定理：$n$ 个互异的点确定一个 $n-1$ 次函数。我们把 $A,B$ 用点值表示法表示就是：

$A=\{ (x_1,\operatorname{f}(x_1)),(x_2,\operatorname{f}(x_2)),(x_3,\operatorname{f}(x_3)),\cdots,(x_n,\operatorname{f}(x_n)) \} \\ B=\{ (x_1,\operatorname{g}(x_1)),(x_2,\operatorname{g}(x_2)),(x_3,\operatorname{g}(x_3)),\cdots,(x_n,\operatorname{g}(x_n)) \}$

然后考虑 $C$ 如何表达。我们举个例子，假设一次多项式 $A=\{(0,-1),(2,5)\}$，一次多项式 $B=\{(0,3),(2,-1)\}$。不难得出系数表示法 $A=3x-1,B=-x+3$，根据系数表示法，得出它们的乘积 $C=-3x^2+10x-3$（假设记为函数 $\operatorname{h}$）。画下来：

我们发现一个现象：对于 $\forall x\in R$，$\operatorname{f}(x)\operatorname{g}(x)=\operatorname{h}(x)$。因此用点值表示法运算两个多项式的乘积的方法是：（后续 $n$ 均表示为两个多项式项数的较大值，且不足 $n$ 项的那个多项式，少的项全部视为用 $0x^k$ 补充成 $n$ 项）

$C=\{ (x_i,\operatorname{f}(x_i)\operatorname{g}(x_i))\mid i\in N^+,i\in[1,n] \}$

时间复杂度是 $\operatorname{O}(n)$。也就是说如果、我们要求两个用系数表示法表示的多项式 $A,B$ 的积，可以把它们转化成点值表示法求积，然后转化回系数表示法。

DFT

复数

形如 $a+bi$ 的数被称为复数，其中 $a$ 为实部，$b$ 为虚部，$i$ 为虚数单位（虚数单位不属于虚部）。当 $b=0$ 时，这个数就是实数，当 $a=0$ 时，它是个纯虚数。任何复数系多项式在复数部一定有根。复数的集合记作 $C$。

初中里面 $\sqrt{-1}$ 无意义，在复数域里，我们记作 $\sqrt{-1}=i$，$i$ 就是前面提到的虚数单位，再比如，$\sqrt{-9}=3i,\sqrt{-5}=\sqrt{5}i$。

将复数的实部和虚部分别用 $x$ 和 $y$ 表示，就得到了一个复平面直角坐标系，比如：

其中点 $A$ 是原点，表示 $(0,0)$，点 $B$ 为 $(1,2i)$，代表复数 $1+2i$。定义复数的模表示向量 $u$ 的大小，即 $\sqrt{a^2+b^2}$。在这里为 $\sqrt{2^2+1^2}=\sqrt{5}$。记复数的辐角（极角）为此处的 $\alpha$。复数的模和辐角构成极坐标，如 $B(\sqrt{5},\alpha)$。

复数满足四则运算，既然复数可以用向量 $u$ 表示，也就说明了复数的加法满足平行四边形法则。复数 $z=a+b_i$ 的 共轭复数 为 $\overline{z}=a-bi$。至于乘法，也就是把括号拆开算了，通式如下：

$(a+bi)\ast(c+di)=(ac-bd)+(ad+bc)i$

在这里这对 DFT 没有什么帮助，但是如果我们举个例子画在坐标系上，记复数 $z_1=2+i,z_2=-1-i$，套式子得到它们的乘积 $z_3=-1-3i$，画下来：

为了避免太乱就没有把三个辐角都标出来了，其中 $u$ 的辐角是 $\alpha$，$v$ 的辐角是优角 $\angle BAC+\alpha$，$w$ 的辐角是优角 $\angle BAD +\alpha$。我们发现 $\alpha=\beta$，同时 $AB\cdot AC=AD$。所以得出另一个结论：

两个复数相乘，模长相乘，辐角相加。

离散傅里叶变换

其实就是 DFT 的全称，考虑如何把一个系数表达式转化为点值表达式。

从这里开始，默认 $n$ 是 $2$ 的正整数次幂，如果不是的话，就把它补成高次项系数都是 $0$ 的形式。

显然可以随便选几个不同的 $x$，然后代入到多项式里面，求出对应的 $\operatorname{f}(x)$，就可以求出点值表达式了，因为要选出 $n$ 个点，多项式是 $n-1$ 次，有 $n$ 项，复杂度是 $n^3$ 级别的。用快速幂的话把复杂度优化为 $\operatorname{O}(n^2\log n)$，预处理 $x^1,x^2,x^3,\cdots,x^n$ 复杂度为 $\operatorname{O}(n^2)$。

傅里叶觉得不行，随便选 $x$ 太随意了，生活要有仪式感。虽然 DFT 复杂度也是 $n^2$，但是他觉得我要是选一些点使得不要算那么多次乘方运算，好歹常数也小点。如果选择出来的一些 $x$，使得 $\exist p,x^p=1$，那么所有的乘方运算常数都可以大幅度地下降。显然 $\pm 1,\pm i$ 都可以满足，但是只有四个 $x$ 太少了，$n$ 大于 $4$ 就没掉了。傅里叶发现了一个性质：如果在复平面直角坐标系上，原点为圆心，$1$ 为半径画个圆，圆上的点都满足。

这样就只需要在圆上面任取 $n$ 个点就可以了。傅里叶把这个圆 $n$ 等分了，假设 $n=8$，给一个图片方便后续的理解：

从 $(1,0)$ 开始，$(1,0)$ 为 $X_0$，然后逆时针标号直到 $X_{n-1}$。我们记 $\omega_n^k$ 表示 $X_k$ 对应的复数的值。根据复数相乘时，模长相乘辐角相加，不难得到 $(\omega_n^1)^k=\omega_n^k$。定义 $\omega_n^1$ 为 $n$ 次单位根。

那么把 $\omega_n^0,\omega_n^2,\omega_n^3,\cdots,\omega_n^{n-1}$ 当作 $x_0,x_1,x_2,\cdots,x_{n-1}$ 带进多项式。因为 $(\omega_n^1)^k=\omega_n^k$，$\omega_n^n=\omega_n^0=1$，就可以转点值了。

傅里叶觉得非常满足，很有仪式感地完成了系数表达转点值表达。总体复杂度 $\operatorname{O}(n^2)$。

如何把点值转化为系数表达式呢？拉格朗日插值 $n^2$ 秒掉了。

FFT

快速傅里叶变换

子标题就是 FFT 的全称。

傅里叶又觉得不行，$n^2$ 复杂度太拉了，于是有了 $n\log n$ 级别的优秀的 FFT。我们发现 $\omega_n^k$ 有以下几个优秀的性质：

• $(\omega_n^p)^k=\omega_n^{pk}$
• $\omega_n^n=\omega_n^0=1+0i=1$
• $\omega_n^k=\omega_{2n}^{2k}$
• $\omega_n^{k+\frac{n}{2}}=-\omega_n^k$
• $\omega_n^k\omega_n^p=\omega_n^{k+p}$

第一个：相当于旋转了 $k$ 次辐角；第二个：相当于旋转一周回来了；第三个：相当于 $n$ 等分变成 $2n$ 等分；第四个：相当于转了 $180^\circ$；第五个：辐角直接加起来。

重新考虑原系数表达式的形式，就分治的思想，按照指数的奇偶性分成两个部分，记原多项式为 $\operatorname{A}(x)$：

$\begin{aligned} \operatorname{A}(x) &= \sum_{i=0}^{n-1} a_ix^i \\ &= (a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+\cdots+a_{n-1}x^{n-1}) \\ &= (a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots+a_{n-1}x^{n-2}) \end{aligned}$

感觉两个部分长得差不多，考虑开两个多项式 $\operatorname{A_1}(x),\operatorname{A_2}(x)$：

$\operatorname{A_1}(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n}{2}-1} \\ \operatorname{A_2}(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n}{2}-1} \\ \therefore \operatorname{A}(x)=\operatorname{A_1}(x^2)+x\operatorname{A_2}(x^2)$

对于 $k<\frac{n}{2}$，把 $x=\omega_n^k$ 代入原式：

$\begin{aligned} \operatorname{A_1}(x^2)+x\operatorname{A_2}(x^2) &= \operatorname{A_1}[(\omega_n^k)^2]+\omega_n^k\operatorname{A_2}[(\omega_n^k)^2] \\ &= \operatorname{A_1}(\omega_n^{2k})+\omega_n^k\operatorname{A_2}(\omega_n^{2k}) \\ &= \operatorname{A_1}(\omega_{\frac{n}{2}}^k)+\omega_n^k\operatorname{A_2}(\omega_{\frac{n}{2}}^k) \end{aligned}$

考虑另一半（也就是 $k+\frac{n}{2}$），把 $x=\omega_n^{k+\frac{n}{2}}$ 代入原式：

$\begin{aligned} \operatorname{A_1}(x^2)+x\operatorname{A_2}(x^2) &= \operatorname{A_1}[(\omega_n^{k+\frac{n}{2}})^2]+\omega_n^{k+\frac{n}{2}}\operatorname{A_2}[(\omega_n^{k+\frac{n}{2}})^2] \\ &= \operatorname{A_1}[(-\omega_n^k)^2]-\omega_n^k\operatorname{A_2}[(-\omega_n^k)^2] \\ &= \operatorname{A_1}(\omega_n^{2k})-\omega_n^k\operatorname{A_2}(\omega_n^{2k}) \\ &= \operatorname{A_1}(\omega_{\frac{n}{2}}^k)-\omega_n^k\operatorname{A_2}(\omega_{\frac{n}{2}}^k) \end{aligned}$

和上面相比，只有第二项的符号不一样。如果我们知道了 $\operatorname{A_1}(\omega_{\frac{n}{2}}^k),\operatorname{A_2}(\omega_{\frac{n}{2}}^k),\omega_n^k$，就可以得到 $\operatorname{A}(\omega_n^k),\operatorname{A}(\omega_n^{k+\frac{n}{2}})$。

考虑递归，每次只需要算出 $k\in[1,\frac{n}{2}]$ 的 $\operatorname{A_1}(\omega_n^k)$ 和 $\operatorname{A_2}(\omega_n^k)$ 以及 $\omega_n^k$，后面的 $k\in(\frac{n}{2},n]$ 的部分就可以直接 $\operatorname{O}(1)$ 得到。

边界条件为 $k=1$，也就变成了单位根，$\operatorname{O}(1)$ 可以计算。换句话说当 $k=1$ 时直接 return，所以最终需要继续计算答案的递归状态是 $n=2$。递归时问题序列每次减半，时间复杂度是 $n\log_2 n$。

快速傅里叶逆变换

子标题是 IFFT 的简写。众所周知我们用 FFT 以 $n\log n$ 的优秀复杂度得到了点值表示法，$\operatorname{\Theta}(n)$ 的更优秀的复杂度算完了乘积多项式的点值表示法，但是转回来还是只能拉格朗日插值 $\operatorname{O}(n^2)$，复杂度级别挂成和暴力 DFT 一样了。

我们不仅仅需要一个快速换过去的方法，还要研究一个快速换回来的方法。为此我们还需要进一步研究 $\omega$ 有关的性质。回到原来的那个图片：

我们现在已经知道了从 $X_0$ 开始逆时针编号，得到的编号为 $k$ 的点对应复数大小是 $\omega_n^k$。在这种方式的定义下，$k$ 一定是大于 $0$ 的。而前面我们已经接受了 $k\geqslant n$ 的情况：相当于转了几圈回来。现在考虑 $k<0$ 的情况。

既然逆时针的编号是 $0,1,2,3,\cdots$，那么编号为负的不就相当于顺时针编号吗：$0,-1,-2,-3,\cdots$。

用另一种方法来看，因为 $\omega_n^k\omega_n^p=\omega_n^{k+p}$，随便代个数：$\omega_n^2\omega_n^{-1}=\omega_n^1$，在图像上面表示，既可以视为是 $X_2$ 顺时针转了一个单位角度（一个单位角度的大小为 $X_1$ 的辐角的大小），也可以视为逆时针转了 $7$ 个单位角度。所以得到一个性质：$\omega_n^{-k}=\omega_n^{n-k}$。

扩展了 $\omega_n^k$ 的定义之后，考虑怎么把一个点值表达式转化为系数表达式。假设我们点值表达式得到的每个点的点值（即 $(x,y)$ 中的 $y$）提出来，组成 $g_i$，系数表达式的每一个系数也提出来，组成 $a_i$。

先上结论，再证明（因为我也不知道怎么推结论只会推式子）

$n\ast a_k=\sum_{p=0}^{n-1} g_p(\omega_n^{-k})^p$

证明：

脑补一下 $g_i$ 都是怎么来的，$g_i$ 是原多项式的点值表示法的点值。点值是 $x=\omega_n^i$ 代入到原多项式得来的。原多项式的系数表示法是 $\sum_{k=0}^{n-1} a_kx^k$。把 $x=\omega_n^i$ 代入：

$g_i=\sum_{k=0}^{n-1} a_k(\omega_n^i)^k$

代入到要证明的式子的等号右边的东西：

$\begin{aligned} \sum_{p=0}^{n-1} g_p(\omega_n^{-k})^p &= \sum_{p=0}^{n-1}\sum_{q=0}^{n-1} a_q(\omega_n^p)^q(\omega_n^{-k})^p \\ &= \sum_{p=0}^{n-1}\sum_{q=0}^{n-1} a_q\omega_n^{pq}\omega_n^{-kp} \\ &= \sum_{p=0}^{n-1}\sum_{q=0}^{n-1} a_q\omega_n^{p(q-k)} \end{aligned}$

分类讨论 $q,k$ 的关系：

1. $q\neq k$ 设 $r=q-k$：

考虑对于一个 $r$，对总和造成的贡献是：

$\begin{aligned} \sum_{p=0}^{n-1} a_{r+k}\omega_n^{pr} &= a_{r+k}\sum_{p=0}^{n-1} \omega_n^p\omega_n^r \\ &= a_{r+k}\omega_n^r \sum_{p=0}^{n-1}\omega_n^p \end{aligned}$

后面的求和就是等比数列求和，已经见很多次了吧，不多赘述了。注意 $\omega_n^n=1$：

$\begin{aligned} \sum_{p=0}^{n-1}\omega_n^p &= \frac{(\omega_n^1)^{n-1+1}-1}{\omega_n^1-1} \\ &= \frac{\omega_n^n-1}{\omega_n^1-1} \\ &= \frac{0}{\omega_n^1-1} \\ &= 0 \end{aligned}$

也就是如果 $q\neq k$，式子对答案没有贡献。

2. $q=k$

$\begin{aligned} \sum_{p=0}^{n-1}\sum_{q=0}^{n-1} a_q\omega_n^{p(q-k)}[q=k] &= \sum_{p=0}^{n-1} a_k\omega_n^0 \\ &= \sum_{p=0}^{n-1} a_k \\ &= na_k \end{aligned}$

因此我们得出，$n\ast a_k=\sum_{p=0}^{n-1} g_p(\omega_n^{-k})^p$。

这个东西的本质就是单位根反演（没错，这玩意居然还能反演……）。

我们再回顾一下怎么用 $a_i$ 求 $g_i$，对比一下怎么用 $g_i$ 求 $a_i$：

$g_i=\sum_{k=0}^{n-1} a_k(\omega_n^i)^k \\ n\ast a_i=\sum_{k=0}^{n-1} g_k(\omega_n^{-i})^k$

不就是原来 FFT 的过程，把代入的数字从 $\omega_n^i$ 换成 $\omega_n^{-i}$，最后结果再除个 $n$ 吗。

因此 IFFT 和 FFT 过程几乎没差多少，一般都写在一个函数里面，开个 $type$ 判断我现在是求 FTT 还是 IFFT。

蝴蝶变换

这是 FFT 的优化。只要你的常数是正常人的常数，没有刻意地卡常，而且写的分治是递归，应该会被模板题卡成 $70\sim 90$，反正过不去。然后你就会发现，题解里面全部不是递归实现的，最优的解法都是循环实现的。

我们考虑分治的时候怎么分治的：我们把下标（从 $0$ 开始标下标）是偶数的提出来放在前面，下标是奇数的提出来放在后面，然后往下递归：

第一行是原数组的下标顺序，最后一行是分治到最终位置后原数组的下标跑到哪里去了。看起来好像顺序挺乱的，但是我们考虑一下第一行和最后一行的二进制：

我们发现每个位置分治后的最终位置为它的二进制倒过来得到的数字对应的位置。

不妨记录好每一个数的最终位置 $rev_i$。$rev_i$ 的求法可以如下考虑：分类讨论 $i$ 的最后一位是不是 $1$：

1. 若最后一位是 $0$，$rev_i$ 相当于 $i$ 除去最后一位剩下的二进制翻过来。
2. 若最后一位是 $1$，$rev_i$ 相当于 $i$ 除去最后一位剩下的二进制翻过来，然后在最高位前面加上 $1$。
3. 递推顺序为从小到大顺序递推。

rev[i]=(rev[i>>1]>>1) | ((i&1)<<(lim-1))
// lim 表示最高位
复制代码

这样子就省去了递归和数组复刻的常数了，时间进一步优化，可以过那个板子题了。

三次变两次优化

本优化默认需要计算的两个多项式系数都是实数。

朴素的 FFT 求两个多项式 $A\ast B$ 的积 $C$（假设题目要求 $C$ 用系数表示法表示），那么就需要先对 $A$ 进行一次 FFT，$B$ 进行一次 FFT，$\operatorname{O}(n)$ 的点值表达式乘法求出 $C$ 后，对 $C$ 进行一次 IFFT。分治的函数（即 FFT 和 IFFT 使用的总次数）一共调用了三次。这玩意怎么优化呢？

回到复数的表示，$z=a+bi$。如果求一下 $z^2$ 会怎样：

$\begin{aligned} (a+bi)^2 &=a^2+(bi)^2+2abi \\ &= a^2-b^2+2abi \end{aligned}$

我们发现后面这个部分有一个 $2ab$，所以如果把 B 的系数当作 A 的系数的虚部，得到的新多项式记作 $A'$，计算一遍 $A'\ast A'$，结果的虚部提出来，也就是 $2ab$，然后除以 $2$，不就得出 $C$ 的系数了吗？

也就是说，把 $B$ 的系数放在 $A$ 的虚部里，得到新的多项式 $A'$，计算 $A'\ast A'$，结果 IFFT 回系数表示法，然后虚部提出来除以 $2$ 就是 $C$ 的系数了。只需要 FFT 一次 $A'$，IFFT 一次 $A'\ast A'$，只需要两次就够了。

常数大概优化为原来的 $\dfrac{2}{3}$。

NTT

与 FFT 的区别

NTT 的复杂度和 FFT 的复杂度都是一样的，都是 $n\log n$，为什么要有这个神奇的 NTT 呢？

（下文加粗部分都是 FFT 的缺点）首先，多项式的问题往往需要对一个数取模（通常是 $998244353$），如果你用的是 FFT，众所周知复数是不能取模的，就会出问题。其次，FFT 用的只能是 double 类型的变量。假设给定的两个多项式，一个系数范围是 $[10^{-6},10^{-5}]$，另一个范围是 $[10^5,10^6]$，中间就跨了 $10^{12}$ 级别的精度，很容易出现精度误差，这种情况基本崩不住。最后，NTT 的运算是基于整数的，FFT 是基于浮点数的，众所周知有四则运算的变量类型里 int 常数最小，因此相比较 NTT，FFT 常数较大。

基于以上三点，人们又想出了一个能取模，没有精度误差，基于整数操作的，系数表示法和点值表示法互转的方法：快速数论变换算法——NTT。

根据名字就可以猜出来，NTT 是基于数论的内容的，别转走，其实没那么多数论内容，而且很多地方和 FFT 是一样的。一般认为 NTT 是 FFT 的可取模升级版。至于效率，小数据下的运行效率和 FFT 相比快了不少，大数据差不太多。

原根

这就是支撑 NTT 的数，它的重要性和 $\omega$ 在 FFT 里面的重要性一样。对于一个给定模数 $p$，如果有一个数 $g\in Z,g^{p-1}\equiv 1\pmod{p}$，其中 $g$ 满足 $\forall 0<i<j<p,g^i\not\equiv g^j\pmod{p}$，则称 $g$ 是 $p$ 的一个原根。

原根的确切的定义如下：

首先定义阶：对于两个数 $a,p$，找到一个最小的整数 $k$ 使得 $a^k\equiv 1\pmod{p}$，则称 $k$ 为 $a$ 在模 $p$ 意义下的阶。用 $\delta_p(a)$（$\delta$ 读作 delta）表示。阶数等于幂的最小循环节（模意义下），不难理解吧。

原根的确切定义是：对于满足 $\delta_p(g)=\varphi(p)$ 的 $g$，我们称 $g$ 为 $p$ 的原根。

只有 $2,4,2p^k,p^k$（其中 $p$ 是奇质数）这些数才有原根。

在 FFT 中，用到了 $\omega_n$ 的以下基本性质：

• $(\omega_n^1) ^k=\omega_n^k$
• $\omega_n^k=\omega_{2n}^{2k}$
• $\omega_n^{k+n}=\omega_n^k$
• $\omega_n^{0\sim(n-1)}$ 互不相同

其他的所有性质都基于这四个性质推出（可以对照前面的 FFT 试着推一下），由第一点、第三点、第四点，综合起来就等价于：$\omega_n^1$ 在模意义下阶恰好为 $n$。我们设题目给定的模数是 $p$，$g$ 是 $p$ 的一个原根。首先，显然 $p-1$ 不一定恰好等于 $n$（$p-1>n$ 看起来好像可以，但实际上无法满足第三点性质，$p-1<n$ 时一定会违背第四点），因此 $g$ 不能直接替代单位根。

至于如何找到阶数为 $n$ 的点，需要用到一个有关阶的性质：若 $p$ 的原根是 $g$，则 $g^k$ 的阶为 $\dfrac{p-1}{\gcd(k,p-1)}$。考虑能否控制 $k$ 的值，使得阶数可以取 $n$：

$\begin{aligned} \frac{p-1}{\gcd(k,p-1)}=n &\Rightarrow \gcd(k,p-1)=\frac{p-1}{n} \\ &\Rightarrow k=\frac{p-1}{n} \end{aligned}$

当 $n \nmid (p-1)$ 时，无解。否则，当 $k=\frac{p-1}{n}$ 时，$g_k$ 的阶恰好为 $n$。既然如此，就说明 $g^k$ 已经满足了 $\omega$ 基本性质的第一、三、四点。还差第二点 $\omega_n^k=\omega_{2n}^{2k}$ 不知道是否满足，考虑证明。

证明：

先考虑 $\omega_n^k=\omega_{2n}^{2k}$ 能否进行一定的转换。我们抛出一个引理：

$\omega_n^k=\omega_{2n}^{2k}$ 是 $(\omega_{2n}^1)^2=\omega_n^1$ 的充要条件。

因为众所周知，$\omega_n^1\cdot \omega_n^1=\omega_n^2$，充分性必要性均显然。

所以证明 $g^k$ 是否满足原性质，即证明其是否满足 $g^k=(g^{\frac{p-1}{2n}})^2$。将 $k=\frac{p-1}{n}$ 代入可得原式成立。证毕。

因此当 $k=\frac{p-1}{n}$ 时，$g^k$ 符合所有的要求，可以用 $g^k$ 代替单位根进行分治。

快速数论变换算法

FFT 之所以可以大大地优化 DFT，是因为傅里叶选的这些数太好了，$\omega$ 有非常多的神奇性质，使得我可以加以利用，进行分治。

但是我们都知道那个圆上面的点，除了四个，别的点的直角坐标都至少有一个是小数，带来了 double 类型尴尬的精度误差。如果我能在整数里面找到一个有相同或类似性质的数，就可以解决这种问题。而刚刚提及的原根就恰好符合这些性质。

在 NTT 下，一般题目都要求我们对一个数取模：$998244353$，而且不会是 $10^7+9,10^9+7$ 之类的这些质数，一定是 $998244353$。事实上，之所以 NTT 的题目模数都是这个数，是因为前面推过 $n\mid (p-1)$ 时，才能用原根构造单位根，而一般分治算法 $n$ 都是 $2$ 的若干次幂，因此我们期望找到一个质数 $p$，使得 $p-1$ 有很多因子 $2$，$998244353=2^{23}\times 7\times 17+1$，很完美地符合条件。再者，这个数的最小原根是 $3$，比较方便。

然后把 FFT 有关单位根的东西全部换成 $g^{\frac{p-1}{n}}$ 即可。注意因为没有了实数和虚数部分，所以没有三变二优化，只有蝴蝶变换结论优化递归常数。

虽然 NTT 的时间复杂度和 FTT 一样，但是没有毒瘤数据的情况下，比 FFT 快了约 $\dfrac{1}{5}$，预处理单位根（减少乘法运算）后又比原来快了约 $\dfrac{1}{5}$。常数已经非常可观了。

upd：封装后实测 NTT 用时为 FFT 的 $70\%$，而且不知道为什么我写的 FFT 空间居然是 NTT 的 $7$ 倍，如图：

MTT

考虑这样一个阴间的毒瘤题。在这里我们需要取模，但是模数并不是很友好的 $998244353$。

这里如果从 NTT 的角度出发，可以选择三个友好的模数，跑九次 NTT，最后的答案用中国剩余定理合并。假设求卷积前值域都小于 $p$，那么卷积后能产生的最大值域在本题就是 $p\times p\times n=10^{26}$，long long 也救不了你，只能手写高精度，不会有人想写吧。

考虑如何用 FFT 的方法去做，常见的方法是拆系数。我们发现值域的范围是 $10^9$，不妨把一个系数 $a$ 拆成 $a=2^{15}p+q$ 的形式（反正 FFT 用的 double，没必要纠结能不能被 $2^{15}$ 整除吧），可以发现这样 $p,q\leqslant 2^{15}$。值域范围缩小为 $p\times n=10^{14}$。可以用 long long 存储，乘法运算时勤快取模防溢出就可以了。

考虑这样子拆掉系数原式会变成什么：

$a_1a_2=(2^{15}p_1+q_1)(2^{15}p_2+q_2)=2^{30}p_1p_2+2^{15}p_1q_2+2^{15}p_2q_1+q_1q_2$

一次 FFT 要拆成四次 FFT，然后没有三到二的 FFT 板子需要跑三次 FFT，拆完之后总的 FFT 和 IFFT 次数就是 $12$ 次，常数关系已经让 $n\log n$ 的 FFT 差点赶上 $n\log^2 n$ 了。

考虑如何优化。现在给定了四个多项式 $A_1,A_2,B_1,B_2$（所有多项式的系数都是实数）。不妨用三变二的思想，开一个多项式 $P$，$P$ 的实部是 $A_1$ 的系数，虚部是 $A_2$ 的系数。同理再开一个多项式 $Q$。

$C_1=P\ast Q=A_1B_1-A_2B_2+A_1B_2i+A_2B_1i$

有四项就很烦，但是考虑到i如果我能算出来一个 $C_2$ 和这里面有几项可以抵消掉，不就省了几次 FFT 吗？考虑再开一个 $R$，实部是 $A_1$，虚部是 $-A_2$。

$C_2=R\ast Q=A_1B_1+A_2B_2+A_1B_2i-A_2B_1i$

$C_1+C_2=2(A_1B_1+A_1B_2i) \\ C_1-C_2=2(A_2B_1i-A_2B_2)$

分别取出这两个结果的实部和虚部就可以得到四个多项式两两的乘积分别是多少。代价是要求出 $P,Q,R$ 三个多项式中的 $P\ast Q,R\ast Q$。八次 FFT（拆完系数的两个多项式乱搞，如果先转可以变成四次）优化成了三次，还行。

把 $A_1,B_1,A_2,B_2$ 当成拆系数后得到的四个多项式，然后利用虚实部转化为 $P,Q,R$ 三个多项式，对这三个求 FFT，得到 $C_1,C_2$，把 $C_1,C_2$ 给 IFFT 回来。然后提取 $a_1b_1,a_1b_2,a_2b_1,a_2b_2$，带到最前面的式子 $2^{30}p_1p_2+2^{15}p_1q_2+2^{15}p_2q_1+q_1q_2$，就解出最后的答案了。

需要进行三次 FFT，两次 IFFT，常数不算那么不可观了。

多项式除法

不是字面意义的多项式除法，字面意义的多项式除法并没有在网络上找到很好的做法。本处的多项式除法参见多项式乘法逆元板子题。

众所周知，除一个数等于乘一个数的逆元。所以多项式除法可以转化为求一个多项式的逆元，然后套多项式乘法板子，FFT,NTT,MTT 复杂度都是 $n\log n$，其中 MTT 的常数是 FFT 的五倍，NTT 的理论常数最小，但是还是具体情况具体分析。

考虑如何求 $A$ 的逆元。

多项式乘法逆元

设 $A$ 在模 $x^n$ 意义下的逆元为 $B_1$，在模 $x^{\left\lceil\frac{n}{2}\right\rceil}$ 意义下的逆元为 $B_2$。根据逆元定义：

$A\ast B_1\equiv 1\pmod{x^n} \\ A\ast B_2\equiv 1\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}$

上面这两个式子就是这样的形式：$ab=xp^2+1,ac=yp+1$，考虑作差：$ab-ac=xp^2-yp=p(xp-y)$，也就是说二者之差是 $p$ 的倍数。即 $a(b-c)\equiv 0\pmod{p}$。

$A\ast(B_1-B_2)\equiv 0\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}} \\ \therefore B_1-B_2\equiv 0\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}$

我们知道，如果 $p$ 可以被 $a$ 整除，那么 $p^2$ 可以被 $a^2$ 整除。换句话说，如果 $p\equiv 0\pmod{a},p^2\equiv 0\pmod{a^2}$。所以：

$(B_1-B_2)^2\equiv 0\pmod{x^n} \\ B_1^2+B_2^2-2B_1B_2\equiv 0\pmod{x^n} \\ A\ast(B_1^2+B_2^2-2B_1B_2)\equiv 0\pmod{x^n}$

因为：$A\ast B_1\equiv 1\pmod{x^n}$，所以乘进去：

$B_1-2B_2+AB_2^2\equiv 0\pmod{x^n} \\ B_1\equiv 2B_2-AB_2^2\equiv 0\pmod{x^n}$

也就是说，对于 $A\ast B_1\equiv 1\pmod{x^n}$ 的求 $B_1$，可以通过求出使 $A\ast B_2\equiv 1\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}$ 的 $B_2$，然后加上多项式乘法来求。

从上往下递归的常数传参要传多项式，非常恶心，不如直接从下往上递推，从 $x^1$ 推到第一个 $x^{2^k}$，满足 $2^k\geqslant n$ 的 $k$ 就行了。

支撑主体的是 NTT 求逆，还用到了倍增。整体的复杂度是 $n\log n$。

代码

多项式 FFT，NTT，MTT，求逆代码：

namespace DXS{ // 封装
	inline int mi(int x,int p){ // 快速幂
		int res=1;
		while(p){
			if(p&1) res*=x,res%=mod;
			x*=x,x%=mod,p>>=1;
		}
		return res;
	}
	struct com{ // 手写复数结构体，支持 +-* 三种运算
		double a,b; // a+bi
		com(double x=0,double y=0){a=x,b=y;}
		com operator + (const com b) const{
			com a=*this,res;
			res.a=a.a+b.a,res.b=a.b+b.b;
			return res;
		}
		com operator - (const com b) const{
			com a=*this,res;
			res.a=a.a-b.a,res.b=a.b-b.b;
			return res;
		}
		com operator * (const com b) const{
			com a=*this,res;
			res.a=a.a*b.a-a.b*b.b,
			res.b=a.a*b.b+a.b*b.a;
			return res;
		}
	};

	int rev[max_n*2]; // 蝴蝶变换
	const int mod=998244353,G=3,IG=332748118; // NTT 模数，原根和原根的逆元
	const double Pi=acos(-1.0); // 圆周率，计算单位根用

	inline void NTT(int *g,int n,bool op){ // 传入数组 g，op=1 表示转点值，0 表示转系数
		for(register int i=0;i<n;++i) if(i<rev[i])
			swap(g[i],g[rev[i]]);
		for(register int mid=1;mid<n;mid<<=1LL){
			int omega=mi(op?G:IG,(mod-1)/(mid+mid));
			for(register int i=0;i<n;i+=mid+mid){
				int w=1;
				for(register int j=0;j<mid;++j,w=(w*omega)%mod){
					int x=g[i+j],y=w*g[i+j+mid]%mod;
					g[i+j]=(x+y)%mod,
					g[i+j+mid]=(x-y+mod)%mod;
				}
			}
		}
		if(!op){
			int iv=mi(n,mod-2);
			for(register int i=0;i<n;++i)
				g[i]=g[i]*iv%mod;
		}
	}
	inline void FFT(com *g,int n,bool op){ // 传入数组 g，op=1 表示转点值，0 表示转系数
		for(register int i=0;i<n;++i) if(i<rev[i])
			swap(g[i],g[rev[i]]);
		for(register int mid=1;mid<n;mid<<=1LL){
			com omega(cos(Pi/mid),sin(Pi/mid)*(op?1:-1));
			for(register int i=0;i<n;i+=mid+mid){
				com w(1,0);
				for(register int j=0;j<mid;++j,w=w*omega){
					com x=g[i+j],y=w*g[i+j+mid];
					g[i+j]=x+y,
					g[i+j+mid]=x-y;
				}
			}
		}
	}
	inline void mulNTT(int *f,int *g,int n,int m){ // 传入系数表示数组 f,g，计算 f=f*g
		int len=1;
		while(len<n+m) len<<=1LL;
		for(register int i=0;i<len;++i)
			rev[i]=(rev[i>>1]>>1)|((i&1)?len>>1:0);
		NTT(f,len,1),NTT(g,len,1);
		for(register int i=0;i<len;++i)
			f[i]*=g[i],f[i]%=mod;
		NTT(f,len,0);
	}
	inline void mulFFT(int *f,int *g,int n,int m){ // 传入系数表示数组 f,g，计算 f=f*g
		static com F[max_n],G[max_n];
		for(register int i=0;i<n;++i)
			F[i].a=f[i],F[i].b=0;
		for(register int i=0;i<m;++i)
			G[i].a=g[i],G[i].b=0;
		int len=1;
		while(len<n+m) len<<=1LL;
		for(register int i=0;i<len;++i)
			rev[i]=(rev[i>>1]>>1)|((i&1)?len>>1:0);
		FFT(F,len,1),FFT(G,len,1);
		for(register int i=0;i<len;++i)
			F[i]=F[i]*G[i];
		FFT(F,len,0);
		for(register int i=0;i<len;++i)
			f[i]=(int)(F[i].a/len+0.5);
	}
	const int K=32768,S=K*K;
	inline void mulMTT(int *f,int *g,int n,int m,int mod){ // 传入系数表示数组 f,g，计算 f=f*g
		static com P[max_n],Q[max_n],R[max_n];
		memset(P,0,sizeof(P)),
		memset(Q,0,sizeof(Q)),
		memset(R,0,sizeof(R));
		for(register int i=0;i<n;++i){
			P[i].a=f[i]/K,
			P[i].b=f[i]%K;
			R[i].a=f[i]/K,
			R[i].b=-f[i]%K;
		}
		for(register int i=0;i<m;++i){
			Q[i].a=g[i]/K,
			Q[i].b=g[i]%K;
		}
		int len=1;
		while(len<n+m) len<<=1LL;
		for(register int i=0;i<len;++i)
			rev[i]=(rev[i>>1]>>1)|((i&1)?len>>1:0);

		FFT(P,len,1),FFT(Q,len,1),FFT(R,len,1);
		for(register int i=0;i<len;++i){
			Q[i].a/=len,
			Q[i].b/=len;
		}
		for(register int i=0;i<len;++i){
			P[i]=P[i]*Q[i],
			R[i]=R[i]*Q[i];
		}
		FFT(P,len,0),FFT(R,len,0);
		for(register int i=0,a1b1=0,a1b2=0,a2b1=0,a2b2=0;i<len;++i){
			a1b1=(int)floor((P[i].a+R[i].a)/2+0.5)%mod,
			a1b2=(int)floor((P[i].b+R[i].b)/2+0.5)%mod,
			a2b1=((int)floor(P[i].b+0.5)-a1b2)%mod,
			a2b2=((int)floor(R[i].a+0.5)-a1b1)%mod;
			f[i]=(a1b1*S%mod+(a1b2+a2b1)%mod*K%mod+a2b2)%mod;
			while(f[i]<0) f[i]+=mod;
			f[i]%=mod;
		}
	}
	inline void invNTT(int *f,int *g,int n){ // 用 NTT 计算逆元
		if(n==1){
			g[0]=mi(f[0],mod-2);
			return;
		}
		invNTT(f,g,(n+1)>>1);
		int len=1;
		while(len<n+n) len<<=1LL;
		for(register int i=0;i<len;++i)
			rev[i]=(rev[i>>1]>>1)|((i&1)?len>>1:0);
		static int h[max_n];memset(h,0,sizeof(h));
		for(register int i=0;i<n;++i)
			h[i]=f[i];
		NTT(h,len,1),NTT(g,len,1);
		for(register int i=0;i<len;++i)
			g[i]*=(2-h[i]*g[i]%mod+mod)%mod,
			g[i]%=mod;
		NTT(g,len,0);
		for(register int i=n;i<len;++i) g[i]=0;
	}
	inline void invMTT(int *f,int *g,int n,int mod){ // 用 MTT 计算逆元
		if(n==1){
			g[0]=mi(f[0],mod-2);
			return;
		}
		invMTT(f,g,(n+1)>>1,mod);
		int len=1;
		while(len<n+n) len<<=1LL;
		static int h[max_n*2];memset(h,0,sizeof(h));
		for(register int i=0;i<len;++i)
			h[i]=f[i];
		mulMTT(h,g,n,n,mod),
		mulMTT(h,g,n,n,mod);
		for(register int i=0;i<n;++i)
			g[i]=(g[i]*2-h[i]+mod)%mod;
		for(register int i=n;i<len;++i) g[i]=0;
	}
}
复制代码