FFT简陋入门
FFT入门
FFT的用途
在\(\,\Theta(n\log{n})\,\)的时间内计算离散傅里叶变化（DFT），通常用来计算多项式乘法
点值表达式
引理1：任何\(\,n-1\,\)次多项式可以由其在\(\,n\,\)个点的取值唯一确定
考虑反证，设\(\,n\,\)个点\(\,a_1,a_2\cdots a_n\)同时被两个\(\,n-1\,\)次多项式函数\(\,A,B\,\)经过，令
与代数基本定理矛盾，假设不成立，得证
那么，可以用点值表达式唯一确定一个多项式，而点值表达式的乘积可以在\(\,\Theta(n)\,\)的时间内算出，就是对应点的点值相乘
关于复数
如果您学过复数，请略过
我们定义复数为形如\(\,a+bi\,\)的数，其中\(\,a\,\)叫实部，\(\,b\,\)叫虚部
定义复数的运算：
不难发现这些运算符合交换律，结合律和分配律
定义复数的模长\(\,|z|\,\)和辐角\(\,\theta\,\)：
定义\(\,z\,\)的共轭复数\(\,z^\prime\,\)：
结合定义，有
把复数\(\,z=a+bi\,\)映射为平面上的点\(\,(a,b)\,\)，可以发现复数\(\,a+bi\,\)和向量\(\,(a,b)\,\)一一对应
考虑两者之间的关联，发现加减规则相同，向量内积的几何意义是平行四边形的对角线长度，而复数乘意为模长相乘，辐角相加，这一点结合复数的运算不难推出
单位根
定义\(\,n\,\)次单位根\(\,\omega\,\)为满足\(\,\omega^{n}=1\,\)的复数\(\,\omega\)，发现\(\,n\,\)次单位根对应平面上单位圆的\(\,n\,\)等分点
设\(\,\omega_{n}^{k}\,\)表示辐角第\(\,k\,\)大的\(\,n\,\)次单位根，联想复数乘法的几何意义，不难发现：
运用欧拉公式\(\,e^{i\theta}=isin\theta+cos\theta\,\)，可得\(w_n^k=e^{2\pi\frac{k}{n}i}\)，由此可以得到下方的两条引理：
引理2（消去引理）：\(w_{mn}^{mk}=w_{n}^k\).
引理3（折半引理）：\(w_n^{k+\frac{n}{2}}=-w_n^k\).
引理2证明如下：
引理3证明如下：
单位根与FFT
考虑到FFT有两个难点，一是系数表达式转点值表达式，二是反向转化
正向变换
先考虑用单位根优化第一步，令偶次多项式\(\,A,B\,\)为
这里如果次数不足，系数用\(\,0\,\)补足
不妨讨论\(\,A\,\)，令
于是有：
将\(\,n\,\)次单位根一一带入\(\,A\,\)有
若\(\,0\leq k \leq \frac{n}{2}-1\)，则
若\(\,\frac{n}{2}\leq k^\prime \leq n-1\)，记\(\,k^\prime=k+\frac{n}{2}\,\)则
显然，知道\(\,A1,A2\,\)的值之后，我们可以\(\,\Theta(1)\,\)求出\(\,A\,\)的取值，这样操作\(\,n\,\)次，我们就得到了\(\,A\,\)的点值表达式
\(A1,A2\)是两个\(\frac{n}{2}\)次的多项式，可以递归求解，于是复杂度为：
反向变换
再考虑将点值表达式转化为系数表达式，这个过程叫做离散傅里叶反变换
不妨记原多项式函数为\(\,F(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}=\sum_{i=0}^{n-1}a_ix^i\)
令\(\,A(x)=\sum_{i=0}^{n-1}d_ix^i\,\)其中\(\,d\,\)是\(\,a\,\)离散傅里叶变换的结果，即我们对\(\,a\,\)进行正向变换从而得到\(\,d\,\)
设
考虑将\(\,d\,\)展开，得
为了利用单位根的性质，考虑交换求和号
令
不难看出这是一个等比数列，而且有\(S(i,i)=n\)，接着考虑\(j\neq k\)的情况
直接运用等比数列求和公式：
可以注意到上式中\(\,w_n^n\equiv1\,\)，那么便不难计算
整理得：
将\(\,S(j,k)\,\)带回\(\,c_k\,\)的表达式，得
考虑用\(\,c\,\)来表达\(F(x)\)，于是有
此时就得到了大致的思路：
正向变换时已经得到了\(\,d\)，就是点值两两相乘的结果
由\(c_k=A(w_n^{-k})=\sum_{i=0}^{n-1}d_i(w_n^{-k})^i\)不难看出\(\,c\,\)是\(\,d\,\)进行一次离散傅里叶变换的结果
最后有
综上，我们实现了在\(\,\Theta(n\log n)\,\)的时间内计算多项式乘法
代码实现
class Complex
{
public:
double r,i;
Complex()=default;
Complex(double R,double I):r(R),i(I){}
friend Complex operator + (Complex a,Complex b){return Complex(a.r+b.r,a.i+b.i);}
friend Complex operator - (Complex a,Complex b){return Complex(a.r-b.r,a.i-b.i);}
friend Complex operator * (Complex a,Complex b){return Complex(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
friend Complex operator / (Complex a,double b){return Complex(a.r/b,a.i/b);}
friend void operator /=(Complex &a,double b){a=a/b;}
friend Complex operator ^ (Complex a,int b)
{
Complex c(1,0);
while(b) c=(b&1?c*a:c),a=a*a,b>>=1;
return c;
}
friend void operator ^=(Complex &a,int b){a=a^b;}
};
Complex res[maxn];
int r[maxn];int len;
int x,y,lim;
const double PI=acos(-1),eps=1e-9;
inline void clear() { for(int i=0;i<lim;++i) res[i]=Complex(),r[i]=0; }
inline void calc() { for(int i=1;i<lim;i++) r[i]=r[i>>1]>>1|(i&1)<<len; }
inline void fft(Complex *tp,int n,int op)
{
for(int i=1;i<n;++i) if(i<r[i]) swap(tp[i],tp[r[i]]);
for(int k=1;k<n;k<<=1)
for(int s=0;s<n;s+=k<<1)
{
Complex now(1,0),rt(cos(PI/k),op*sin(PI/k)),ev,od;
for(int i=s;i<s+k;++i,now=now*rt)
ev=tp[i],
od=tp[i+k],
tp[i]=ev+now*od,
tp[i+k]=ev-now*od;
}
if(op==-1) for(int i=0;i<n;i++) tp[i]/=n;
}
