四大变换学习笔记（FFT，NTT，FMT，FWT）
FFT
介绍 FFT 之前，我们先介绍一下多项式。
多项式
多项式的概念
多项式是一个形如 \(A(x)=a_0+a_1x+a_2x^2+a_3x^3+\cdots+a_nx^n\) 的式子，我们称此时的 \(A(x)\) 为一个 \(n\) 次多项式，用更数学的语言来讲就是：
\[A(x)=\sum_{i=0}^na_ix^i \]可以发现 \(A(x)\) 的本质是一个函数，所以我们也称 \(A(x)\) 为一个 \(n\) 次函数，即 \(n\) 次多项式函数，我们发现当 \(n\) 取 \(0\) 时该函数是一个常量，而当 \(n\) 取 \(1\) 或 \(2\) 时则是我们熟悉的一次函数和二次函数。
其中我们约定 \(A(x)_i\) 表示多项式函数 \(A(x)\) 的第 \(i\) 项，即上面的 \(a_i\)，并且约定如果 \(i>n\)，则 \(n\) 次多项式的第 \(i\) 项为 \(0\)。
多项式的四则运算
多项式的加法与减法非常简单，两个多项式的和差就是相对应的项相减的结果，即若 \(C(x)=A(x)\pm B(x)\)，则：
\[C(x)_i=A(x)_i\pm B(x)_i \]来看多项式的乘法，其实与我们初中对多项式的乘法处理相似，回忆一下：
\[\begin{aligned} &(a_0+a_1x+a_2x^2)(b_0+b_1x)\\ =&a_0b_0+a_0b_1x+a_1b_0x+a_1b_1x^2+a_2b_0x^2+a_2b_1x^3\\ =&a_0b_0+(a_0b_1+a_1b_0)x+(a_1b_1+a_2b_0)x^2+a_2b_1x^3 \end{aligned} \]总结一下就是若 \(C(x)=A(x)\times B(x)\)，则：
\[C(x)_i=\sum_{j=0}^iA(x)_jB(x)_{i-j} \]我们也将多项式的乘法称为多项式卷积。
多项式的除法涉及到多项式的乘法逆元，这里暂且跳过。
FFT
接下来讲 FFT，我们发现对于多项式加减法来说，时间复杂度都是 \(\mathcal O\left(n\right)\) 的，是一个比较优的复杂度，相对而言，多项式卷积的时间复杂度就偏劣，是 \(\mathcal O\left(n^2\right)\) 的，令人无法接受，考虑有没有更优的办法，其实是有的，这里引入一个点值表示法的概念。
初中时有用两点确定一条直线以及三点确定一条抛物线的情况，即用两点确定一个一次函数，三点确定一个二次函数，相对的我们想到可不可以用 \(n+1\) 个点确定一个 \(n\) 次函数呢？
这实际上是可以的，因为一个 \(n\) 次函数本质上就是由 \(n+1\) 个系数组成，所以相当于列了 \(n+1\) 个线性方程，所以用高斯消元是可以解的。
所以我们也可以用 \(n+1\) 个点来表示一个 \(n\) 次多项式，我们称这种用点来表示多项式的办法叫点值表示法，相对的我们将之前用 \(n+1\) 个系数表示的方法叫做系数表示法。
我们发现如果我们想要知道两个多项式的乘积，我们其实可以找出若干点，用两个多项式分别算出这些点对应的点值，然后将这两个点值乘起来，我们就得到了乘积多项式的点值表示，我们再把它转回系数表示就好了。
这样的话我们求乘积的时间复杂度就变成 \(\mathcal O\left(n\right)\) 的了，但是我们并没有什么好的办法将多项式从系数表示转为点值表示，暴力求解时间复杂度还是 \(\mathcal O\left(n^2\right)\) 的，得不偿失。
所以我们需要这种算法—— FFT。
FFT 能够在 \(\mathcal O\left(n\log n\right)\) 时间复杂度内将一个系数表示转为点值表示，是多项式转化的利器。
我们发现如果我们随意的选择我们求值的点，肯定是不行的，要考虑取一些点使它们满足某种优美的性质。
首先，我们约定我们选取的点值为 \(x_0\sim x_n\)，令 \(\operatorname{FFT}(A)_i=A(x_i)=\sum_{j=0}^na_jx_i^j\)。
我们可以随意将 \(A(x)\) 的系数奇偶分开，即找到任意一个 \(m\ge\dfrac n2\)，然后令：
\[\begin{aligned} &A_0(x)=\sum_{i=0}^mA(x)_{2i}x\\ &A_1(x)=\sum_{i=0}^mA(x)_{2i+1}x\\ \end{aligned} \]那么前面的 \(\operatorname{FFT}(A)_i\) 就可以表示为：
\[\operatorname{FFT}(A)_i=A(x_i)=A_0(x_i^2)+A_1(x_i^2)x_i \]发现每次如果我们将一个多项式分成两半，多项式次数规模就减少一半，如果我们把待求点值数也减半，因为我们发现上式的合并对于每个点值是 \(\mathcal O\left(1\right)\)，也就是说总时间复杂度是 \(\mathcal O\left(\text{len}\right)\) 的，那么我们的总时间复杂度就可以控制在 \(\mathcal O\left(n\log n\right)\)。
我们发现点值从 \(A\) 到 \(A_0\) 和 \(A_1\) 时从 \(x\) 变成了 \(x^2\)，也就是说如果我们找到一类数，它们的平方的数量只有原数量的一半，我们就成功了，第一个想法是相反数，但是取完一次就不行了，看来实数域没有这种东西，我们可以考虑从复数域找。
然后有一个东西惊人的与我们的目的吻合——单位根，我们可以先来认识一下这位救星。
单位根 \(\large\omega\)
再讲单位根之前得先讲一讲复数。
我们发现在实数域里面对负数开方是无解的，这究其原因都是因为 \(\sqrt{-1}\) 是未被定义的，因为任何一个满足 \(a>0\) 的 \(\sqrt{-a}\) 都可以写成 \(\sqrt a\cdot\sqrt{-1}\)，所以我们可以强制规定 \(i=\sqrt{-1}\)，我们一般称这个 \(i\) 为虚数单位，然后我们就可以将任何一个实数开方了。
我们可以将任意一个复数写成形如 \(a+bi\) 的形式，其中 \(a\) 叫实部，\(b\) 叫虚部，且因为这里的 \(a\)，\(b\) 都是实数，所以实数的运算法则依然成立，所以复数的加减乘可以优美地定义如下：
\[(a+bi)\pm(c+di)=(a\pm c)+(b\pm d)i\\ (a+bi)(c+di)=ac+adi+bci+bdi^2=(ac-bd)+(ad+bc)i \]我们之前发现一个事实——相反数的平方相等，其中的一个特例更为优美，就是 \(1\) 和 \(-1\)，它们的平方都是 \(1\)，仅仅是减少了一半的数，我们可以看成是在解 \(x^2=1\) 这样一个方程，它的解是 \(x_1=1,x_2=-1\)，那么拓展延伸就可以得到在我们尝试将奇偶分离第二轮时我们就是在求一个方程 \(x^4=1\)，然后发现在实数域的解还是只有 \(1\) 和 \(-1\)，但是我们可以换一个角度思考，我们等于说是求 \((x^2)^2=1\)，那么就是 \(x^2=\pm1\)，也就是说我们还有两个根被舍了就是 \(x^2=-1\) 的两个根，就是虚数单位 \(i=\sqrt{-1}\)，即 \(x_1=i,x_2=-i\)，然后其实我们可以在虚数域上对这个数再开根，可以无限次数开下去。
然后发现这些数一定存在一个 \(n\) 满足它的 \(n\) 次方为 \(1\)，所以我们称满足 \(\omega_n^n=1\) 的所有 \(\omega_n\) 为 \(n\) 次单位根。
我们发现找高次单位根是一件很麻烦的事情，考虑有没有更简单的办法，就要引入复数的三角表示。
数形结合一直是数学研究的利器，当我们研究实数时，我们利用数轴这一工具，并且以从原点出发的有向线段的长度表示一个数，相对的，我们发现其实每个虚数都可以表示成两个实数的组合，所以我们可以考虑将虚数表示在平面直角坐标系上，用 \(x\) 轴表示实部， \(y\) 轴表示虚部，然后用向量来表示一个复数，即用 \(\overrightarrow a=(a,b)\) 来表示复数 \(a+bi\)，这样的话，称该向量与 \(x\) 轴的夹角为辐角 \(\theta\)，可以得到 \(a+bi=|\overrightarrow a|(\cos\theta+\sin\theta i)\)。
这种表示法的其中一个优点是可以直接将复数加减转化为向量加减，并且乘法也非常简单。
令 \(a_1+b_1i\)，\(a_2+b_2i\) 对应的向量分别为 \(\overrightarrow a_1\)，\(\overrightarrow a_2\)，辐角分别为 \(\theta_1\)，\(\theta_2\)，\((a_1+b_1i)(a_2+b_2i)\) 的乘积表示成：
\[\begin{aligned} &|\overrightarrow a_1||\overrightarrow a_2|(\cos\theta_1+\sin\theta_1i)(\cos\theta_2+\sin\theta_2i)\\ =&|\overrightarrow a_1||\overrightarrow a_2|((\cos\theta_1\cos\theta_2-\sin\theta_1\sin\theta_2)+(\cos\theta_1\sin\theta_2+\sin\theta_1\cos\theta_2)i)\\ =&|\overrightarrow a_1||\overrightarrow a_2|(\cos(\theta_1+\theta_2)+\sin(\theta_1+\theta_2)i) \end{aligned} \]讲人话就是模长相乘，辐角相加。
我们发现这种方式可以很方便的求一个复数的高次幂，将对应向量的模长简写为 \(r\) ，则一个复数 \((r(\cos\theta+\sin\theta i))^k\) 为：
\[(r(\cos\theta+\sin\theta i))^k=r^k(\cos(k\theta)+\sin(k\theta)i) \]这样我们就可以非常方便地求 \(\omega_n\)，即 \(\omega_n=(\cos\dfrac{2k\pi}n+\sin\dfrac{2k\pi}ni),k\in\mathbb N\cap[0,n)\)。
这里有很多性质需要记一下：
- \(\omega_n^2=\omega_\frac n2\)，而且因为 \(n\) 次单位根有 \(n\) 个，所以肯定是每两个 \(\omega_n^2\) 对应一个 \(\omega_\frac n2\)。
- \(\forall n>1,\ \omega_n^\frac n2=-1\)。
- 所有的 \(\omega_n\) 沿单位圆均匀分布。
- \(\omega_n^k\) 和 \(\omega_n^{-k}\) 的虚部互为相反数
- \(\omega_n^k\) 的值存在长度为 \(n\) 的循环节，即 \(\omega_n^k=\omega_n^{k\bmod n}\)。
DFT/IDFT
因为我们每次都是将范围减半，所以我们可以将 \(n\) 扩大到最近的 \(2^m-1\)，然后我们令 \(n=2^m-1\)，这样我们可以假设它是一个 \(n-1\) 次多项式。
我们可以令 \(x_i=\omega_n^i\) ，再带回前面的 FFT 的式子中得到：
\[\begin{aligned} \operatorname{FFT}(A)_i&=A(\omega_n^i)\\ &=A_0(\omega_n^{2i})+A_1(\omega_n^{2i})\omega_n^i\\ &=A_0(\omega_\frac n2^i)+A_1(\omega_\frac n2^i)\omega_n^i \end{aligned} \]然后发现存在 \(i<\dfrac n2\) 和 \(i\ge\dfrac n2\) 两种情况，我们钦定 \(i<\dfrac n2\)，令后一种情况为 \(i+\dfrac n2\)，我们把后一个情况拿出来考虑：
\[\begin{aligned} \operatorname{FFT}(A)_{\frac n2+i}&=A_0(\omega_\frac n2^i)+A_1(\omega_\frac n2^i)\omega_n^{\frac n2+i}\\ &=A_0(\omega_\frac n2^i)-A_1(\omega_\frac n2^i)\omega_n^i \end{aligned} \]这样就可以每次用一个点值和奇偶两边的计算结果计算出两个点值，至此我们就成功的完成了被称为 DFT 的过程，就是从系数表示到点值表示的过程。
那么我们就已经离胜利只剩一步了，就是把点值表示转回系数表示。
这里想要一步步退出结果还是很烦的，所以假设我们已经通过某种东方的神秘力量得到了答案：
\[A(x)_i=\dfrac1n\sum_{j=0}^{n-1}\operatorname{FFT}(A)_j\omega_n^{-ij} \]我们可以试着去证明一下这个东西：
\[\begin{aligned} A(x)_i&=\dfrac1n\sum_{j=0}^{n-1}\operatorname{FFT}(A)_j\omega_n^{-ij}\\ &=\dfrac1n\sum_{j=0}^{n-1}\sum_{k=0}^{n-1}A(x)_k\omega_n^{jk}\omega_n^{-ij}\\ &=\dfrac1n\sum_{k=0}^{n-1}A(x)_k\sum_{j=0}^{n-1}\omega_n^{j(k-i)} \end{aligned} \]我们可以显然地发现，上式中当 \(k=i\) 时满足 \(\sum_{j=0}^{n-1}\omega_n^{j(k-i)}=\sum{j=0}^{n-1}\omega_n^0=n\)，所以就只要证明当 \(k\neq i\) 时，\(\sum_{j=0}^{n-1}\omega_n^{j(k-i)}=0\) 就好了。
这其实也很简单，因为我们利用上面说道单位根在单位圆上均匀分布的性质，可以发现 \(k-i\) 取便所有 \(\omega_n^1\) 到 \(\omega_n^{n-1}\)，正好两两抵消。
证明完成，我们观察我们 DFT 和 IDFT 的两个式子：
\[\begin{aligned} \operatorname{FFT}(A)_i&=\sum_{j=0}^{n-1}A(x)_j\omega_n^{ij}\\ A(x)_i&=\dfrac1n\sum_{j=0}^{n-1}\operatorname{FFT}(A)_j\omega_n^{-ij} \end{aligned} \]我们发现 IDFT 的过程也可以像 DFT 一样用 FFT 求解，实际上在代码实现上更为简单，因为 \(\omega_n^k\) 和 \(\omega_n^{-k}\) 的差别只在于虚部的正负，只要传一个标记进去，最后再乘上 \(\dfrac1n\) 就好了。
这样就完成了整个 FFT 求解多项式卷积的过程。
例题是P3803 【模板】多项式乘法（FFT），代码实现：
FFT 多项式卷积
#include
using namespace std;
#define Reimu inline void // 灵梦赛高
#define Marisa inline int // 魔理沙赛高
#define Sanae inline bool // 早苗赛高
typedef long long LL;
typedef unsigned long long ULL;
typedef pair Pii;
typedef tuple Tiii;
#define fi first
#define se second
struct Complex {
double r, i;
inline Complex operator+(const Complex &o) const { return {r + o.r, i + o.i}; }
inline Complex operator-(const Complex &o) const { return {r - o.r, i - o.i}; }
inline Complex operator*(const Complex &o) const { return {r * o.r - i * o.i, r * o.i + i * o.r}; }
friend istream &operator>>(istream &cin, Complex &x) { return cin >> x.r; }
friend ostream &operator<<(ostream &cout, const Complex &x) { return cout << '(' << x.r << ',' << x.i << ')'; }
};
const int N = 3000010;
const double Pi = acos(-1);
int n, m, lim;
Complex A[N], B[N], t[N];
Reimu FFT(Complex *F, int len, int ty = 1) {
if (len == 1) return;
int half = len >> 1;
for (int i = 0; i < half; ++i) t[i] = F[i << 1], t[half + i] = F[i << 1 | 1];
memcpy(F, t, sizeof(Complex) * len);
FFT(F, half, ty); FFT(F + half, half, ty);
Complex w = {cos(Pi / half), sin(Pi / half) * ty}, wi = {1};
for (int i = 0; i < half; ++i, wi = wi * w) {
Complex F0 = F[i], F1 = F[half + i] * wi;
F[i] = F0 + F1; F[half + i] = F0 - F1;
}
}
int main() {
ios::sync_with_stdio(false); cin.tie(nullptr); cout.tie(nullptr);
cin >> n >> m;
for (int i = 0; i <= n; ++i) cin >> A[i];
for (int i = 0; i <= m; ++i) cin >> B[i];
lim = __lg(n + m) + 1;
FFT(A, 1 << lim); FFT(B, 1 << lim);
for (int i = 0; i < 1 << lim; ++i) A[i] = A[i] * B[i];
FFT(A, 1 << lim, -1);
for (int i = 0; i <= n + m; ++i) cout << int(A[i].r / (1 << lim) + .5) << ' ';
return 0;
}
然后 FFT 有一个优化，就是利用复数乘法带来的复数平方的性质：
\[(a+bi)^2=(a^2-b^2,2ab) \]我们发现如果我们将需要相乘的两个实数多项式放在同一个复数多项式的实部和虚部，平方完后将虚部就是我们需要的结果，这样我们就可以只做两遍 FFT。
还有一个节省代码的优化，就是关于 IDFT 中的式子，因为 \(\omega_n^{-i}=\omega_n^{n-i}\)，所以实际上我们只要把点值表示再扔进 DFT 中出来的结果除了第 \(0\) 项前后翻转就好了，就得到一个更简洁更快的代码：
FFT plus
#include
using namespace std;
#define Reimu inline void // 灵梦赛高
#define Marisa inline int // 魔理沙赛高
#define Sanae inline bool // 早苗赛高
typedef long long LL;
typedef unsigned long long ULL;
typedef pair Pii;
typedef tuple Tiii;
#define fi first
#define se second
struct Complex {
double r, i;
inline Complex operator+(const Complex &o) const { return {r + o.r, i + o.i}; }
inline Complex operator-(const Complex &o) const { return {r - o.r, i - o.i}; }
inline Complex operator*(const Complex &o) const { return {r * o.r - i * o.i, r * o.i + i * o.r}; }
friend istream &operator>>(istream &cin, Complex &x) { return cin >> x.r; }
friend ostream &operator<<(ostream &cout, const Complex &x) { return cout << '(' << x.r << ',' << x.i << ')'; }
};
const int N = 3000010;
const double Pi = acos(-1);
int n, m, lim;
Complex A[N], t[N];
Reimu FFT(Complex *F, int len) {
if (len == 1) return;
int half = len >> 1;
for (int i = 0; i < half; ++i) t[i] = F[i << 1], t[half + i] = F[i << 1 | 1];
memcpy(F, t, sizeof(Complex) * len);
FFT(F, half); FFT(F + half, half);
Complex w = {cos(Pi / half), sin(Pi / half)}, wi = {1};
for (int i = 0; i < half; ++i, wi = wi * w) {
Complex F0 = F[i], F1 = F[half + i] * wi;
F[i] = F0 + F1; F[half + i] = F0 - F1;
}
}
int main() {
ios::sync_with_stdio(false); cin.tie(nullptr); cout.tie(nullptr);
cin >> n >> m;
for (int i = 0; i <= n; ++i) cin >> A[i].r;
for (int i = 0; i <= m; ++i) cin >> A[i].i;
lim = __lg(n + m) + 1;
FFT(A, 1 << lim);
for (int i = 0; i < 1 << lim; ++i) A[i] = A[i] * A[i];
FFT(A, 1 << lim); reverse(A + 1, A + (1 << lim));
for (int i = 0; i <= n + m; ++i) cout << int(A[i].i / (1 << lim + 1) + .5) << ' ';
return 0;
}
但这个速度其实还是很慢，其重点在于递归的巨大额外开销以及每次递归都要进行内存拷贝，时间上有许多无用开销，可以考虑进行优化，所以就请到我们的重量级优化——迭代 FFT。
我们考虑类似手模递归，FFT 的本质是递归向下走时将数列奇偶分类，向上回溯时每次将两个区间合并，我们考虑模拟这个过程。
合并难度不大，阻碍我们的就是奇偶分类部分，每次奇偶分类都会把数列“打乱”，这给我们造成了一定的麻烦，我们考虑把这个过程写出来：
\[\begin{aligned} \text{len}=8&:0\ 1\ 2\ 3\ 4\ 5\ 6\ 7\\ \text{len}=4&:0\ 2\ 4\ 6\ 1\ 3\ 5\ 7\\ \text{len}=2&:0\ 4\ 2\ 6\ 1\ 5\ 3\ 7\\ \text{len}=1&:0\ 4\ 2\ 6\ 1\ 5\ 3\ 7\\ \end{aligned} \]好像注意到什么，我们可以把二进制表示写出来：
\[\begin{aligned} \text{begin}&:0\ 1\ 2\ 3\ 4\ 5\ 6\ 7&000\ 001\ 010\ 011\ 100\ 101\ 110\ 111\\ \text{end}&:0\ 4\ 2\ 6\ 1\ 5\ 3\ 7&000\ 100\ 010\ 110\ 001\ 101\ 011\ 111 \end{aligned} \]这样看应该很明显了，就是它的位置和它的原位置的二进制表示下正好是翻转，其实也很好理解，你第一轮时你按照最低位分类，将 \(0\) 放前面，\(1\) 放后面，之后按照从低到高为一点点分类，等于是按照二进制翻转后排序，所以我们可以预处理翻转后的位置，当要做 FFT 时先把式子转成递归至底层时的样子，在自下而上循环求解，答案就出来了，而且跑得飞快。
还有一个很小的优化，既然已经预处理了，就把单位根一块与处理出来，因为理论上单位根每次求三角函数也非常慢。
迭代 FFT
#include
using namespace std;
#define Reimu inline void // 灵梦赛高
#define Marisa inline int // 魔理沙赛高
#define Sanae inline bool // 早苗赛高
typedef long long LL;
typedef unsigned long long ULL;
typedef pair Pii;
typedef tuple Tiii;
#define fi first
#define se second
struct Complex {
double r, i;
inline Complex operator+(const Complex &o) const { return {r + o.r, i + o.i}; }
inline Complex operator-(const Complex &o) const { return {r - o.r, i - o.i}; }
inline Complex operator*(const Complex &o) const { return {r * o.r - i * o.i, r * o.i + i * o.r}; }
friend istream &operator>>(istream &cin, Complex &x) { return cin >> x.r; }
friend ostream &operator<<(ostream &cout, const Complex &x) { return cout << '(' << x.r << ',' << x.i << ')'; }
};
const int N = 3000010;
const double Pi = acos(-1);
int n, m, lim;
int rev[N];
Complex A[N], w[N];
Reimu init(int len) {
lim = __lg(len) + 1;
for (int i = 1; i < 1 << lim; ++i) rev[i] = rev[i >> 1] >> 1 | (i & 1) << lim - 1;
w[0] = {1}; w[1] = {cos(Pi / (1 << lim - 1)), sin(Pi / (1 << lim - 1))}; for (int i = 2; i < 1 << lim; ++i) w[i] = w[i - 1] * w[1];
}
Reimu FFT(Complex *F) {
for (int i = 1; i < 1 << lim; ++i) if (i < rev[i]) swap(F[i], F[rev[i]]);
for (int len = 1; len <= lim; ++len) {
int half = 1 << len - 1;
for (int i = 0; i < 1 << lim; i += 1 << len) {
for (int j = 0; j < half; ++j) {
Complex F0 = F[i + j], F1 = F[i + half + j] * w[j << lim - len];
F[i + j] = F0 + F1; F[i + half + j] = F0 - F1;
}
}
}
}
int main() {
ios::sync_with_stdio(false); cin.tie(nullptr); cout.tie(nullptr);
cin >> n >> m;
for (int i = 0; i <= n; ++i) cin >> A[i].r;
for (int i = 0; i <= m; ++i) cin >> A[i].i;
init(n + m);
FFT(A);
for (int i = 0; i < 1 << lim; ++i) A[i] = A[i] * A[i];
FFT(A); reverse(A + 1, A + (1 << lim));
for (int i = 0; i <= n + m; ++i) cout << int(A[i].i / (1 << lim + 1) + .5) << ' ';
return 0;
}
这样的话 FFT 最最基本的应用就结束了。
NTT
NTT 是基于 FFT 基础上的一种多项式卷积算法，其实本质上和 FFT 是一样的，只是因为复数有许多不优美的性质，比如说运算速度慢，计算有较大的精度误差，或者说实数域算法都有这个缺点，所以我们考虑在更优美的一种数域里进行多项式卷积——取模域下。
取模的用途在于将整个有理数域以及部分无理数域和虚数域有重映射到一个有限自然数域中，并且存在许多优美的性质，我们可以考虑在这个数域里去做多项式卷积。
模意义下多项式乘法
延续 FFT 的思路，我们是要找到一种数至少有 \(n\) 个，每两个数的平方相等，平方结果再平方仍满足前面性质，知道只剩一个数。
正好在取模域下就有这种优秀的数——原根。
先讲讲原根的定义：
模 \(p\) 意义下的原根 \(g\) 满足 \(g^{\varphi(p)}=1\) 且 \(\forall x<\varphi(p),\ g^x\neq1\)。
具体怎么找原根这里不多做赘述，在这里我们只需要运用原根这个优美的性质就好了，用类似上面 FFT 中的语言来说就是模 \(p\) 意义下存在 \(\varphi(p)\) 次单位根，当然任何 \(\varphi(p)\) 的因数次的单位根也存在，但是如果 \(\varphi(p)\) 的质因数中 \(2\) 的个数太少也不行，因为这样的话我们就只能奇偶分类没几次就停下了。
其中又有一个得天独厚的模数 \(998244353\)，这个模数有多优美呢？\(998244353=2^{23}\times7\times17+1\)！
也就是说我们可以奇偶分类 \(23\) 次，这是完全够的，\(2^23>10^6\)，而 \(10^6\) 也基本是 \(\mathcal O\left(n\log n\right)\) 能处理的极限了。
作为优秀的 NTT 模数，我们应该记住它的原根是 \(3\) 和 \(114514\)（有点臭）。
其实接下来就很简单了，把 FFT 模板改一下，把单位根改成原根，然后注意原根没有上面单位根的对称性，也不存在实部虚部，所以不能后面的优化都不能加，唯一的优化是将递归改成迭代，因为这个优化本身和数的性质没关系，上面有一个 \(\omega^-k\)，在 NTT 中就直接改成原根的逆元的 \(k\) 次方就好了。
例题依然是P3803 【模板】多项式乘法（FFT），这里就只放迭代 NTT 的代码了。
迭代 NTT
#include
using namespace std;
#define Reimu inline void // 灵梦赛高
#define Marisa inline int // 魔理沙赛高
#define Sanae inline bool // 早苗赛高
typedef long long LL;
typedef unsigned long long ULL;
typedef pair Pii;
typedef tuple Tiii;
#define fi first
#define se second
const int P = 998244353, G = 3, Gi = 332748118;
Marisa add(int x, int y) { return (x += y) < P ? x : x - P; }
Marisa&Add(int&x, int y) { return (x += y) < P ? x : x -= P; }
Marisa sub(int x, int y) { return (x -= y) < 0 ? x + P : x; }
Marisa&Sub(int&x, int y) { return (x -= y) < 0 ? x += P : x; }
Marisa mul(int x, int y) { return 1LL * x * y % P; }
Marisa&Mul(int&x, int y) { return x = 1LL * x * y % P; }
Marisa qpow(int x, int y) { int res = 1; for (; y; y >>= 1, Mul(x, x)) if (y & 1) Mul(res, x); return res; }
const int N = 3000010;
int n, m, lim;
int A[N], B[N], rev[N], g[N], gi[N];
Reimu init(int len) {
lim = __lg(len) + 1;
for (int i = 1; i < 1 << lim; ++i) rev[i] = rev[i >> 1] >> 1 | (i & 1) << lim - 1;
g[0] = gi[0] = 1; g[1] = qpow(G, P - 1 >> lim); gi[1] = qpow(Gi, P - 1 >> lim);
for (int i = 2; i < 1 << lim; ++i) g[i] = mul(g[i - 1], g[1]), gi[i] = mul(gi[i - 1], gi[1]);
}
Reimu NTT(int *F, int *G = g) {
for (int i = 1; i < 1 << lim; ++i) if (i < rev[i]) swap(F[i], F[rev[i]]);
for (int len = 1; len <= lim; ++len) {
int half = 1 << len - 1;
for (int i = 0; i < 1 << lim; i += 1 << len) {
for (int j = 0; j < half; ++j) {
int F0 = F[i + j], F1 = mul(F[i + half + j], G[j << lim - len]);
F[i + j] = add(F0, F1); F[i + half + j] = sub(F0, F1);
}
}
}
}
int main() {
ios::sync_with_stdio(false); cin.tie(nullptr); cout.tie(nullptr);
cin >> n >> m;
for (int i = 0; i <= n; ++i) cin >> A[i];
for (int i = 0; i <= m; ++i) cin >> B[i];
init(n + m);
NTT(A); NTT(B);
for (int i = 0; i < 1 << lim; ++i) Mul(A[i], B[i]);
NTT(A, gi);
for (int i = 0, inv = qpow(1 << lim, P - 2); i <= n + m; ++i) cout << mul(A[i], inv) << ' ';
return 0;
}
毕竟前面 FFT 已经讲这么多了，也为 NTT 铺下了良好的基础，普通的多项式乘法这里就讲得简略了些，接下来才是重头戏。
多项式除法
前面讲到多项式除法的逆元，我们发现除法常常是不友好的，直接将两个多项式相除显然也不现实，所以可以考虑使用一些技巧。
回想取模意义下的除法，我们发现在模意义下不能直接进行除法，但是我们知道除法在模意义下是可以进行的，我们需要使用乘法逆元来解决这个问题，模意义下的一个数的乘法逆元很好求，那么现在的问题转化成要求一个多项式的乘法逆元。
如果在实数域里求乘法逆元，可以想到会出现 \(x^{-1}\) 这样的项，而这样显然对多项式是不友好的，就像分数对整数不太友好一样，所以考虑在取模意义下做，把分数转化为乘法逆元，现在就是要把 \(x^{-1}\) 转化成一个 \(x^a\) 满足 \(a\ge0\) 才好，所以可以想到模一个 \(x^n\) 这样的数，感觉上就挺可求的。
那么让我们再重新声明一下我们的问题，已知一个多项式 \(F(x)\)，求一个多项式 \(G(x)\) 满足：
\[F(x)G(x)\equiv1\pmod{x^n} \]这个 \(G(x)\) 能不能直接求呢？好像难度比较大，我们可以考虑使用倍增法。
倍增法是多项式中常用的一种思想，因为有一个不明原因的性质，就是很多情况下如果 \(a\le n\)，那么模 \(x^a\) 下的解是模 \(x^n\) 解的一部分，也就是说真实的解其实是一个无限长的式子，截取到哪里就能算到哪里。
首先设置倍增的起点，我们可以简单地发现在模 \(x\) 意义下的 \(F(x)\) 逆元就是 \(F(x)_0^{-1}\)，就是 \(F(x)_0\) 的乘法逆元，这个边界是非常明显的。
接下来考虑：
\[F(x)G'(x)\equiv1\pmod{x^n}\\ F(x)G(x)\equiv1\pmod{x^{2n}} \]已知 \(G'(x)\)，求 \(G(x)\)，那么接下来可以看推导过程：
\[\begin{aligned} F(x)G'(x)&\equiv1 &\pmod{x^n}\\ F(x)G(x)&\equiv1 &\pmod{x^{2n}}\\ F(x)(G(x)-G'(x))&\equiv0 &\pmod{x^n}\\ G(x)-G'(x)&\equiv0 &\pmod{x^n}\\ (G(x)-G'(x))^2&\equiv0 &\pmod{x^{2n}}\\ G(x)^2+G'(x)^2-2G(x)G'(x)&\equiv0 &\pmod{x^{2n}}\\ G(x)+F(x)G'(x)^2-2G'(x)&\equiv0 &\pmod{x^{2n}}\\ G(x)&\equiv2G'(x)-F(x)G'(x)^2&\pmod{x^{2n}}\\ \end{aligned} \]其中可能存在问题的主要是两边同时平方的部分，下面给出证明：
假设 \(F(x)\) 在模 \(x^n\) 意义下为 \(0\)，即满足 \(\forall i\in[0,n)\cap\mathbb N,\ F(x)_i=0\)，我们发现令 \(G(x)=F(x)^2\)，则 \(G(x)_i=\sum_{j=0}^{i}F(x)_jF(x)_{i-j}\)，其中因为 \(i<2n\)，所以 \(j\) 和 \(i-j\) 必然存在一个小于 \(n\)，所以在满足 \(\forall i\in[0,2n)\cap\mathbb N,\ G(x)_i=0\)。
这样就结束了，直接倍增，每次用上面最后一个式子算就好了，求模 \(x^n\) 意义下就倍增到 \(x^{2^k}\)，满足 \(2^k\ge n\) 就好了。
例题P4238 【模板】多项式乘法逆 代码：
多项式求逆
#include
using namespace std;
#define Reimu inline void // 灵梦赛高
#define Marisa inline int // 魔理沙赛高
#define Sanae inline bool // 早苗赛高
typedef long long LL;
typedef unsigned long long ULL;
typedef pair Pii;
typedef tuple Tiii;
#define fi first
#define se second
const int P = 998244353, G = 3, Gi = 332748118;
Marisa add(int x, int y) { return (x += y) < P ? x : x - P; }
Marisa&Add(int&x, int y) { return (x += y) < P ? x : x -= P; }
Marisa sub(int x, int y) { return (x -= y) < 0 ? x + P : x; }
Marisa&Sub(int&x, int y) { return (x -= y) < 0 ? x += P : x; }
Marisa mul(int x, int y) { return 1LL * x * y % P; }
Marisa&Mul(int&x, int y) { return x = 1LL * x * y % P; }
Marisa qpow(int x, int y) { int res = 1; for (; y; y >>= 1, Mul(x, x)) if (y & 1) Mul(res, x); return res; }
const int N = 400010;
int n, lim;
int A[N], B[N], A_[N], rev[N], g[N], gi[N];
Reimu init(int len) {
lim = __lg(len) + 1;
for (int i = 1; i < 1 << lim; ++i) rev[i] = rev[i >> 1] >> 1 | (i & 1) << lim - 1;
g[0] = gi[0] = 1; g[1] = qpow(G, P - 1 >> lim); gi[1] = qpow(Gi, P - 1 >> lim);
for (int i = 2; i < 1 << lim; ++i) g[i] = mul(g[i - 1], g[1]), gi[i] = mul(gi[i - 1], gi[1]);
}
Reimu NTT(int *F, int *G = g) {
for (int i = 1; i < 1 << lim; ++i) if (i < rev[i]) swap(F[i], F[rev[i]]);
for (int len = 1; len <= lim; ++len) {
int half = 1 << len - 1;
for (int i = 0; i < 1 << lim; i += 1 << len) {
for (int j = 0; j < half; ++j) {
int F0 = F[i + j], F1 = mul(F[i + half + j], G[j << lim - len]);
F[i + j] = add(F0, F1); F[i + half + j] = sub(F0, F1);
}
}
}
}
int main() {
ios::sync_with_stdio(false); cin.tie(nullptr); cout.tie(nullptr);
cin >> n;
for (int i = 0; i < n; ++i) cin >> A[i];
B[0] = qpow(A[0], P - 2);
for (int len = 1; len < n; len <<= 1) {
init((len << 2) - 1);
memcpy(A_, A, sizeof(int) * len << 1); memset(A_ + (len << 1), 0, sizeof(int) * ((1 << lim) - (len << 1)));
NTT(A_); NTT(B);
for (int i = 0; i < 1 << lim; ++i) B[i] = sub(add(B[i], B[i]), mul(A_[i], mul(B[i], B[i])));
NTT(B, gi);
for (int i = 0, inv = qpow(1 << lim, P - 2); i < len << 1; ++i) Mul(B[i], inv);
memset(B + (len << 1), 0, sizeof(int) * ((1 << lim) - (len << 1)));
}
for (int i = 0; i < n; ++i) cout << B[i] << ' ';
return 0;
}
要注意实现上不能每次都算一个长为 \(\mathcal O\left(n\right)\) 的多项式乘法，这样时间复杂度是 \(\mathcal O\left(n\log^2n\right)\) 的，时间复杂度就退化了，如果每次从 \(k\) 递推到为 \(2k\) 时只乘到 \(\mathcal O\left(k\right)\) 的话就是 \(\mathcal T\left(n\right)=\mathcal T\left(\dfrac n2\right)+\mathcal O\left(n\log n\right)\)，可以发现，这里 \(\mathcal T\left(n\right)\) 中的 \(n\) 取遍 \(1,2,\cdots 2^k\)，所以总和也只有 \(2n\) 左右，总时间复杂度 \(\mathcal O\left(n\log n\right)\)，当然也可以直接套主定理，时间复杂度算出来是一样的。
MTT
任意模数 NTT，也称为 MTT。一般 MTT 有两种实现方法，其实原理都比较简单，就是根本不管模数，算完再说，这个办法确实管用。
FFT 实现 MTT
首先想到如果说不管模数第一个想法就是上 FFT，这个办法可行，但是 FFT 是基于虚数的，想要用虚数就得有 double 这个爆精王，最多存 \(15\) 位，实际上 \(13\) 位以上就有点慌兮兮的，两个多项式系数全是 \(10^9\)，乘在一起肯定爆精度，long double 会好一些，但也不指望能好到哪里去，手写高精浮点数就是脑子被驴踢了，我们得想一种高明的办法。
发现 \(10^9\) 确实会爆，我们可以使用一些奇技淫巧，比如把数拆成两半，比如因为 \(2^{30}>10^9\)，所以直接把 \(a\) 拆成 \(\left\lfloor\dfrac a{2^{15}}\right\rfloor\cdot2^{15}+a\bmod2^{15}\)，同理我们可以啊多项式也给拆开，即 \(F(x)=2^{15}F'(x)+F''(x)\)，这样的话就可以得到：
\[F(x)G(x)=(2^{15}F'(x)+F''(x))(2^{15}G'(x)+G''(x))=2^{30}F'(x)G'(x)+2^{15}(F'(x)G''(x)+F''(x)G'(x))+F''(x)G''(x) \]这样就是 \(8\) 遍 FFT，然后最后再把答案求出来，取模就好了。
显然，这样常数太大了，我们需要一个更优的解法。
实际上，我们回想 FFT 的优化，发现这道题的可优化空间也并不是没有，因为虚部并没有用上，之前在做两个数乘法时的优化这里没有用，因为每个多项式都被用到了两次，那个等于白搭，所以我们考虑用另一种方式来利用虚部。
首先要介绍一个概念——共轭复数。