A 签到题
B FFT/NTT快速傅立叶
P3803 【模板】多项式乘法(FFT) & P1919 【模板】高精度乘法 | A*B Problem 升级版
参考资料:
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)\) 就行了)
标签:limits,卷积,Contest5385,FFT,2i,omega,sum From: https://www.cnblogs.com/AugustLight/p/18345691