Appearance
在 快速傅里叶变换中,FFT 算法使用分治思想,每次都将系数数组
每次拆分时,都要把偶数项拷贝到新数组
如果我们一开始就将
并自下而上地进行运算
这样就能避免愚蠢的拷贝操作,节省计算机的算力。
这个优化过程被称为蝶形优化。
蝶形优化
简言之,我们需要的是这样的变换:
这个变换的策略出人意料地简单:将每个数转为二进制(在最高位补零补到一样长),然后反转它们,就能得到蝶形优化的结果。
举个例子:
是不是很神奇。
证明
我们仅保留数组
这张图所展示的递归策略可以被简要地概括为:
在每一层,如果数
需要特别强调「在其序列中」这个点。例如在第二层
其中数字
问题来了:在某一层,如何判断数
首先我们需要归纳
容易看出:
- 在第一层:数
排在第 位; - 在第二层:数
排在第 位; - 在第三层:数
排在第 位; - ...
- 在第
层,数 排在第 位。
我们直接将
可能你已经发现了
也就是说,只需判断
考察数字
- 在第一层,
的第一位是 ,说明它在奇位,因此被划分到右侧; - 在第二层,
的第二位是 ,说明它在偶位,因此被划分到左侧; - 在第三层,
的第三位是 ,说明它在偶位,因此被划分到左侧。
思考:数
进一步地,我们可以归纳出:
- 在第一层:
- 被划分到左侧
最高位为 ; - 被划分到右侧
最高位为 ;
- 被划分到左侧
- 在第二层:
- 被划分到左侧
次高位为 ; - 被划分到右侧
次高位为 ;
- 被划分到左侧
- ...
现在重新审视数字
- 在第一层,
的第一位是 ,被划分到右侧,变换后它的位序的最高位是 ; - 在第二层,
的第二位是 ,被划分到左侧,变换后它的位序的次高位是 ; - 在第三层,
的第三位是 ,被划分到左侧,变换后它的位序的第三高位是 。
即
同理,
那么,原先是
证毕。
模板
cpp
#include <bits/stdc++.h>
using namespace std;
const double PI = acos(-1);
typedef complex<double> Comp;
int reverseBits(int n, int log2n) {
int reversed = 0;
for (int i = 0; i < log2n; i++) {
if (n & (1 << i)) {
reversed |= 1 << (log2n - 1 - i);
}
}
return reversed;
}
void bit_reverse_swap(vector<Comp>& a) {
int n = a.size();
int log2n = 0;
while ((1 << log2n) < n) log2n++;
for (int i = 0; i < n; i++) {
int reversed = reverseBits(i, log2n);
if (i < reversed) {
swap(a[i], a[reversed]);
}
}
}
void FFT(vector<Comp>& a, bool invert) {
int n = a.size();
bit_reverse_swap(a);
for (int len = 2; len <= n; len <<= 1) {
double theta = 2 * PI / len * (invert ? -1 : 1);
Comp wn(cos(theta), sin(theta));
for (int i = 0; i < n; i += len) {
Comp w(1);
for (int j = 0; j < len / 2; ++j) {
Comp u = a[i + j];
Comp v = a[i + j + len / 2] * w;
a[i + j] = u + v;
a[i + j + len / 2] = u - v;
w *= wn;
}
}
}
if (invert) {
for (Comp& x : a)
x /= n;
}
}
vector<Comp> multiply(vector<Comp> A, vector<Comp> B) {
int n = 1;
while (n < A.size() + B.size())
n *= 2;
A.resize(n);
B.resize(n);
FFT(A, false);
FFT(B, false);
vector<Comp> C(n);
for (int i = 0; i < n; i++)
C[i] = A[i] * B[i];
FFT(C, true);
for (int i = 0; i < n; i ++)
C[i] = round(C[i].real());
while (C.size() && ! C.back().real())
C.pop_back();
return C;
}
int main() {
vector<Comp> A = {1, 2, 3}; // Represents the polynomial 1 + 2x + 3x^2
vector<Comp> B = {4, 5}; // Represents the polynomial 4 + 5x
vector<Comp> C = multiply(A, B);
for (auto i : C)
cout << i.real() << ' ';
cout << endl;
return 0;
}