多项式卷积 FFT 和 NTT

FFT

多项式

定义

多项式是形如 $A(x)={\sum }_{i=0}^n{a}_i\times x^i$ 的表达式，其中 $x$ 为不定元，这是多项式的系数表示法。

将多项式中不定元的最大次数称为多项式的次数，记做 $\partial (A(x))$。

运算

$$\begin{array}{c}A(x)={\sum }_{i=0}^n{a}_i{x}^i,B(x)={\sum }_{i=0}^m{b}_i{x}^i\\ \\ A(x)\pm B(x)={\sum }_{i=0}^{\max \{n,m\}}\left(a_i\pm b_i\right)x^i\\ \\ A(x)\times B(x)={\sum }_{i=0}^n{\sum }_{j=0}^m\left(a_i\times b_j\right)x^{i+j}\\ ={\sum }_{k=0}^{n+m}{\sum }_{i=0}^n\left(a_i\times b_{k-i}\right)x^k\end{array}$$

多项式乘法满足交换律，结合律和分配率。

多项式的除法通常指带余除法。

$$\begin{gathered} A(x)=q(x)B(x)+r(x)\partial (r(x))<\partial (B(x)) \\ A(x)\equiv r(x)(\mathrm{m}\mathrm{o}\mathrm{d}B)(x) \end{gathered}$$

$q(x)$ 称为 $B(x)$ 除 $A(x)$ 的商式，$r(x)$ 称为 $B(x)$ 除 $A(x)$ 的余式。若 $r(x)=0$，则称 $B(x)$ 整除 $A(x)$，记作 $B(x)|A(x)$。

点值表示法

$n$ 个点可以确定唯一的 $n-1$ 次多项式。如果多项式 $f(x)$ 和 $g(x)$ 的次数不超过 $n$，且对于 $n+1$ 个不同的数都有相同的对应值，则 $f(x)=g(x)$。

如果选取 $n+1$ 个不同的数，$x_0\sim x_n$ 对多项式求值，得到 $A(x_0),A(x_1),\cdots ,A(x_n)$。那么称 $\{(x_i,A(x_i)|0\le i\le n,i\in \mathbb{N}\}$ 为多项式 $A(x)$ 的点值表示。通过系数表示求点值表示的过程叫做求值，通过点值表示求系数表示的过程叫做插值。

设 $A(x)$ 和 $B(x)$ 是两个多项式，它们在 $x_0\sim s_n$ 上的取值分别为 $c_0\sim c_n$ 和 $D_0\sim D_n$。那么 $A(x)\times B(x)$ 在 $x_0,x_1,\dots ,x_n$ 上的取值分别为 $c_0{D}_0\sim c_n{D}_n$。

复数

• 定义

虚数单位 $i=\sqrt{-1}$。

形如 $a+b_i(a,b\in \mathbb{R})$ 的数称为复数，$a$ 称为实部，$b$ 称为虚部。全体复数的集合记为 $\mathbb{C}$。

• 运算

设 $z=a+bi,w=c+di$，有

$$\begin{gathered} z\pm w=(a\pm c)+(b\pm d)i \\ z\times w=ac+adi+bci+bd{i}^2=(ac-bd)+(ad+bc)i \\ z÷w=\frac{a+bi}{c+di}=\frac{(a+bi)(c-di)}{(c+di)(c-di)}=\frac{(ac+bd)+(bc-ad)i}{c^2+d^2} \end{gathered}$$

复数的相加满足平行四边形法则。

共轭复数：$\overline{z}=a-bi$，相当于对虚部取反。

• 几何意义

令横轴为实轴，纵轴为虚轴。复数 $a+bi$ 就对应坐标为 $(a,b)$ 的点或向量。

复数对应的点到原点的距离称为复数的模长，记为 $r=|z|$。 由勾股定理 $|z|=\sqrt{a^2+b^2}=\sqrt{z\overline{z}}$

复数与实轴正半轴形成的角称为复数的辐角，记为 $\phi =\arg z$。

复数也可以用模长和辐角表示，称为复数的极坐标形式。由极坐标与直角坐标的转化关系，$z=r(\cos \phi +i\sin \phi )$。

复数的乘法，按照模长相乘，辐角相加的运算规则。除法为乘法的逆运算。

复数的共轭，则为复数对应向量与实轴轴对称的结果。

单位根

欧拉公式

对于任意实数 $x$，都有

$$e^{ix}=\cos x+i\sin x$$

证明略。这个公式可以说明当 $x$ 为实数时，函数 $f(x)=e^{ix}$ 可以在复数平面表述一个单位元，且 $x$ 为平面上的一个幅角。

定义

将方程 $z^n=1$ 的所有复数根称为 $n$ 次单位根，这样的单位根共有 $n$ 个。

$$e^{\frac{2\pi ki}{n}}(k=0,1,\cdots ,n-1)$$

这些单位根分布在复数平面的单位圆上，构成了一个正 $n$ 边形，它们把单位圆等分成 $n$ 个部分。

简单的，定义 $\omega =e^{\frac{2\pi ki}{n}}$，若有 $n$ 与 $k$ 互质，则 $\omega$ 称为 $n$ 次本原单位根。$w^x(x=0,1,\cdots ,n-1)$ 构成了所有的 $n$ 次单位根。特殊的，当 $k=1$ 时，本原单位根记为 ${\omega }_n$。

性质

$$\begin{gathered} w_n^n=1 \\ {\omega }_n^0,{\omega }_n^1,\cdots ,{\omega }_n^{n-1}互不相同 \\ {\omega }_n^k=\cos \frac{2\pi k}{n}+i\sin \frac{2\pi k}{n} \\ {\omega }_n^k={\omega }_{2n}^{2k} \\ {\omega }_n^{\frac{n}{2}}=e^{i\pi }=-1 \\ {\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k \\ 单位根的共轭复数为其倒数 \end{gathered}$$

Cooley-Tukey 算法

离散傅里叶变换 DFT

多项式有两个表示方法——系数表示法和点值表示法。它们各有千秋，系数表示法适用范围广，但乘法操作为 $O(n^2)$；点值表示法适用范围窄，但乘法操作为线性。那么能不能在两个表示法相互转换呢？

假设现在有一个 $n-1$ 次多项式，$A(x)={\sum }_{i=0}^n{a}_i\times x_i$。为了方便，设 $n=2^m,m\in \mathbb{N}$，不足在高位补 $0$。

将 $n$ 个 $n$ 次单位根 ${\omega }_n^0,{\omega }_n^1,\dots ,{\omega }_n^{n-1}$ 代入 $A(x)$ 将其转换成点值表达 $A\left({\omega }_n^k\right)={\sum }_{i=0}^{n-1}{a}_i{\omega }^{ki}$。 点值向量 $\vec{y}=\left(A\left({\omega }_n^0\right),A\left({\omega }_n^1\right),\dots ,A\left({\omega }_n^{n-1}\right)\right)$ 称为系数向量 $\vec{a}=\left(a_0,a_1,\dots ,a_{n-1}\right)$ 的离散傅里叶变换，记为 $\vec{y}={\mathrm{DFT}}_n(\vec{a})$。

可以发现离散傅里叶变换的复杂度为 $O(n^2)$ 的。Cooley-Tukey 算法对 DFT 进行了优化，它将每一项进行了奇偶分类。当 $k<\frac{n}{2}$ 时

根据性质 ${\omega }_n^k={\omega }_{2n}^{2k}$

$$\begin{gathered} A\left({\omega }_n^k\right)=\sum _{i=0}^{n-1}{a}_i{\omega }_n^{ki} \\ =\sum _{i=0}^{n-1}{a}_{2i}{\omega }_n^{2ki}+{\omega }_n^k\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{\omega }_n^{2ki} \\ =\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i}{\omega }_{\frac{n}{2}}^{ki}+{\omega }_n^k\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{\omega }_{\frac{n}{2}}^{ki} \end{gathered}$$
根据性质 $w_n^n=1$ 和 ${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$

$$\begin{array}{rl}A\left({\omega }_n^{k+\frac{n}{2}}\right)& =\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i}{\omega }_{\frac{n}{2}}^{ki+\frac{n}{2}\times i}+{\omega }_n^{k+\frac{n}{2}}\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{\omega }_{\frac{n}{2}}^{ki}\\ & =\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i}{\omega }_{\frac{n}{2}}^{ki}-{\omega }_n^k\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{\omega }_{\frac{n}{2}}^{ki}\end{array}$$

可以发现左右两个式子相似。令左式的系数表达式为 $A_1$，右式的系数表达式为 $A_2$。

$$\begin{gathered} A\left({\omega }_n^k\right)=A_1({\omega }_{\frac{n}{2}}^k)+A_2({\omega }_{\frac{n}{2}}^k) \\ A\left({\omega }_n^{k+\frac{n}{2}}\right)=A_1({\omega }_{\frac{n}{2}}^k)-A_2({\omega }_{\frac{n}{2}}^k) \end{gathered}$$

这意味着，如果能求出 $A_1({\omega }_{\frac{n}{2}}^k)$ 和 $A_2({\omega }_{\frac{n}{2}}^k)$，那么可以直接求出 $A\left({\omega }_n^k\right)$ 和 $A\left({\omega }_n^{k+\frac{n}{2}}\right)$。

于是，我们只要算出 $a_0,a_2,\dots ,a_{\frac{n}{2}}$ 和 $a_1,a_3,\dots ,a_{\frac{n}{2}-1}$ 的 DFT 就可以算出 $a_0,a_1,\dots a_n$ 的 DFT，把问题规模缩小了一半，递归下去即可。时间复杂度 $O\left(n{\log }_{2}n\right)$。

傅里叶逆变换 IDFT

刚刚计算的是 $\vec{y}={\mathrm{DFT}}_n(\vec{a}),$ 可以将多项式转化成点值表示，现在为了将点值表示转化成系数表示，需要计算 IDFT，它是 DFT 的逆运算，也就是插值的过程。提起插值，可能先想起拉格朗日插值，不过 FFT 不需要它。

考虑解一个方程组

$$\{\begin{array}{ccccc}{a}_0{\left({\omega }_n^0\right)}^0& +\cdots & +a_{n-2}{\left({\omega }_n^0\right)}^{n-2}& +a_{n-1}{\left({\omega }_n^0\right)}^{n-1}& =A\left({\omega }_n^0\right)\\ a_0{\left({\omega }_n^1\right)}^0& +\cdots & +a_{n-2}{\left({\omega }_n^1\right)}^{n-2}& +a_{n-1}{\left({\omega }_n^1\right)}^{n-1}& =A\left({\omega }_n^1\right)\\ ⋮& ⋮& ⋮& ⋮& \\ a_0{\left({\omega }_n^{n-1}\right)}^0& +\cdots & +a_{n-2}{\left({\omega }_n^{n-1}\right)}^{n-2}& +a_{n-1}{\left({\omega }_n^{n-1}\right)}^{n-1}& =A\left({\omega }_n^{n-1}\right)\end{array}$$

写成矩阵的形式

$$\left[\begin{array}{cccc}{\left({\omega }_n^0\right)}^0& {\left({\omega }_n^0\right)}^1& \cdots & {\left({\omega }_n^0\right)}^{n-1}\\ {\left({\omega }_n^1\right)}^0& {\left({\omega }_n^1\right)}^1& \cdots & {\left({\omega }_n^1\right)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ {\left({\omega }_n^{n-1}\right)}^0& {\left({\omega }_n^{n-1}\right)}^1& \cdots & {\left({\omega }_n^{n-1}\right)}^{n-1}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}A\left({\omega }_n^0\right)\\ A\left({\omega }_n^1\right)\\ ⋮\\ A\left({\omega }_n^{n-1}\right)\end{array}\right]$$

左矩阵即为系数矩阵 $V$，要求它的逆矩阵 $V^{-1}$。可以发现，这个矩阵可以使用范德蒙德矩阵的逆矩阵公式，令 $D$ 为 $V$ 的伴随矩阵。

$$V^{-1}=\frac{1}{n}D=\frac{1}{n}\left[\begin{array}{cccc}{\left({\omega }_n^{-0}\right)}^0& {\left({\omega }_n^{-0}\right)}^1& \cdots & {\left({\omega }_n^{-0}\right)}^{n-1}\\ {\left({\omega }_n^{-1}\right)}^0& {\left({\omega }_n^{-1}\right)}^1& \cdots & {\left({\omega }_n^{-1}\right)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ {\left({\omega }_n^{-(n-1)}\right)}^0& {\left({\omega }_n^{-(n-1)}\right)}^1& \cdots & {\left({\omega }_n^{-(n-1)}\right)}^{n-1}\end{array}\right]$$

证明较为复杂，我们简单验证一下，设 $E=V\times D$，则有

$$E_{i,j}=\sum _{k=0}^{n-1}{D}_{ik}{V}_{kj}=\sum _{k=0}^{n-1}{\omega }_n^{-ik}{\omega }_n^{kj}=\sum _{k=0}^{n-1}{\omega }_n^{k(j-i)}$$

当 $i=j$ 时，显然 $E_{i,j}=n$。当 $i\ne j$ 时，有

$$E_{i,j}=\sum _{k=0}^{n-1}{\left({\omega }_n^{j-i}\right)}^k=\frac{1-{\left({\omega }_n^{j-i}\right)}^n}{1-{\omega }_n^{j-i}}=0$$

因此 $D$ 的主对角线都是 $n$，其它位置都是 $0$，除以 $n$ 后即为单位矩阵，证毕。

综上所述，可以得到

$$\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}\right]=\frac{1}{n}\left[\begin{array}{cccc}{\left({\omega }_n^{-0}\right)}^0& {\left({\omega }_n^{-0}\right)}^1& \cdots & {\left({\omega }_n^{-0}\right)}^{n-1}\\ {\left({\omega }_n^{-1}\right)}^0& {\left({\omega }_n^{-1}\right)}^1& \cdots & {\left({\omega }_n^{-1}\right)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ -(n-1)& {\left({\omega }_n^{-(n-1)}\right)}^1& \cdots & {\left({\omega }_n^{-(n-1)}\right)}^{n-1}\end{array}\right]\left[\begin{array}{c}A\left({\omega }_n^0\right)\\ A\left({\omega }_n^1\right)\\ ⋮\\ A\left({\omega }_n^{n-1}\right)\end{array}\right]$$

IDFT 就相当于把 DFT 过程中的 $w_n^i$ 换成 $w_n^{-i}$，即为共轭复数。然后做一次 DFT，将结果除以 $n$ 即可。

递归实现

#include <bits/stdc++.h>
using namespace std;
const int N = 5e6 + 5, INF = 0x3f3f3f3f; //一定要把 N 开到 5 倍 
inline int read() {
	int x = 0, f = 0; char ch = 0;
	while (!isdigit(ch)) f |= ch == '-', ch = getchar();
	while (isdigit(ch)) x = (x << 3) + (x << 1) + (ch ^ 48), ch = getchar();
	return f ? -x : x;
}
typedef complex <double> E;
const double pi = acos(-1); //π
void fft(E *a, int mid, int flg) { //单位根共轭的 flg 
	if (mid == 0) return ;
	E a1[mid], a2[mid]; //奇偶分类 
	for (int i = 0; i < mid; i++) a1[i] = a[i << 1], a2[i] = a[i << 1 | 1]; //必须 < mid 
	fft(a1, mid >> 1, flg); fft(a2, mid >> 1, flg);  //递归调用数组次数太多
	for (int i = 0; i < mid; i++) {
		E w(cos(pi * i / mid), flg * sin(pi * i / mid)); //单位根的 2 可以被 n / 2 = mid 约去 
		a[i] = a1[i] + w * a2[i], a[i + mid] = a1[i] - w * a2[i];
	} 
}
E a[N], b[N]; 
int main() {
	int n = read(), m = read();
	for (int i = 0; i <= n; i++) a[i] = read(); //实部 
	for (int i = 0; i <= m; i++) b[i] = read();
	for (m += n, n = 1; n <= m; n <<= 1); // 把长度补到 2 的幂
	fft(a, n >> 1, 1); fft(b, n >> 1, 1); //DFT
	for (int i = 0; i < n; i++) a[i] *= b[i]; //n 个点值
	fft(a, n >> 1, -1); //IDFT
	for (int i = 0; i <= m; i++) { //循环到 m 
		printf("%.0f ", fabs(a[i].real()) / n); //除以 n 要小心 -0 
	}
	return 0;
}

迭代实现

观察递归的观察，可以发现，最后的运算顺序，是按二进制位反转之后从小到大排序的。因此我们只需要实现一个二进制反转然后排序，这样差为 $\frac{n}{2}$ 的元素就在一组了。

具体的，模拟从高位到低位的二进制加法。维护当前下标 $i$ 和反相加的下标 $j$，表示反转后下标 $i$ 应放反转前下标 $j$ 的东西。如果 $i>j$，就将 $i$ 和 $j$ 存的东西交换，这样就能得到所需的数组。

#include <bits/stdc++.h>
using namespace std;
const int N = 5e6 + 5, INF = 0x3f3f3f3f; //一定要把 N 开到 5 倍 
inline int read() {
	int x = 0, f = 0; char ch = 0;
	while (!isdigit(ch)) f |= ch == '-', ch = getchar();
	while (isdigit(ch)) x = (x << 3) + (x << 1) + (ch ^ 48), ch = getchar();
	return f ? -x : x;
}
typedef complex <double> E; //手写会更快
const double pi = acos(-1); //π
void fft(E *a, int n, int flg) {
	for (int i = 0, j = 0; i < n; i++) {
		if (i < j) swap(a[i], a[j]);
		for (int k = n >> 1; (j ^= k) < k; k >>= 1);
	} 
	for (int i = 2; i <= n; i <<= 1)  //待合并的区间长度
		for (int j = 0, mid = (i >> 1); j < n; j += i) {
			E wn(cos(pi / mid), flg * sin(pi / mid)) /*本原单位根*/, w(1, 0) /*幂*/; 
			for (int k = 0; k < mid; k++, w *= wn) {
				//E w(cos(pi * k / mid), flg * sin(pi * k / mid)); n log n 次 cos 和 sin 也会超时 
				E x = a[j + k + mid] * w;
				a[j + k + mid] = a[j + k] - x, a[j + k] = a[j + k] + x; //合并的顺序一定不能换 
			}	
		}		
}
E a[N], b[N]; 
int main() {
	int n = read(), m = read();
	for (int i = 0; i <= n; i++) a[i] = read(); //实部 
	for (int i = 0; i <= m; i++) b[i] = read();
	for (m += n, n = 1; n <= m; n <<= 1); // 把长度补到 2 的幂
	fft(a, n, 1); fft(b, n, 1); //DFT
	for (int i = 0; i < n; i++) a[i] *= b[i]; //n 个点值
	fft(a, n, -1); //IDFT
	for (int i = 0; i <= m; i++) { //循环到 m 
		printf("%.0f ", fabs(a[i].real()) / n); //除以 n 要小心 -0 
	}
	return 0;
}

NTT

原根

review 一下 FFT 中单位根的一些性质

$$\begin{gathered} w_n^n=1 \\ {\omega }_n^0,{\omega }_n^1,\cdots ,{\omega }_n^{n-1}互不相同 \\ {\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k \end{gathered}$$

根据费马小定理，任意一个质数 $p$，有

$$a^{p-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$$

定义 $p$ 的原根 $g$ 为使得 $g^0,g^1,\cdots ,g^{p-2}\mathrm{m}\mathrm{o}\mathrm{d}p$ 互不相同的数。

如果一个质数 $p$ 能表示为 $p=k\times 2^n+1$，并找到它的原根 $g$。令 $g_n=g^k\mathrm{m}\mathrm{o}\mathrm{d}p$，这样可以使得$g_n^0,g_n^1,\cdots ,g_n^{p-2}\mathrm{m}\mathrm{o}\mathrm{d}p$ 互不相同，且有 $g_n^n\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$。又因为 $g_n^k$ 互不相同，所以 $g_n^{\frac{n}{2}}\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}p)$

所以在 NTT 中，可以用原根代替单位根，精度更高常数更小。不足之处是对模数有限制。这一点可以使用任意模数 NTT，因为过于毒瘤，本文不作讲解。

$$\begin{array}{l}998244353=119\cdot 2^{23}+1,g=3\\ 1004535809=479\cdot 2^{21}+1,g=3\end{array}$$

实现