struct Complex
{
double r,i;
Complex(double _r = 0,double _i = 0)
{
r = _r; i = _i;
}
Complex operator +(const Complex &b)
{
return Complex(r+b.r,i+b.i);
}
Complex operator -(const Complex &b)
{
return Complex(r-b.r,i-b.i);
}
Complex operator *(const Complex &b)
{
return Complex(r*b.r-i*b.i,r*b.i+i*b.r);
}
};
void change(Complex y[],ll len) //二进制平摊反转置换 O(logn)
{
ll i,j,k;
for(i = 1, j = len/2;i < len-1;i++)
{
if(i < j)swap(y[i],y[j]);
k = len/2; // Rader算法 二进制反转
while( j >= k)
{
j -= k;
k /= 2;
}
/*if(j < k)*/j += k;
}
}
/*
* *len必须为2^k形式,不够的话,就补0
* * on为1的话就是DFT,为-1就是IDFT
*/
void fft(Complex y[],ll len,ll on)
{
change(y,len); //调用反转置换
for(ll h = 2;h <= len;h <<= 1)
{
Complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h)); //初始化单位复根
for(ll j = 0;j < len;j += h)
{
Complex w(1,0);
for(ll 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(ll i = 0;i < len;i++)
y[i].r /= len; //IDFT
}