有用的学习链接&书籍
维基百科看英文的更多内容&有趣的图
快速傅立叶变化-百度百科,注意其中的图!
组合数学(第4版) Page 287~291(讲得挺详细)
FFT/DFT是个什么东西?
说实话,我也不知道,不过根据维基百科上面的图,就可以略窥一二了:
傅里叶变换将函数的时域(红色)与频域(蓝色)相关联。频谱中的不同成分频率在频域中以峰值形式表示。
——“Fourier transform time and frequency domains (small)”作者Lucas V. Barbosa - 自己的作品。来自维基共享资源 - http://commons.wikimedia.org/wiki/File:Fourier_transform_time_and_frequency_domains_(small).gif#mediaviewer/File:Fourier_transform_time_and_frequency_domains_(small).gif根据公有领域授权
(注,下面的话都是我乱编的,不具有很强的科学性)
它的意思是,红色的周期函数可以通过若干的波(图中蓝色部分)来表示。
FFT的作用就是通过红色的那一部分求出蓝色那一部分。在oi中几乎都是离散函数(你可以认为数组的下标就是定义域),所以oi中用的几乎都是DFT。
它的最基本函数是(摘自维基百科):
它们的过程十分相似,我们只要会了DFT,就能套入IDFT了。
注:\(e^{i\theta} = cos \theta + i sin \theta\),没错,这是复数运算,c++中有complex<double>
。
如果我们按照上面的公式计算,时间复杂度为\(O(n^2)\)。
加快计算
令\(\omega = e^{2 \pi i}\),\(\omega_n^p= e^{ \frac {2 \pi i \cdot p} {n}}\)。
组合数学一书中通过举例子,发现了计算规律(通过化简式子,比如\(\omega^6_4=\omega^2_4\)),发现规律的过程就不展开了。
下面的图片摘自百度百科:
最左边的\(8\)个黑点是原数列,最右边的\(8\)个是输出,没错,输出的下标是有顺序的,但是原数列是无序的(但有规律,看第\(4\)行,左边是\(6=1100_2\),右边是\(3=0011_2\),二进制的字符顺序恰好是反的)。
计算过程可以从代码中看出:
C表示,\(1\)为DFT,\(-1\)为IDFT,MyType
即complex<double>
,cnt
是数列大小,unit
是\(\omega_{cnt}\)。
void FFT(MyType A[], int cnt, int C) {
if (cnt == 1) return;
MyType* L = newMemory(cnt >> 1);
MyType* R = newMemory(cnt >> 1);
for (int i = 0; i < cnt; i += 2) {
L[i >> 1] = A[i];
R[i >> 1] = A[i + 1];
}
FFT(L, cnt >> 1, C);
FFT(R, cnt >> 1, C);
MyType w(1, 0);
MyType unit(cos(C * PI * 2 / cnt), sin(PI * 2 / cnt * C));
//also unit(cos( PI * 2 / cnt), sin(PI * 2 / cnt ) * C)
for (int i = 0; i < cnt >> 1; i ++) {
A[i] = L[i] + w * R[i];
A[i + (cnt >> 1)] = L[i] - w * R[i];
w *= unit;
}
}
时间复杂度\(O(n * logn)\)。
oi里用DFT来求两个多项式的乘积
算法(此处仅有算法,我不知道它的正确性如何证明):
对于两个多项式\(A\)和\(B\),我们要算出结果\(C\):
- 对\(A\)进行DFT得到\(\hat{A}\)
- 对\(B\)进行DFT得到\(\hat{B}\)
- 令\(\hat C_i = \hat A_i * \hat B_i\)
- 对\(\hat C\)进行逆DFT(IDFT)得到\(C\)
后
这篇随笔写得不清不楚的,其实我是想我自己以后看了自己明白。
听说可以写模意义下的DFT(没有误差),找个时间看看??
lyl所说的“总结”图:
附上求两个多项式乘积的代码,注complex的两个值调用分别是real和imag,看这里:
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <complex>
#include <assert.h>
using namespace std;
const int MAXLOGN = 18;
const int MAXN = 1 << MAXLOGN;
const int MAXMEMORY = MAXLOGN * MAXN;
const double PI = acos(-1.);
typedef complex<double> MyType;
int n;
MyType mem[MAXMEMORY];
MyType* cur_mem;
MyType* newMemory(int size) {
cur_mem += size;
return cur_mem - size;
}
void FFT(MyType A[], int cnt, int C) {
if (cnt == 1) return;
MyType* L = newMemory(cnt >> 1);
MyType* R = newMemory(cnt >> 1);
for (int i = 0; i < cnt; i += 2) {
L[i >> 1] = A[i];
R[i >> 1] = A[i + 1];
}
FFT(L, cnt >> 1, C);
FFT(R, cnt >> 1, C);
MyType w(1, 0);
MyType unit(cos(C * PI * 2 / cnt), sin(PI * 2 / cnt * C));
for (int i = 0; i < cnt >> 1; i ++) {
A[i] = L[i] + w * R[i];
A[i + (cnt >> 1)] = L[i] - w * R[i];
w *= unit;
}
}
int main() {
static MyType A[MAXN], B[MAXN];
int cntA, cntB;
scanf("%d%d", &cntA, &cntB);
cntA ++, cntB ++;
for (int i = 0; i < cntA; i ++)
scanf("%lf", &A[i].real());
for (int i = 0; i < cntB; i ++)
scanf("%lf", &B[i].real());
n = 1;
while (n < cntA + cntB - 1) n <<= 1;
cur_mem = mem; FFT(A, n, 1);
cur_mem = mem; FFT(B, n, 1);
for (int i = 0; i < n; i ++) A[i] *= B[i];
cur_mem = mem; FFT(A, n, -1);
for (int i = 0; i < cntA + cntB - 1; i ++)
printf("%d ", (int) floor(A[i].real() / n + 0.1));
printf("\n");
return 0;
}