前前后后鸽了好久，终于抄完了。
单位根
概念
在复数下，满足 \(x^n = 1\) 的 \(x\) 称为 \(n\) 次单位根。
\(n\) 次单位根一共有 \(n\) 个。
将所有的单位根按照辐角大小排列，第 \(k\) 个（\(0 \leq k < n\)）个 \(n\) 次单位根为：
\(x_k = e^{i \frac{2 k \pi}{n}}\)
所有的单位根模都是 \(1\)，\(n\) 个单位根平分单位圆。
本原单位根：\(0\) 到 \(n - 1\) 次方的值能生成所有 \(n\) 次单位根的单位根称为 \(n\) 次本原单位根。
欧拉公式：\(e^{i \pi} = 1\)
\(x_1 = e^{i \frac{2 \pi}{n}}\) 是一个本原单位根。
记 \(n\) 次本原单位根为 \(\omega_n = e^{i \frac{2 \pi}{n}} = \cos \frac{2 \pi}{n} + i \sin \frac{2 \pi}{n}\)
性质
设 \(n\) 是一个正偶数，则：
\((\omega_n^k)^2 = w_{\frac{n}{2}}^k\)
\(\omega_{n}^{\frac{n}{2} + k} = -\omega_n^k\)
算法
求出一个 \(n\) 次多项式在每个 \(n\) 次单位根下的点值称为 离散傅里叶变换（DFT），而重新从这 \(n\) 个点值得到 \(n\) 次多项式的过程称为 离散傅里叶逆变换（IDFT）
对于要进行 FFT 的多项式，先将它的次数补到 \(2\) 的整次幂，记其次数为 \(n - 1\)，令 \(m = \frac{n}{2}\)
DFT
考虑求一个长度为 \(n\) 的数列 \(b\)，其中：
\(\sum\limits_{i = 0}^{n - 1} b_i = \sum\limits_{i = 0}^{n - 1} a_i \times \omega_n^{ik}\)
即这个数列的第 \(k\) 项是原多项式在 \(n\) 次单位根的 \(k\) 次幂处的取值。
FFT
FFT 是对 \(O(n^2)\) 的朴素 DFT 的优化。
考虑把上面的和式按照下标分类，即：
\(A(x) = \sum\limits_{i = 0}^{n - 1} a_i \cdot x_i = \sum\limits_{i = 0}^{m - 1} a_{2i} \cdot x^{2i} + \sum\limits_{i = 0}^{m - 1} a_{2i + 1} \cdot x^{2i + 1}\)
后半部分提出一个 \(x\)，得：
\(A(x) = \sum\limits_{i = 0}^{n - 1} a_i \cdot x_i = \sum\limits_{i = 0}^{m - 1} a_{2i} \cdot x^{2i} + x \sum\limits_{i = 0}^{m - 1} a_{2i + 1} \cdot x^{2i} = \sum\limits_{i = 0}^{m - 1} a_{2i} \cdot (x^2)^i + x \sum\limits_{i = 0}^{m - 1} a_{2i + 1} \cdot (x^2)^i\)
设 \(A_0(x), A_1(x)\) 是两个 \(m - 1\) 次多项式，使得 \(A_0(x) = \sum\limits_{i = 0}^{m - 1} a_{2i} x^i, A_1(x) = \sum\limits_{i = 0}^{m - 1} a_{2i + 1} x^i\)
那么 \(A(x) = A_0(x^2) + x A_1(x^2)\)
这说明只要求出 \(A_0, A_1\) 在各处的点值，我们就可以在 \(O(n)\) 时间内计算 \(A\) 在各处的点值！并且子问题的结构是类似地，完全可以递归处理！
但是这样递归下去的复杂度是 \(O(n^2)\)
所以我们可以考虑一下单位根的性质：
\((\omega_n^k)^2 = \omega_m^k\)
\(\omega_{n}^{m + k} = -\omega_n^k\)
考虑 \(m\) 次以内的点值和高于 \(m\) 次的点值之间的关系：
\(\forall 0 \leq k < m, A(\omega_n^k) = A_0((\omega_n^k)^2) + \omega_n^k A_1((\omega_n^k)^2)\)
根据第一个式子化简得到：
\(A(\omega_n^k) = A_0(\omega_m^k) + \omega_n^k A_1(\omega_m^k)\)
对于大于等于 \(m\) 次的点值：
\(A(\omega_n^{m + k}) = A_0((\omega_n^{m + k})^2) + \omega_n^{m + k} A_1((\omega_n^{m + k})^2)\)
根据第二个式子化简得到：
\(A(\omega_n^{m + k}) = A_0((\omega_n^k)^2) - \omega_n^k A_1((\omega_n^k)^2)\)
再根据第一个式子得到：
\(A(\omega_n^{m + k}) = A_0(\omega_m^k) - \omega_n^k A_1(\omega_m^k)\)
比对一下：\(A(\omega_n^k) = A_0(\omega_m^k) + \omega_n^k A_1(\omega_m^k)\)
差别只有后半部分的系数！！！1
所以我们只需要把小于 \(m\) 次的部分递归处理出来，然后再处理后半部分即可。
时间复杂度优化到 \(O(n \log n)\)
上面的推导称为 蝴蝶操作。
IDFT
不太会证明，暂且先放个结论，等以后再回来填坑。
假设当前已知一个 \(n - 1\) 次多项式 \(A(x) = \sum\limits_{i = 0}^{n - 1} a_i \cdot x^i\) 经过 DFT 后得到的点值序列 \(b\)，则：
\(b_k = \sum\limits_{i = 0}^{n - 1} a_i \cdot \omega_n^{ik}\)
结论是：\(a_k = \frac{1}{n} \sum\limits_{i = 0}^{n - 1} b_i \omega_n^{-ki}\)
这里有一个消去引理的知识，之后再补回来。
IFFT
如何快速求 \(B\) 在 \(\omega_n^{-ki}\) 处的取值？
根据几何意义易知：\(\omega_n^{-k} = \omega_n^{n - k}\)
所以只需要正常求 \(B\) 在 \(\omega_n^k\) 处的取值，然后对数组取反即可，magic！
优化
递归改成迭代会快很多。
考虑画一下递归调用的原数组下标：
stage 1: 0 1 2 3 4 5 6 7
stage 2: 0 2 4 6|1 3 5 7
stage 3: 0 4|2 6|1 5|3 7
stage 4: 0|4|2|6|1|5|3|7
人类智慧：stage 4 二进制写法反过来刚好是：
000, 001, 010, 011, 100, 110, 101, 111
也就是 01234567，magic！
所以只需要把原本的 \(a\) 数组按照二进制下标反转后的大小排序，倒数第二层的蝴蝶操作就是合并 \(a\) 中相邻的项。
用类似数位 dp 的方法可以快速求出每个值二进制反转后的值：
rev[i] = (rev[i >> 1] >> 1 | (i & 1 ? k >> 1 : 0));
因为 \(rev(rev(i)) = i\)，所以对于一个下标 \(p\)，若 \(rev(p) = q\)，那么 \(rev(q) = p\)。所以只需要对于 \(rev(i) > i\) 的位置交换一下就能排序力，这个称为位逆序置换。
代码
#include <cstdio>
#include <cmath>
#include <complex>
#include <iostream>
#include <algorithm>
using namespace std;
const int sz = 5e6 + 5;
const double PI = acos(-1);
int n, m;
int rev[sz];
complex<double> F[sz], G[sz];
void calc_rev(int k) { for (int i = 0; i < k; i++) rev[i] = (rev[i >> 1] >> 1 | (i & 1 ? k >> 1 : 0)); }
void FFT(complex<double> *A, int n)
{
for (int i = 1; i < n; i++)
if (rev[i] > i) swap(A[i], A[rev[i]]);
for (int len = 2, m = 1; len <= n; m = len, len <<= 1)
{
complex<double> W(cos(PI / m), sin(PI / m)), w(1.0, 0.0);
for (int l = 0, r = len - 1; r <= n; l += len, r += len)
{
auto w0 = w;
for (int p = l; p < l + m; p++)
{
auto x = A[p] + w0 * A[p + m], y = A[p] - w0 * A[p + m];
A[p] = x, A[p + m] = y;
w0 *= W;
}
}
}
}
void IFFT(complex<double> *A, int n)
{
FFT(A, n);
reverse(A + 1, A + n);
}
int main()
{
scanf("%d%d", &n, &m);
for (int i = 0, v; i <= n; i++) scanf("%d", &v), F[i] = v;
for (int i = 0, v; i <= m; i++) scanf("%d", &v), G[i] = v;
int len = n + m, k = 1;
while (k <= len) k <<= 1;
calc_rev(k);
FFT(F, k), FFT(G, k);
for (int i = 0; i < k; i++) F[i] *= G[i];
IFFT(F, k);
for (int i = 0; i <= len; i++) printf("%d ", (int)(F[i].real() / k + 0.5));
return 0;
}