简介
FFT,快速傅里叶变换,用于快速求多项式。首先,我们知道任何一个任何一个周期函数都可以表示成为傅里叶级数的形式。通过对应关系我们可以求出周期函数系数或者傅里叶级数的系数。
建议具体过程看复变函数或者高数。
FFT过程:
两个次数界为n的多项式a(x)和b(x)相乘,输入输出均采用系数表示法。(假定n为2的幂)
1)使次数界增加一倍:a(x)和b(x)扩充为次数界为2n的多项式,并构造起系数表示。
2)求值:两次应用2n阶FFT,计算出a(x)和b(x)的长度为2n的点值表示。
3)点乘:计算多项式的点值表示。
4)插值:对2n个点值对应用一次FFT计算出其逆DFT,就可以构造出多项式c(x)的系数表示。
第1步和第3步需要时间为O(n),第2步和第4步运用FFT需要时间为O(nlogn)。
模板:
const double PI = acos(-1.0);
struct Complex{
double x, y; // 实部和虚部 x+yi
Complex(double _x = 0.0, double _y = 0.0) {
x = _x;
y = _y;
}
Complex operator-(const Complex &b) const {
return Complex(x - b.x, y - b.y);
}
Complex operator+(const Complex &b) const {
return Complex(x + b.x, y + b.y);
}
Complex operator*(const Complex &b) const {
return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
}
};
/*
* 进行 FFT 和 IFFT 前的反转变换。
* 位置 i 和 (i 二进制反转后位置)互换
* len 必须去 2 的幂
*/
void change(Complex y[], int len){
for (int i = 1, j = len / 2,k;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 -= k;
k >>= 1;
}
if (j < k) j += k;
}
}
/*
* 做 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(-on * 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;
w = w * wn;
}
}
}
if (on == -1)
for (int i = 0; i < len; i++) y[i].x /= len;
}
void Conv(Complex a[],Complex b[],int len) //求卷积
{
fft(a,len,1);
fft(b,len,1);
for(int i = 0;i < len;++i) a[i] = a[i] * b[i];
fft(a,len,-1);
}