数学——快速傅里叶变换(FFT)

Shan xizeng

基础知识

快速傅里叶变换，用来求出两个多项式相乘，如果暴力相乘，时间复杂度为\(O(n^2 )\)，使用快速傅里叶变换，可以优化到\(O(n log n)\)。

准备知识：

多项式：设\(A(x)\)表示一个n-1次多项式，则\(A(x)=a_0x^0+a_1x^1+\cdots+a_nx^n\)

多项式的表示方法：

一是用系数表示法，表示为\(\sum_{i=0}^na_ix^i\)，二是点值表示法，表示为对于几个具体的\(x\)对应的\(A(x)\)的值，最少需要\(n\)个不同的点就能表示唯一一个\(n-1\)次多项式。

多项式运算：

加法：

如果使用系数表示法，则将各个系数相加，复杂度为\(O(n)\)；

如果使用点值表示法，则将横坐标相同点的纵坐标相加，复杂度相同。

乘法：

如果使用系数表示法，则设得到的多项式为\(\sum_{i=0}^nc_ix^i\)，其中，\(c_i=\sum_{j+k=i,0\leq j,k\leq n}a_jb_k\)，很显然，时间复杂度为\(O(n^2)\)；

如果使用点值表示法，则将横坐标相同点的纵坐标相乘，复杂度仍不变，为\(O(n)\)。

向量：

物理、几何学意义：同时具有大小和方向的量。向量相加满足平行四边形定则。

复数：

分为实部与虚部，形如\(a+bi\),。

复数单位根：

将复数\(a+bi\)中a看做横坐标，b看做纵坐标，则复数单位根为到原点距离为1的点。这些点构成的集合为一个圆，称为单位圆。单位圆的n等分点称为n次单位根，将幅角为正且最小的数设为\(\omega_n\)，则n次单位根分别为\(\omega_n^2,\omega_n^3,\dots,\omega_n^n\)。\(\omega_n=\omega_n^n=1\)。根据欧拉公式\(e^{\theta i}=\cos \theta+i\ \sin \theta\)，\(\omega_n^k=e^{\frac{2\pi k i}{n}}=\cos k \cdot \frac{2\pi}{n}+i\sin k \cdot \frac{2\pi}{n}\)。

易知，\(\omega_n^{k+\frac{n}{2}}=-\omega_n^k,\omega_n^2=\omega_{\frac{n}{2}}\)。(可以画图试试)

快速傅里叶变换

快速傅里叶变换的思想主要是分治。

进行快速傅里叶变换，就是将n个n次单位根分别对多项式求出对应的值。

对于多项式A(x),将其n次单位根带入，则可得到\(A(\omega_n^k)=\sum_{i=0}^na_i\omega_n^{ki}\)

我们将其按照奇偶性分组，得到：\(A(\omega_n^k)=\sum_{i=0}^{n/2-1}a_{2i}\omega_n^{2ki}+\omega_n^k\sum_{i=0}^{n/2-1}a_{2i+1}\omega_n^{2ki}\)

由上面的公式得：\(A(\omega_n^k)=\sum_{i=0}^{n/2-1}a_{2i}\omega_{n/2}^{ki}+\omega_n^k\sum_{i=0}^{n/2-1}a_{2i+1}\omega_{n/2}^{ki}\)

并且

\[\begin{eqnarray*} A(\omega_n^{k+\frac{n}{2}}) &=& \sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^{k+\frac{n}{2}}\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki} \\ &=&\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}-\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki} \end{eqnarray*}\]

这样需要带入的值就减少了一半，时间复杂度为\(T(n)=2T(n/2)+O(n)=O(n\log n)\)。

当然，我们还需要将点值表示转化为系数表示，具体实现就是傅里叶逆变换，就是把原来的傅里叶变换里的\(\omega_n^k\)变成\(-\omega_n^k\)带上然后把系数除以n就好了。证明：~~我不会。。~~

实现的快速方法：二进制优化，背板子就好了QAQ

#include
    
     #include
     
      #include
      
       #include
       
        using namespace std;const int Maxn=1100000;const double Pi=3.14159265358979323846;int n,m,r[Maxn],limit=1,l;struct complex {    double x,y;}a[Maxn],b[Maxn];complex operator + (complex a,complex b) {    return (complex) {a.x+b.x,a.y+b.y};}complex operator - (complex a,complex b) {    return (complex) {a.x-b.x,a.y-b.y};}complex operator * (complex a,complex b) {    return (complex) {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}void fft(complex *a,int type) {    for(int i=0;i
        
         >1]>>1|((i&1)<