Contest5385 - FFT

Contest

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 礼物

P3723 [AH2017/HNOI2017] 礼物

注意到在一个手环上 \(+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)\) 就行了）

posted @ 2024-08-06 17:38 August_Light 阅读(88) 评论(4) 收藏举报

刷新页面返回顶部