快速傅立叶变换(FFT)
快速傅立叶变换
根据oiwiki整理
函数的两种表示方法——点值表示法和系数表示法
系数表示法就是常见的函数表示形式: \(f(x) = a_0 + a_1x^1 + a_2x^2 + \cdots + a_nx^n\)
点值表达式就是从函数 \(f(x)\)中取出一些点 \((x,f(x))\)， 并且可以求出对应系数表达式的相应系数，对于 \(n\)个系数的函数，最少取出 \(n+1\)个点即可通过高斯消元求出相应系数。
通俗地说，多项式由系数表示法转为点值表示法的过程，就是 \(DFT\)(离散傅立叶变换) 的过程。相对地，把一个多项式的点值表示法转化为系数表示法的过程，就是 \(IDFT\)。而 \(FFT\) 就是通过取某些特殊的 \(x\) 的点值来加速 \(DFT\) 和 \(IDFT\) 的过程。
单位复根
对于 \(f(x)\)的表达式，当 \(x = 1\)和\(x = -1\)时 \(x^i\) 很好计算，但是我们需要 \((n+1)\)个点才可以求出 \(f(x)\)的表达式，所以这里引入了单位复根的概念。
我们称. \(x^n = 1\)在复数意义下的解是 \(n\)次复根，并且有 \(n\)个解，设 \(\omega_n = e^{\frac {2\pi i} n}\)，则 \(x^n = 1\)的解集为：\(\{\omega_n^k|k = 0, 1, 2, \cdots, (n-1)\}\)，其中 \(\omega_n\)为 \(n\)次复根，同时由欧拉公式可得: \(\omega_n = e^{\frac {2\pi i} n} = cos(\frac {2\pi}n) + i * sin(\frac {2\pi}n)\)
\(n\)次复根就是复平面把单位圆 \(n\)等分的第一个角所对应的向量。
例如当 \(n = 4\)时， \(\omega_4 = i\)，所以 \(i\)就是 \(4\)次复根。
单位复根的性质：
- \(\omega_n^n = 1\)
- \(\omega_n^k = \omega_{2n}^{2k}\)
- \(\omega_{2n}^{k+n} = -\omega_{2n}^k\)
\(f(x)\) 的求解
假如 \(f(x)\)的系数只有7项且如下所示：
将奇偶幂次的项分开：
令: \(G(x) = a_0 + a_2x^1 + a_4x^2 + a_6x^2, H(x) = a_1 + a_3x^1 + a_5x^2 + a_7x^3\)
则
根据 \(DFT\)的性质计算：
因此可以递归求解 \(f(x)\)
Code(递归版)
#include <cmath>
#include <complex>
typedef std::complex<double> Comp; // STL complex
const Comp I(0, 1); // i
const int MAX_N = 1 << 20;
Comp tmp[MAX_N];
void DFT(Comp *f, int n, int rev) { // rev=1,DFT; rev=-1,IDFT
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];
}
Comp *g = f, *h = f + n / 2;
DFT(g, n / 2, rev), DFT(h, n / 2, rev); // 递归 DFT
Comp cur(1, 0), step(cos(2 * M_PI / n), sin(2 * M_PI * rev / n));
// Comp step=exp(I*(2*M_PI/n*rev)); // 两个 step 定义是等价的
for (int k = 0; k < n / 2; ++k) {
tmp[k] = g[k] + cur * h[k];
tmp[k + n / 2] = g[k] - cur * h[k];
cur *= step;
}
for (int i = 0; i < n; ++i) f[i] = tmp[i];
}
位逆序置换优化
这个算法还可以从“分治”的角度继续优化。我们每一次都会把整个多项式的奇数次项和偶数次项系数分开，一直分到只剩下一个系数。但是，这个递归的过程需要更多的内存。因此，我们可以先“模仿递归”把这些系数在原数组中“拆分”，然后再“倍增”地去合并这些算出来的值。
Code(O(nlogn))
/*
* 进行 FFT 和 IFFT 前的反置变换
* 位置 i 和 i 的二进制反转后的位置互换
*len 必须为 2 的幂
*/
void change(Complex y[], int len) {
int i, j, k;
for (int i = 1, j = len / 2; i < len - 1; i++) {
if (i < j) swap(y[i], y[j]);
// 交换互为小标反转的元素，i<j 保证交换一次
// i 做正常的 + 1，j 做反转类型的 + 1，始终保持 i 和 j 是反转的
k = len / 2;
while (j >= k) {
j = j - k;
k = k / 2;
}
if (j < k) j += k;
}
}
Code(O(n))
// 同样需要保证 len 是 2 的幂
// 记 rev[i] 为 i 翻转后的值
void change(Complex y[], int len) {
for (int i = 0; i < len; ++i) {
rev[i] = rev[i >> 1] >> 1;
if (i & 1) { // 如果最后一位是 1，则翻转成 len/2
rev[i] |= len >> 1;
}
}
for (int i = 0; i < len; ++i) {
if (i < rev[i]) { // 保证每对数只翻转一次
swap(y[i], y[rev[i]]);
}
}
return;
}
FFT板子:
/*
* 做 FFT
* len 必须是 2^k 形式
* on == 1 时是 DFT，on == -1 时是 IDFT
*/
void fft(Complex y[], int len, int on) {
change(y, len);
for (int h = 2; h <= len; h <<= 1) { // 模拟合并过程
Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h)); // 计算当前单位复根
for (int j = 0; j < len; j += h) {
Complex w(1, 0); // 计算当前单位复根
for (int k = j; k < j + h / 2; k++) {
Complex u = y[k];
Complex t = w * y[k + h / 2];
y[k] = u + t; // 这就是把两部分分治的结果加起来
y[k + h / 2] = u - t;
// 后半个 “step” 中的ω一定和 “前半个” 中的成相反数
// “红圈”上的点转一整圈“转回来”，转半圈正好转成相反数
// 一个数相反数的平方与这个数自身的平方相等
w = w * wn;
}
}
}
if (on == -1) {
for (int i = 0; i < len; i++) {
y[i].x /= len;
}
}
}
大概就是这些了，别的学不会了
