Contest5385 - FFT
Published:
A 签到题
B FFT/NTT快速傅立叶
| P3803 【模板】多项式乘法（FFT） & [P1919 【模板】高精度乘法 | A*B Problem 升级版](https://www.luogu.com.cn/problem/P1919) |
参考资料：
FFT 总体思路
FFT 处理循环卷积问题，而卷积问题通用的处理方法为：
- $A = FFT(a)$
- $B = FFT(b)$
- $C = A \cdot B$
- $c = IFFT(C)$
其中 $IFFT$ 是 $FFT$ 的逆变换，即 $a = IFFT(FFT(a))$。
在 OI 中遇到卷积问题一般都这样四步走解决，包括但不限于 $\gcd$ 卷积，$\text{bitand}$ 卷积等。
小细节：循环卷积
FFT 处理的是循环卷积问题，形如求数列 $c$ 使得：
\[c_r = \sum\limits_{p,q} [(p + q) \bmod n = r] a_p b_q\]而我们平时的多项式乘法，指的是求数列 $c$ 使得：
\[c_r = \sum\limits_{p,q} [p + q = r] a_p b_q\]但是只要 $n$ 开的足够大（比两个多项式度数之和更大），两个东西就是一样的。
推公式：FFT & IFFT 与反演
卷积问题往往和反演有关。
令 $\omega_n^k$ 为 $n$ 次单位根，我们有单位根反演：
\[\frac 1 n \sum\limits_{i=0}^{n-1} (\omega_n^k)^i = [k \bmod n = 0]\]套一下：
\[\begin{aligned} & [(p + q) \bmod n = r] \\ =& [(p + q - r) \bmod n = 0] \\ =& \frac 1 n \sum\limits_{i=0}^{n-1} (\omega_n^{p+q-r})^i \\ =& \frac 1 n \sum\limits_{i=0}^{n-1} (\omega_n^p)^i (\omega_n^q)^i (\omega_n^{-r})^i \\ \end{aligned}\]套一下：
\[\begin{aligned} c_r =& \sum\limits_{p} \sum\limits_{q} [(p + q) \bmod n = r] a_p b_q \\ =& \sum\limits_{p} \sum\limits_{q} \frac 1 n \sum\limits_{i=0}^{n-1} (\omega_n^p)^i (\omega_n^q)^i (\omega_n^{-r})^i a_p b_q \\ =& \frac 1 n \sum\limits_{i=0}^{n-1} (\omega_n^{-r})^i \left( \sum\limits_{p} (\omega_n^p)^i a_p \right) \left( \sum\limits_{q} (\omega_n^q)^i b_q \right) \end{aligned}\]我们就得到了 FFT 的核心公式：
FFT:
\[f_k = \sum\limits_{i=0}^{n-1} \omega_n^{ki} g_i\]IFFT:
\[g_k = \frac 1 n \sum\limits_{i=0}^{n-1} \omega_n^{-ki} f_i\]从这个公式可以看出，IFFT 可以方便地套用 FFT 的代码：
void IFFT(Complex a[], int n) {
FFT(a, n);
reverse(a+1, a+n);
for (int i = 0; i < n; i++)
a[i] = a[i].real() / n; // 这里只保留了实部
}
接下来我们就不用为 IFFT 操心了，专心于 FFT 的实现即可。
FFT 的实现
\[f_k = \sum\limits_{i=0}^{n-1} \omega_n^{ki} g_i\]我们要在 $O(n \log n)$ 完成这个式子的全部计算。不难想到分治。
递归版
考虑按下标 $i$ 的奇偶性分类，令 $m = \frac n 2$：
\[\begin{aligned} f_k =& \sum\limits_{i=0}^{m - 1} a_{2i} (\omega_n^k)^{2i} + \sum\limits_{i=0}^{m-1} a_{2i+1} (\omega_n^k)^{2i+1} \\ =& \sum\limits_{i=0}^{m - 1} a_{2i} (\omega_n^k)^{2i} + \omega_n^k \sum\limits_{i=0}^{m-1} a_{2i+1} (\omega_n^k)^{2i} \\ =& \sum\limits_{i=0}^{m - 1} a_{2i} (\omega_m^k)^{i} + \omega_n^k \sum\limits_{i=0}^{m-1} a_{2i+1} (\omega_m^k)^{i} \\ \end{aligned}\]显然可以递归计算。
void FFT(Complex a[], int n) {
if (n == 1) return;
int m = n / 2;
Complex a0[m], a1[m];
for (int i = 0; i < m; i++) {
a0[i] = a[i << 1];
a1[i] = a[i << 1 | 1];
}
FFT(a0, m), FFT(a1, m);
auto W = Complex(cos(PI / m), sin(PI / m)); // n 次单位根
auto w = Complex(1, 0); for (int i = 0; i < m; i++) {
a[i] = a0[i] + w * a1[i];
a[i + m] = a0[i] - w * a1[i];
w = w * W;
}
}
注意：分治时需要保证 $n$ 是 $2$ 的幂，所以应该在 FFT 前先将 $n$ 补全到比两个多项式度数之和更大的 $2$ 的幂。
迭代版
递归 FFT 的时间复杂度虽然是正确的，但是因为递归和动态开空间，所以常数十分巨大。如果有一道题要做十次甚至九次 FFT，可能会导致超时。改为迭代可以获得更高效率。
观察可得，递归到最底层的时候，$a_i$ 的值为原先的 $a_{rev(i)}$，其中 $rev(i)$ 为 $i$ 的二进制表示反转后的结果，可以这样 $O(n)$ 递推：
rep(i, 0, tot-1)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (__lg(tot) - 1));
那么我们可以写出以下代码：
void FFT(Complex a[], int n) {
for (int i = 0; i < n; ++i)
if (rev[i] > i) // 只能换一次，换两次等于白换
swap(a[i], a[rev[i]]);
for (int len = 2, m = 1; len <= n; m <<= 1, len <<= 1) { // m 始终是 len 的一半
auto W = Complex(cos(PI / m), sin(PI / m)); // len 次单位根
for (int l = 0, r = len-1; r < n; l += len, r += len) {
auto w = Complex(1, 0);
for (int p = l; p < l + m; p++) {
Complex x = a[p] + w * a[p + m], y = a[p] - w * a[p + m];
a[p] = x, a[p + m] = y;
w = w * W;
}
}
}
}
至此，我们可以在 $\Theta(n \log n)$ 的时间复杂度下，以不大的常数，完成一次 $\Theta(n)$ 的多项式乘法。
C 礼物
注意到在一个手环上 $+c$ 等价于在另一个手环上 $-c$，所以就不用细分是加在哪条上。
$c$ 范围很小可以枚举，把 $+c$ 带进去展开之后，发现决定答案的是一个 $\sum\limits_{i=0}^{n-1} x_{(i+k) \bmod n} y_i$ 的结构。
断环为链，翻转 $x$ 之后就变成了一个卷积，直接 FFT 就行了。
时间复杂度 $O(n \log n + m)$。（精细一点的话 $+m$ 都不需要，直接用数学算出 $c$，$O(n \log n)$ 就行了）