快速傅里叶/数论变换学习简记

多项式部分是 OI 中一个比较 math 和 naive 的部分，各种多项式题声(chou)名(ming)远(zhao)扬(zhu)。
因为作者数学较弱，因此这篇文章会比较注重数学部分的理解。
同时感谢所有提供借鉴的博文和帮助作者的人！

接下来进入正文

1.多项式乘法

现在我们有两个多项式 $F(x),G(x)$，次数分别为 $n,m$，将两个多项式相乘，得到一个新的多项式 $H(x)$，新的多项式应该是 $n+m$ 次的。记多项式 $F(x)$ 的 $n$ 次项系数为 $[x^n]F(x)$。朴素地算 $H(x)$ ，

$$[x^i]H(x)=\sum _{j+k=i}[x^j]F(x)\cdot [x^k]G(x)$$

时间复杂度为 $\mathcal{O}(nm)$。
这是由系数表示法引出来的计算方法。
这种方法由于每个点的系数固定，没有什么想法优化。
考虑初中的时候，我们用三个点的坐标求二次函数的表达式，受此启发，我们用 $n$ 个点来表示 $n-1$ 次多项式。将 $n$ 个互不相同的 $x$ 带入多项式，会得到 $n$ 个不同的取值 $y$。
该多项式被 $n$ 个点唯一确定。
其中 $y_i=\sum _{j=1}^n{a}_j\cdot x_i^j$。
当我们知道 $n$ 个点和 $m$ 个点表示出来的多项式点值时，我们做乘法只需要 $\mathcal{O}(n+m)$ 的时间复杂度（因为要 $n+m+1$ 个点才能确定 $n+m$ 次多项式）。
我们发现这简直太快了，因此我们需要找到一个将系数转换为点值，再重新转换为系数的算法。
这就是 FFT。
FFT 分为两个部分：a.系数转点值（DFT）b.点值转系数（IDFT）
在介绍 FFT 前，我们先介绍以下复数。

2.复数

DFT 的关键之处在于，只要点值足够，代入是可以自己决定的。
我们在实数范围内似乎找不到 $n$ 个有很好性质的点让我们将时间复杂度降低。于是我们将其与数学的另一个毒瘤复数结合。
接下来是关于复数部分的论述。
前置知识：

• 圆的弧度制
• 向量的运算法则
  引出复数定义：设 $a,b\in \mathbb{R},i^2=-1$，形如 $a+bi$ 的数称为复数。其中 $i$ 为虚数单位。
  复平面直角坐标系：$x$ 轴表示实数，$y$ 轴（除原点）代表虚数，$(a,b)$ 表示 $a+bi$。
  模长： $\sqrt{a^2+b^2}$。
  幅角：以逆时针为正方向，从 $x$ 轴正半轴到已知向量的转角的有向角。
  共轭复数：实部相同，虚部相反。
  运算法则：
  加法：在复平面中，复数可以表示为向量，因此等同于向量。
  乘法：模长相乘，幅角相加。

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

除法：分子分母同时乘分母的共轭复数。

$$\frac{a+bi}{c+di}=\frac{(a+bi)(c-di)}{(c+di)(c-di)}=\frac{ac+bd+(bc-ad)i}{c^2+d^2}=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i$$

单位根：
在复平面上，以原点为中心，$1$ 为半径作圆，所作的圆为单位圆。
以原点为起点，圆的 $n$ 等分点为终点做 $n$ 个向量，从原点开始依次标号为 ${\omega }_n^0,{\omega }_n^1,w_n^2,\dots ,{\omega }_n^{n-1},{\omega }_n^n$，其中 ${\omega }_n^0={\omega }_n^n=1$，对于 $k\ge n,{\omega }_n^k={\omega }_n^{k\mathrm{\%}n}$。
单位根的定义非常重要 ${\omega }_n^k=({\omega }_n^1{)}^k$。
这个公式直接导致了下文我们选择求单位根的方法。

3.多项式乘法加速版本

于是先得出一个结论：单位根 ${\omega }_n^0\sim {\omega }_n^{n-1}$ 有不错的性质，可以分治多项式乘法降低时间复杂度到 $\mathcal{O}(n\log n)$。
以下是推导：
先假设 $n=2^k,k\in {\mathbb{Z}}_{\mathbb{+}}$。
设多项式 $F(x)$ 的系数为 $f_0,f_1,f_2,\dots ,f_{n-1}$。
按下标奇偶性分类：

$$F(x)=(f_0+f_2\cdot x^2+f_4\cdot x^4+\cdots +f_{n-2}\cdot x^{n-2})+(f_1\cdot x+f_3\cdot x^3+f_5\cdot x^5+\cdots +f_{n-1}\cdot x^{n-1})$$

设

$$F_1(x)=f_0+f_2\cdot x+f_4\cdot x^2+\cdots +a_{n-2}\cdot x^{\frac{n}{2}-1}$$

$$F_2(x)=f_1+f_3\cdot x+f_5\cdot x^2+\cdots +f_{n-1}\cdot x^{\frac{n}{2}-1}$$

将两式用 $x^2$ 进行一个代换？
$F_1(x^2)+F_2(x^2)$ 好像看不出什么规律。
给 $F_2(x^2)$ 乘个 $x$？
似乎 $F(x)=F_1(x^2)+x{F}_2(x^2)$。
$\mathrm{\forall }k\in [0,\frac{n}{2})$，将 ${\omega }_n^k$ 代入得：

$$\begin{array}{rl}F({\omega }_n^k)& =F_1({\omega }_n^{2k})+{\omega }_n^k\cdot F_2({\omega }_n^{2k})\\ & =F_1({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k\cdot F_2({\omega }_{\frac{n}{2}}^k)\end{array}$$

分治惯用套路，将 ${\omega }_n^{k+\frac{n}{2}}$ 代入得：

$$\begin{array}{rl}F({\omega }_n^{k+\frac{n}{2}})& =F_1({\omega }_n^{2k+n})+{\omega }_n^{k+\frac{n}{2}}\cdot F_2({\omega }_n^{2k+n})\\ & =F_1({\omega }_n^{2k}\cdot {\omega }_n^n)-{\omega }_n^k\cdot F_2{\omega }_n^{2k}\cdot {\omega }_n^n\\ & =F_1({\omega }_n^{2k})-{\omega }_n^k\cdot F_2({\omega }_n^{2k})\\ & =F_1({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k\cdot F_2({\omega }_{\frac{n}{2}}^k)\end{array}$$

我们欣喜地发现，$F({\omega }_n^k)$ 和 $F({\omega }_n^{k+\frac{n}{2}})$ 之间只有一个正负号的区别，因此我们在处理第 $k$ 项时，第 $k+\frac{n}{2}$ 项也可以同时处理出来。
也就是说，当我们知道 $F_1(x)$ 和 $F_2(x)$ 分别在 ${\omega }_{\frac{n}{2}}^i,i\in [0,\frac{n}{2}-1]$ 的点值时，我们就可以 $\mathcal{O}(n)$ 求出 $F(x)$ 在 ${\omega }_n^i,i\in [0,n-1]$ 的点值表示。
于是问题的规模缩小了一半。
接下来的问题在于怎么快速求出 $F_1(x)$ 和 $F_2(x)$ 的点值表示。
嘶~似乎 $F_1(x)$ 和 $F_2(x)$ 的性质与 $F(x)$ 一样？
那我们继续递归？
因为 $n$ 是 $2$ 的正整数次幂，所以不存在会分的不均匀，所以递归正确。
实际情况下，一般通过高次补 $0$ 的办法补齐。
至此我们证明了 FFT 的实现是如同线段树一样的递归实现，时间复杂度 $\mathcal{O}(n\log n)$。
现在我们可以开始写代码了。

4.FFT的代码实现

4.1 复数的四则运算

STL 中封装了 complex 类模板，但实现较为复杂，常数巨大。
因此考虑手写。

struct complex {
    double x, y;
    complex () {}
    complex (double _x, double _y) : x(_x), y(_y) {}
    complex operator+ (const complex& _r) noexcept { return complex(x + _r.x, y + _r.y); }
    complex operator- (const complex& _r) noexcept { return complex(x - _r.x, y - _r.y); }
    complex operator* (const complex& _r) noexcept { 
	return complex(x * _r.x - y * _r.y, x * _r.y + y * _r.x); 
    }
    complex operator/ (const complex& _r) noexcept {
	double t = _r.x * _r.x + _r.y * _r.y;
	return complex((x * _r.x + y * _r.y) / t, (y * _r.x - x * _r.y) / t);
    }
};

均由上述复数运算法则推导。

4.2 预处理单位根

在上文我们已经讲出了 ${\omega }_n^k$ 与 ${\omega }_n^1$ 的关系，因此我们只需要求出 ${\omega }_n^1$ 即可。
${\omega }_n^1$ 的值用高中数学的知识就可求出 $(\cos (\frac{2\pi }{n}),\sin (\frac{2\pi }{n}))$。

complex unit(cos(2*Pi/n), sin(2*Pi/n)), cur(1, 0);
vector<complex> w(n);
for (int i = 0; i < n; ++i) {
    w[i] = cur; cur = cur * unit;
}

4.3 常数超大的写法

我们在上文用的是递归方式思考的 FFT，那么递归自然是一种写法

void dft(complex *f, int n) { // dft -> 系数转点值
    if (n == 1) return;
    complex *f1 = f, *f2 = f + (n >> 1);
    vector<complex> g(n);
    // 指针指向的是数组头
    for (int k = 0; k < n; ++k) g[k] = f[k];
	for (int k = 0; k < (n >> 1); ++k) 
	    f1[k] = g[k << 1], f2[k]  = g[k << 1 | 1];
    fft(f1, n >> 1); fft(f2, n >> 1);
    // 递归求解
    // 每次的单位根都是重新对圆进行划分，所以无法预处理
    complex unit(cos(2 * Pi / n), sin(2 * Pi / n)), cur(1, 0);
    for (int k = 0; k < (n >> 1); ++k) {
	g[k] = f1[k] + cur * f2[k];
	g[k + (n >> 1)] = f1[k] - cur * f2[k];
	cur = cur * unit;
    }
    for (int i = 0; i < n; ++i) f[i] = g[i];
}

接下来的问题是，目前我们得到的是点值表示，但我们要求的是系数，所以要用 IDFT。

4.4 IDFT

将 DFT 中 ${\omega }_n^1$ 替换为 ${\omega }_n^{-1}$，最后除以 $n$。
这个的证明很有价值，可以推出单位根反演。
我们现在已经知道了点值的各种性质，而 DFT 是求点值，因此我们只需要知道怎么将点值还原即可。
设 $F(x)$ 的点值数列为 $G$。
即 $G=DFT(F)$。
回忆我们如何求点值 $g_k=\sum _{i=0}^{n-1}({\omega }_n^k{)}^i\cdot f_i$
先给出结论： $n\cdot f_k=\sum _{i=0}^{n-1}({\omega }_n^{-k}{)}^i\cdot g_i$
这个式子可以将点值转换为系数。

证明：
右边 $=\sum _{i=0}^{n-1}({\omega }_n^{-k}{)}^i\cdot g_i$
将 $g_i$ 代入可得：
$=\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{\omega }_n^{ij}{\omega }_n^{-ik}\cdot f_j$
$=\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{\omega }_n^{i(j-k)}\cdot f_j$
对 $j,k$ 的关系进行分类讨论：

1. $j=k$ 贡献 $\sum _{i=0}^{n-1}{\omega }_n^0\cdot f_k=n\cdot f_k$
2. $j\ne k$，假设 $p=j-k$，贡献 $\sum _{i=0}^{n-1}{\omega }_n^{ip}{f}_{k+p}={\omega }_n^p(\sum _{i=0}^{n-1}{\omega }_n^i)f_{k+p}$
   由等比数列求和公式：$\sum _{i=0}^{n-1}{\omega }_n^i={\omega }_n^0\frac{{\omega }_n^n-1}{{\omega }_n^1-1}=0$，这一部分没有贡献。

于是我们就证明了结论。
再看这个式子，发现其实和我们求点值的表达式非常像，因此相当于将 $G$ 数列当作系数再求一遍点值。
但是我们需要求的东西变成了 ${\omega }_n^{-i},i\in [1,n-1]$。
然后 ${\omega }_n^{-1}=(\cos (\frac{2\pi }{n}),-\sin (\frac{2\pi }{n}))$。
这一段和 DFT 长的特别像。

void idft(Complex *f, int n) {
    if (n == 1) return;
    Complex *f1 = f, *f2 = f + (n >> 1);
    vector<Complex> g(n);
    for (int k = 0; k < n; ++k) g[k] = f[k];
    for (int k = 0; k < (n >> 1); ++k) 
	f1[k] = g[k << 1], f2[k] = g[k << 1 | 1];
    idft(f1, n >> 1); idft(f2, n >> 1);
    Complex unit(cos(2 * Pi / n), -sin(2 * Pi / n)), cur(1, 0);
    // 和上文相比，我们只有这里一个符号有区别，这为我们化简代码提供了思路。
    for (int k = 0; k < (n >> 1); ++k) {
	g[k] = f1[k] + cur * f2[k];
	g[k + (n >> 1)] = f1[k] - cur * f2[k];
	cur = cur * unit;
    }
    for (int i = 0; i < n; ++i) f[i] = g[i];
}

但是我们都知道，递归和数组复制的常数超级大，我们需要一个常数更小的办法。

4.5 基于二进制的优化策略

我们对分治的数组下标进行如下考虑。

这是原来的数组
0 1 2 3 4 5 6 7
变换一次后
0 2 4 6|1 3 5 7
变换两次后
0 4|2 6|1 5|3 7
变换三次后
0|4|2|6|1|5|3|7

我们将最后的数列和初始数组下标转换成二进制。

000 001 010 011 100 101 110 111
000 100 010 110 001 101 011 111

嗯？好像最后每个位置上的数是原位置的二进制翻转？
但是怎么快速求？常规办法 $\mathcal{O}(n\log n)$ 有点慢，不适合卡常。
那我们尝试一下同样使用与递推有关的方式解决？
考虑一个类似 DP 的转移，当前二进制数的翻转等于将最后一位数提取出来放在开头，和剩下的数的翻转合并在一起。
考虑 r[i] 表示 $i$ 的二进制下翻转数。

for (int i = 0;i < n; ++i)
    r[i] = (r[i >> 1] >> 1) | ((i & 1) ? n >> 1 : 0);

这里有一个细节，我们在处理 i>>1 的时候人为地为 r[i >> 1] 后补了 $0$，因为我们要将 i>>1 补全位数，此时我们要剔除这个 $0$ 的干扰。因此要左移处理。
然后我们可以采用某种迭代的方式避免递归。
所以我们拥有了一个飞快的 FFT。

#include <bits/stdc++.h>

using namespace std;

const int N = 1.5e6 + 5;

struct Complex {
	double x, y;
	Complex () {}
	Complex (double _x, double _y) : x(_x), y(_y) {}
	Complex operator+ (const Complex& _r) noexcept { return Complex(x + _r.x, y + _r.y); }
	Complex operator- (const Complex& _r) noexcept { return Complex(x - _r.x, y - _r.y); }
	Complex operator* (const Complex& _r) noexcept { 
		return Complex(x * _r.x - y * _r.y, x * _r.y + y * _r.x); 
	}
	Complex operator/ (const Complex& _r) noexcept {
		double t = _r.x * _r.x + _r.y * _r.y;
		return Complex((x * _r.x + y * _r.y) / t, (y * _r.x - x * _r.y) / t);
	}
}f[N << 1], g[N << 1];

const double Pi = acos(-1.0);

int n, m, r[N << 1], lim = 1;

void fft(Complex* f, int flag) {
    for (int i = 0; i < lim; ++i) if (i < r[i]) swap(f[i], f[r[i]]);
    for (int p = 2; p <= lim; p <<= 1) {
    // 按照递归，我们当前从小区间合并到大区间
        int len = p >> 1;
        Complex unit(cos(2 * Pi / p), flag * sin(2 * Pi / p));
        for (int k = 0; k < lim; k += p) {
            // 当前区间为 k ~ k + p - 1
            // 左半部分为 k ~ k + len - 1，右半部分为 k + len ~ k + p - 1
            Complex cur(1, 0);
	    for (int l = k; l < k + len; ++l) {
	        Complex t = cur * f[len + l];
	        f[l + len] = f[l] - t;
	        f[l] = f[l] + t;
	        cur = cur * unit;
            }
	}
    }
}

int main() {
    ios::sync_with_stdio(false);
    cin >> n >> m;
    for (int i = 0; i <= n; ++i) cin >> f[i].x;
    for (int i = 0; i <= m; ++i) cin >> g[i].x;
    while (lim <= (n + m)) lim <<= 1;
    for (int i = 0; i < lim; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) ? lim >> 1 : 0);
    fft(f, 1); fft(g, 1);
    for (int i = 0; i < lim; ++i) f[i] = f[i] * g[i];
    fft(f, -1);
    for (int i = 0; i <= n + m; ++i) cout << (int)(f[i].x / lim + 0.5) << " ";
    return 0;
}

至此，FFT的部分结束。多项式的基础部分代码不是很繁琐。

---

以上是快速傅里叶变换内容，但是复数的精度收到限制，且不好取模。
那么在一些多项式题目中 what should we do?
于是一个新东西出来了。数论变换是离散傅里叶变换（DFT）在数论基础上的实现。快速数论变换是快速傅里叶变换在数论基础上的实现。
NTT 解决的是多项式乘法带模数的情况。
接下来介绍前置知识：

• 离散傅里叶变换
  见上文
• 生成子群
• 原根
• 离散对数

这一段其实挺离谱的，因为我真的找不到什么特别大的关系。
我们定义 $P$ 为素数，$g$ 为 $P$ 的原根。
原根有以下与单位根相同的性质：

1.$\mathrm{\forall }i,j\in [0,n-1],i\ne j,{\omega }_n^i\ne {\omega }_n^j$

2.${\omega }_{2n}^{2k}={\omega }_n^k$ 等价于 $({\omega }_{2n}{)}^2={\omega }_n$

3.${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$

4.$\sum _0^{n-1}({\omega }_n^k{)}^i=0$

于是我们有一个结论：${\omega }_n\equiv g^{\frac{p-1}{n}}(\mathrm{mod}p)$。
于是 FFT 中的 ${\omega }_n$ 都可以替换掉。
常用 $p=998244353,g=3,invg=332748118$。
接下来的部分同 FFT 啦。