【算法学习笔记】快速傅里叶变换

快速傅里叶变换

快速傅里叶变换（Fast Fourier Transform, FTT）在ACM/OI中最主要的应用是计算多项式乘法。

多项式的系数表示和点值表示

假设\(f(x)\)为\(x\)的\(n\)阶多项式，则其可以表示为：

\[f(x)=\sum_{i=0}^na_ix^i \]

这里的\(n+1\)个系数\(\{a_0,a_1,\cdots,a_n\}\)就称为多项式\(f(x)\)的系数表示。

另一方面，我们也可以把\(f(x)\)看成是一个关于\(x\)的函数，我们可以取\(n+1\)个不同的\(x_i\)，用\(\{(x_0,f(x_0)),(x_1,f(x_1)),\cdots(x_n,f(x_n))\}\)这\(n+1\)个数值对来唯一确定\(f(x)\)，这种表示形式就称为多项式\(f(x)\)的点值表示。

点值表示与多项式乘法的关系

假设我们现在要求的是\(F(x)=f(x)\cdot g(x)\)，如果我们已知\(f(x)\)和\(g(x)\)的点值表示，那么我们可以非常容易地得到\(F(x)\)的点值表示为

\[\{(x_0,f(x_0)g(x_0)),(x_1,f(x_1)g(x_1)),\cdots,(x_n,f(x_n)g(x_n))\} \]

注意这里的\(n\)实际上要取到\(f(x)\)和\(g(x)\)的阶数之和。

现在的关键问题是，如何快速将这一点值表示转换为系数表示。

FFT的实现

为了解决这一问题，我们首先考虑其逆问题，也即：如何从系数表示快速计算点值表示。

FFT

暴力计算\(n\)对点值的总时间复杂度为\(O(n^2)\)。如何优化呢？我们希望我们选择的\(n\)个\(x_i\)之间存在一定的关系，使得我们可以复用\(x_i^k\)的计算结果。那么，应该如何选择呢？

前人的经验告诉我们，可以选择单位复根\(\omega_n^i\)。它有三个重要的性质：

\[\omega_n^n=1 \]

\[\omega_n^i=\omega_{2n}^{2i} \]

\[\omega_{2n}^{n+i}=-\omega_{2n}^i \]

利用上述这三个性质，我们可以实现计算过程的简化。

不妨考虑一个最高阶为7阶的多项式

\[f(x)=a_0+a_1x^1+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7 \]

可以把奇偶项分别处理

\[\begin{aligned} f(x) &=(a_0+a_2x^2+a_4x^4+a_6x^6)+x(a_1+a_3x^2+a_5x^4+a_7x^6) \\ &=G(x^2)+xH(x^2) \end{aligned} \]

从而

\[\text{DFT}(f(x))=\text{DFT}(G(x^2))+x\text{DFT}(H(x^2)) \]

这时把单位复根\(\omega_n^k\)（\(k<n/2\)）代入，可以得到

\[\begin{aligned} \text{DFT}(f(\omega_n^k))&=\text{DFT}(G(\omega_n^{2k}))+\omega_n^k\text{DFT}(H(\omega_n^{2k})) \\ &=\text{DFT}(G(\omega_{n/2}^k))+\omega_n^k\text{DFT}(H(\omega_{n/2}^k)) \end{aligned} \]

而另一方面，代入\(\omega_n^{k+n/2}\)可以得到

\[\begin{aligned} \text{DFT}(f(\omega_n^{k+n/2}))&=\text{DFT}(G(\omega_n^{2k+n}))+\omega_n^k\text{DFT}(H(\omega_n^{2k+n})) \\ &=\text{DFT}(G(\omega_{n/2}^{k+n/2}))+\omega_n^{k+n/2}\text{DFT}(H(\omega_{n/2}^{k+n/2})) \\ &=\text{DFT}(G(\omega_{n/2}^k))-\omega_n^k\text{DFT}(H(\omega_{n/2}^k)) \end{aligned} \]

因此，我们只要求得\(\text{DFT}(G(\omega_{n/2}^k))\)和\(\text{DFT}(H(\omega_{n/2}^k))\)，就可以同时求得\(\text{DFT}(f(\omega_n^k))\)和\(\text{DFT}(f(\omega_n^{k+n/2}))\)，这样就把问题规模缩小了一半。

使用同样的方法对\(\text{DFT}(G(\omega_{n/2}^k))\)和\(\text{DFT}(H(\omega_{n/2}^k))\)进行递归求解，我们有

\[T(n)=2T(n/2) \]

可知总的时间复杂度为\(O(n\log n)\)。

在这一过程中，我们默认\(n/2\)总是整数，因此我们需要\(n=2^k\)。所以在计算之前，我们要先对系数补0，使得总的项数变为2的幂次。

逆FFT

将FFT的运算过程看做一个矩阵乘法，逆FFT，也即从点值表示求取系数表示的过程，可以视为左乘逆矩阵。在点值表示的点选取为\(\omega_n^k\)时，FFT矩阵\(\mathbb{A}(\omega_n^k)\)的逆矩阵恰好为\(\frac{1}{n}\mathbb{A}(\omega_n^{-k})\)，因此可以复用FFT的计算过程，只需要加上一个标志变量来表示当前是在进行FFT还是IFFT。

模板题：洛谷 P3803 - 多项式乘法（FFT）

下面给出了本题的递归实现。

Code(C++)

#include <cmath>
#include <complex>
#include <iostream>
#define MAXN (1 << 22)
using namespace std;
typedef complex cd;
const cd I{0, 1};
cd tmp[MAXN], a[MAXN], b[MAXN];
void fft(cd *f, int n, int rev) {
    if (n == 1) return;
    for (int i = 0; i < n; ++i) tmp[i] = f[i];
    for (int i = 0; i < n; ++i) {
        if (i & 1) f[n / 2 + i / 2] = tmp[i];
        else
            f[i / 2] = tmp[i];
    }
    cd *g = f, *h = f + n / 2;
    fft(g, n / 2, rev), fft(h, n / 2, rev);
    cd omega = exp(I * (2 * M_PI / n * rev)), now = 1;
    for (int k = 0; k < n / 2; ++k) {
        tmp[k] = g[k] + now * h[k];
        tmp[k + n / 2] = g[k] - now * h[k];
        now *= omega;
    }
    for (int i = 0; i < n; ++i) f[i] = tmp[i];
}
int main() {
    int n, m;
    cin >> n >> m;
    int k = 1 << (32 - __builtin_clz(n + m + 1));
    for (int i = 0; i <= n; ++i) cin >> a[i];
    for (int j = 0; j <= m; ++j) cin >> b[j];
    fft(a, k, 1);
    fft(b, k, 1);
    for (int i = 0; i < k; ++i) a[i] *= b[i];
    fft(a, k, -1);
    for (int i = 0; i < k; ++i) a[i] /= k;
    for (int i = 0; i < n + m + 1; ++i) cout << (int)round(a[i].real()) << " ";
}

上述递归方法的常数较大，不能通过洛谷P3803的最后两个测试点。

为了改写非递归方法，我们引入蝴蝶变换的概念。

蝴蝶变换

继续使用前面的例子，经过第一步分治，将原来的系数分为两组：

\[\{a_0,a_2,a_4,a_6\},\{a_1,a_3,a_5,a_7\} \]

继续进行第二步分治，得到四组系数：

\[\{a_0,a_4\},\{a_2,a_6\},\{a_1,a_5\},\{a_3,a_7\} \]

最后一步分治，得到八组系数：

\[\{a_0\},\{a_4\},\{a_2\},\{a_6\},\{a_1\},\{a_5\},\{a_3\},\{a_7\} \]

所谓蝴蝶变换，指的就是从\({a_0,a_1,\cdots,a_{n-1}}\)这一原始系数序列，变换得到最后一步分治后的系数序列。

观察后可以发现，在蝴蝶变换的最终结果中，系数下标的二进制表示恰好是其所在位置二进制表示的逆序，因此，可以利用这一规律来求取蝴蝶变换的结果。

直接利用规律来计算的复杂度是\(O(n\log n)\)，如果从小到大递推实现，复杂度则为\(O(n)\)。

FFT的非递归实现

下面给出了洛谷P3803的非递归实现。

Code(C++)

#include <cmath>
#include <complex>
#include <iostream>
#define MAXN (1 << 22)
using namespace std;
typedef complex cd;
const cd I{0, 1};
cd a[MAXN], b[MAXN];
void change(cd *f, int n) {
    int i, j, k;
    for (int i = 1, j = n / 2; i < n - 1; i++) {
        if (i < j) swap(f[i], f[j]);
        k = n / 2;
        while (j >= k) {
            j = j - k;
            k = k / 2;
        }
        if (j < k) j += k;
    }
}
void fft(cd *f, int n, int rev) {
    change(f, n);
    for (int len = 2; len <= n; len <<= 1) {
        cd omega = exp(I * (2 * M_PI / len * rev));
        for (int j = 0; j < n; j += len) {
            cd now = 1;
            for (int k = j; k < j + len / 2; ++k) {
                cd g = f[k], h = now * f[k + len / 2];
                f[k] = g + h, f[k + len / 2] = g - h;
                now *= omega;
            }
        }
    }
    if (rev == -1)
        for (int i = 0; i < n; ++i) f[i] /= n;
}
int main() {
    int n, m;
    cin >> n >> m;
    int k = 1 << (32 - __builtin_clz(n + m + 1));
    for (int i = 0; i <= n; ++i) cin >> a[i];
    for (int j = 0; j <= m; ++j) cin >> b[j];
    fft(a, k, 1);
    fft(b, k, 1);
    for (int i = 0; i < k; ++i) a[i] *= b[i];
    fft(a, k, -1);
    for (int i = 0; i < n + m + 1; ++i) cout << (int)round(a[i].real()) << " ";
}

学习资源

Matters Computational

第二十一章快速傅里叶变换

练习题

裸FFT并不可怕，本身FFT的码量并不算大，背一背也不是多大的事，关键是如何看出一道题目是FFT。

SPOJ - ADAMATCH

如果暴力枚举子串，时间复杂度为\(O(|r|^2)\)，显然不行。如何降低复杂度呢？

提示一

首先考虑字母'A'。不妨把字符串为'A'的位置设为\(1\)，其余位置设为\(0\)。看起来似乎可以进行多项式乘法，但乘法的结果似乎没有明显的意义。

提示二

如果把r串逆序呢？看看此时乘积的每一项有怎样的含义。

参考代码（C++）

#include <cmath>
#include <complex>
#include <cstring>
#include <iostream>
#include <vector>
#define MAXN (1 << 22)

using namespace std;
typedef complex<double> cd;
const cd I{0, 1};
cd a[MAXN], b[MAXN];
void change(cd *f, int n) {
    for (int i = 1, j = n / 2; i < n - 1; i++) {
        if (i < j) swap(f[i], f[j]);
        int k = n / 2;
        while (j >= k) {
            j = j - k;
            k = k / 2;
        }
        if (j < k) j += k;
    }
}
void fft(cd *f, int n, int rev) {
    change(f, n);
    for (int len = 2; len <= n; len <<= 1) {
        cd omega = exp(I * (2 * M_PI / len * rev));
        for (int j = 0; j < n; j += len) {
            cd now = 1;
            for (int k = j; k < j + len / 2; ++k) {
                cd g = f[k], h = now * f[k + len / 2];
                f[k] = g + h, f[k + len / 2] = g - h;
                now *= omega;
            }
        }
    }
    if (rev == -1)
        for (int i = 0; i < n; ++i) f[i] /= n;
}
int main() {
    string s, r;
    cin >> s >> r;
    int n = s.size(), m = r.size();
    int k = 1 << (32 - __builtin_clz(n + m + 1));
    vector<int> cnt(k);
    for (char c : "ACGT") {
        memset(a, 0, sizeof(a));
        memset(b, 0, sizeof(b));
        for (int i = 0; i < n; ++i) a[i] = s[i] == c;
        for (int i = 0; i < m; ++i) b[i] = r[m - i - 1] == c;
        fft(a, k, 1);
        fft(b, k, 1);
        for (int i = 0; i < k; ++i) a[i] *= b[i];
        fft(a, k, -1);
        for (int i = 0; i < k; ++i) cnt[i] += (int)round(a[i].real());
    }
    int ans = m;
    for (int i = m - 1; i < n; ++i) ans = min(ans, m - cnt[i]);
    cout << ans;
}

SPOJ - TSUM

如果暴力枚举，时间复杂度为\(O(n^3)\)，显然不行。如何降低复杂度呢？

提示一

加法可以变为多项式的乘法。

提示二

如何去除包含重复元素的项？

参考代码（C++）

#include <cmath>
#include <complex>
#include <iostream>
#include <vector>

#define MAXN 131072
#define OFFSET 20000

using namespace std;
typedef complex<double> cd;
const cd I{0, 1};

void change(vector<cd> &f, int n) {
    for (int i = 1, j = n / 2; i < n - 1; i++) {
        if (i < j) swap(f[i], f[j]);
        int k = n / 2;
        while (j >= k) {
            j = j - k;
            k = k / 2;
        }
        if (j < k) j += k;
    }
}

void fft(vector<cd> &f, int n, int rev) {
    change(f, n);
    for (int len = 2; len <= n; len <<= 1) {
        cd omega = exp(I * (2 * M_PI / len * rev));
        for (int j = 0; j < n; j += len) {
            cd now = 1;
            for (int k = j; k < j + len / 2; ++k) {
                cd g = f[k], h = now * f[k + len / 2];
                f[k] = g + h, f[k + len / 2] = g - h;
                now *= omega;
            }
        }
    }
    if (rev == -1)
        for (int i = 0; i < n; ++i) f[i] /= n;
}

int main() {
    int n;
    cin >> n;
    vector<cd> a(MAXN), a2(MAXN);
    vector<int> a3(MAXN);
    for (int i = 0; i < n; ++i) {
        int m;
        cin >> m;
        a[m + OFFSET] = cd{1, 0};
        a2[(m + OFFSET) << 1] = cd{1, 0};
        a3[(m + OFFSET) * 3] = 1;
    }
    vector<cd> tot(a), b(a);
    fft(tot, MAXN, 1);
    fft(b, MAXN, 1);
    fft(a2, MAXN, 1);
    for (int i = 0; i < MAXN; ++i) tot[i] *= b[i] * b[i], a2[i] *= b[i];
    fft(tot, MAXN, -1);
    fft(a2, MAXN, -1);
    for (int i = 0; i < MAXN; ++i) {
        int cnt1 = round(tot[i].real());  // ABC, with permutation
        int cnt2 = round(a2[i].real());   // AAB, no permutation
        int cnt3 = a3[i];                 // AAA
        int cnt = (cnt1 - cnt2 * 3 + cnt3 * 2) / 6;
        if (cnt > 0) cout << i - OFFSET * 3 << " : " << cnt << endl;
    }
}

posted @ 2021-03-15 21:22 Koshkaaa 阅读(556) 评论(1) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

Koshkaaa の Site

Bug writer. & Cpp-Developer、Rustacean. Compilers.

【算法学习笔记】快速傅里叶变换

快速傅里叶变换

多项式的系数表示和点值表示

点值表示与多项式乘法的关系

FFT的实现

FFT

逆FFT

模板题：洛谷 P3803 - 多项式乘法（FFT）

蝴蝶变换

FFT的非递归实现

学习资源

Matters Computational

练习题

SPOJ - ADAMATCH

提示一

提示二

参考代码（C++）

SPOJ - TSUM

提示一

提示二

参考代码（C++）