快速傅里叶变换 学习笔记

快速傅里叶变换 学习笔记

• 前言：这篇学习笔记以( ext{Dew})的意会yy为主，有些地方会比较简略，不过该有的证明应该还是会有的。

多项式的表示法

• 系数表示法

  (f(x)=a_0+a_1x+a_2x^2+dots+a_nx^n)

• 点值表示法

  (n+1)个不同的点可以确定一个(n)次的多项式。

  即一个(n)次多项式可以被((x_1,y_1),(x_2,y_2),dots,(x_n,y_n))全部表示出来

• 考虑如何求出点值表示法的所有点

复数与单位根

• 定义复数运算

  • (a+bi+c+di=a+c+(b+d)i)

    (a+bi-(c+di)=a-c+(b-d)i)

    几何定义同向量

  • ((a+bi)(c+di)=ac-bd+(ad+bc)i)

    几何定义：模长相乘，幅角相加

  struct complex
  {
       double x,y;
       complex(){}
       complex(double x,double y){this->x=x,this->y=y;}
       complex friend operator +(complex n1,complex n2){return complex(n1.x+n2.x,n1.y+n2.y);}
       complex friend operator -(complex n1,complex n2){return complex(n1.x-n2.x,n1.y-n2.y);}
       complex friend operator *(complex n1,complex n2){return complex(n1.x*n2.x-n1.y*n2.y,n1.x*n2.y+n1.y*n2.x);}
  };
  

• 单位根

  • 按照乘法的几何定义，我们发现在模长为(1)的单位圆上，两个负数相乘的模长不变，于是考虑单位圆上的表示法。

  • 在半径为(1)的单位圆上，向量的坐标为((cos heta,isin heta))，我们可以稍稍证明一下幅角相加了。

    设有负数((cos heta_1,isin heta_1)),((cos heta_2,isin heta_2))，那么它们相乘有

    (cos heta_1cos heta_2-sin heta_1sin heta_2+i(cos heta_1sin heta_2+sin heta_1cos heta_2))

    (=cos( heta_1+ heta_2)+isin( heta_1+ heta_2))

  • 将复平面的单位圆划分成(n)份，定义第(k)份的向量对应的复数为(w_n^k)。特殊的(w_n^0=w_n^n=1)。

    为了方便和必要，以下的(n)均认为是(2)的正整次数幂。

  • 欧拉公式：

    (w_n^k)的坐标表示为((cosfrac{2kpi}{n},isinfrac{2kpi}{n}))

  • 性质

    1. ((w_n^k)^2=w_n^{2k}=w_{frac{n}{2}}^k)
    2. (w_n^{k+frac{n}{2}}=-w_n^k)

快速傅里叶变换

对多项式(F(x)=a_0+a_1x+a_2x^2+dots+a_{n-1}x^{n-1})

(F(x)=a_0+a_2x^2+dots+a_{n-2}x^{n-2}+x(a_1+a_3x^2+dots a_{n-1}x^{n-2}))

设(FL(x)=a_0+a_2x+dots+a_{n-2}x^{frac{n}{2}-1}),(FR(x)=a_1+a_3x+dots+a_{n-1}x^{frac{n}{2}-1})

则(F(x)=FL(x^2)+xFR(x^2))

取(k<frac{n}{2})，并代入(x=w_n^k)

(F(w_n^k)=FL(w^k_{frac{n}{2}})+w_n^kFR(w^k_{frac{n}{2}}))

(F(w_n^{k+frac{n}{2}})=FL(w^k_frac{n}{2})-w_n^kFR(w_frac{n}{2}^k))

显然可以子问题求解，于是我们得到了(O(nlogn))求点值表示法的方法了。

void FFT(int len,int *a,int typ)
{
    if(len==1) return;
    complex a1[len>>1],a2[len>>1];
    for(int i=0;i<n;i+=2)
        a1[i>>1]=a[i],a2[i>>1]=a[i+1];
    FFT(len>>1,a1,typ);
    FFT(len>>1,a2,typ);
    complex wn=complex(cos(2*pi/len),typ*sin(2*pi/len)),w=(1,0);
    for(int i=0;i<len>>1;i++.w=w*wn)
    {
        a[i]=a1[i]+w*a2[i];
        a[i+(len>>1)]=a1[i]-w*a2[i];
    }
}

快速傅里叶逆变换

显然我们可以有
(egin{bmatrix}1&1&1&cdots&1\1 &w_n^1&w_n^2&cdots&w_n^{n-1}\vdots& vdots &vdots&ddots&vdots\1& w_n^{n-1}&w_n^{2(n-1)}&cdots&w_n^{(n-1)(n-1)}\end{bmatrix} imes egin{bmatrix} a_0 \ a_1 \ vdots \a_{n-1}\ end{bmatrix}=egin{bmatrix} y_0 \ y_1 \ vdots \y_{n-1}\end{bmatrix})

现在我们有(y)和(w)那些东西，要求(a)，按道理我们要找到一个(w)那个的逆矩阵然后乘在右边，然后我们有结论

(egin{bmatrix}1&1&1&cdots&1\1 &w_n^1&w_n^2&cdots&w_n^{n-1}\vdots& vdots &vdots&ddots&vdots\1& w_n^{n-1}&w_n^{2(n-1)}&cdots&w_n^{(n-1)(n-1)}\end{bmatrix} imes frac{1}{n}egin{bmatrix}1&1&1&cdots&1\1&w_n^{-1}&w_n^{-2}&cdots&w_n^{-(n-1)}\vdots &vdots&vdots&ddots&vdots\1& w_n^{-(n-1)}&w_n^{-2(n-1)}&cdots&w_n^{-(n-1)(n-1)}\end{bmatrix}=I)

这个结合前面单位根的性质很容易意会到的。

然后发现这个也可以看做一个系数表示法转点值表示法的过程，而且和前面的系数矩阵很像，所以只需要改个参数就搞到啦。

(Code:)(递归版)

#include <cstdio>
#include <cmath>
const int N=(1<<21)+10;
struct complex
{
     double x,y;
     complex(){}
     complex(double x,double y){this->x=x,this->y=y;}
     complex friend operator +(complex n1,complex n2){return complex(n1.x+n2.x,n1.y+n2.y);}
     complex friend operator -(complex n1,complex n2){return complex(n1.x-n2.x,n1.y-n2.y);}
     complex friend operator *(complex n1,complex n2){return complex(n1.x*n2.x-n1.y*n2.y,n1.x*n2.y+n1.y*n2.x);}
}a[N],b[N];
int n,m;
const double pi=3.1415926535897;
void FFT(int len,complex *a,int typ)
{
    if(len==1) return;
    complex a1[len>>1],a2[len>>1];
    for(int i=0;i<len;i+=2)
        a1[i>>1]=a[i],a2[i>>1]=a[i+1];
    FFT(len>>1,a1,typ);
    FFT(len>>1,a2,typ);
    complex wn=complex(cos(2*pi/len),typ*sin(2*pi/len)),w=complex(1,0);
    for(int i=0;i<len>>1;i++,w=w*wn)
    {
        a[i]=a1[i]+w*a2[i];
        a[i+(len>>1)]=a1[i]-w*a2[i];
    }
}
int main()
{
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++) scanf("%lf",&a[i].x);
    for(int i=0;i<=m;i++) scanf("%lf",&b[i].x);
    int len=1;while(len<=n+m) len<<=1;
    FFT(len,a,1),FFT(len,b,1);
    for(int i=0;i<=len;i++) a[i]=a[i]*b[i];
    FFT(len,a,-1);
    for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/len+0.5));
    return 0;
}

(Butterfly)优化

• 递归版一看就好慢的。

如上图，

原数列数字	0	1	2	3	4	5	6	7
二进制	000	001	010	011	100	101	110	111

最底层数列数字	0	4	2	6	1	5	3	7
二进制	000	100	010	110	001	101	011	111

然后我们发现二进制表示反过来了。

于是我们可以得到一个转换方法。

for(int i=0;i<len;i++) turn[i]=turn[i>>1]>>1|((i&1)<<L);

然后直接从下往上做就不需要进行递归啦

Code:(非递归版)

#include <cstdio>
#include <algorithm>
#include <cmath>
const int N=(1<<21)+10;
struct complex
{
     double x,y;
     complex(){}
     complex(double x,double y){this->x=x,this->y=y;}
     complex friend operator +(complex n1,complex n2){return complex(n1.x+n2.x,n1.y+n2.y);}
     complex friend operator -(complex n1,complex n2){return complex(n1.x-n2.x,n1.y-n2.y);}
     complex friend operator *(complex n1,complex n2){return complex(n1.x*n2.x-n1.y*n2.y,n1.x*n2.y+n1.y*n2.x);}
}a[N],b[N],tmpx,tmpy,wn,w;
int n,m,len=1,L=-1,turn[N];
const double pi=3.1415926535897;
void FFT(complex *a,int typ)
{
    for(int i=0;i<len;i++)
        if(turn[i]>i) std::swap(a[i],a[turn[i]]);
    for(int le=1;le<len;le<<=1)
    {
        wn=complex(cos(pi/le),typ*sin(pi/le));
        for(int dl=le<<1,p=0;p<len;p+=dl)
        {
            w=complex(1,0);
            for(int k=0;k<le;k++,w=w*wn)
            {
                tmpx=a[p+k],tmpy=w*a[p+k+le];
                a[p+k]=tmpx+tmpy;
                a[p+k+le]=tmpx-tmpy;
            }
        }
    }
}
int main()
{
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++) scanf("%lf",&a[i].x);
    for(int i=0;i<=m;i++) scanf("%lf",&b[i].x);
    while(len<=n+m) len<<=1,++L;
    for(int i=0;i<len;i++) turn[i]=turn[i>>1]>>1|((i&1)<<L);
    FFT(a,1),FFT(b,1);
    for(int i=0;i<=len;i++) a[i]=a[i]*b[i];
    FFT(a,-1);
    for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/len+0.5));
    return 0;
}

---

参考资料：

自为风月马前卒

[Flower_pks]