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

多项式部分是 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\limits_{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\limits_{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\%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_+}\)。
设多项式 \(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) + xF_2(x^2)\)。
\(\forall k \in [0,\frac{n}{2})\)，将 \(\omega_{n}^{k}\) 代入得：

\[\begin{split} 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{split} \]

分治惯用套路，将 \(\omega_{n}^{k+\frac{n}{2}}\) 代入得：

\[\begin{split} 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{split} \]

我们欣喜地发现，\(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\) 的值用高中数学的知识就可求出 \(\big(\cos(\frac{2\pi}{n}),\sin(\frac{2\pi}{n})\big)\)。

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\limits_{i=0}^{n-1}(\omega_n^k)^i\cdot f_i\)
先给出结论： \(n\cdot f_k = \sum\limits_{i=0}^{n-1}(\omega_n^{-k})^i \cdot g_i\)
这个式子可以将点值转换为系数。

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

1. \(j=k\) 贡献 \(\sum\limits_{i=0}^{n-1}\omega_n^0\cdot f_k = n \cdot f_k\)
2. \(j \not = k\)，假设 \(p = j - k\)，贡献 \(\sum\limits_{i=0}^{n-1}\omega_n^{ip}f_{k + p}=\omega_{n}^{p}(\sum\limits_{i=0}^{n-1}\omega_n^i)f_{k+p}\)
   由等比数列求和公式：\(\sum\limits_{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} = \big(\cos (\frac{2\pi}{n}), -\sin (\frac{2\pi}{n})\big)\)。
这一段和 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.\(\forall i,j \in [0, n - 1], i \not = j, \omega_n^i \not = \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\limits_{0}^{n-1}(\omega_{n}^{k})^i=0\)

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