快速傅裏葉變換 學習筆記

• 前言：這篇學習筆記以\(\text{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 \theta,i\sin \theta)\)，我們可以稍稍證明一下幅角相加了。

    設有負數\((\cos \theta_1,i\sin \theta_1)\),\((\cos \theta_2,i\sin \theta_2)\)，那麽它們相乘有

    \(\cos \theta_1\cos \theta_2-\sin\theta_1\sin\theta_2+i(\cos\theta_1\sin\theta_2+\sin\theta_1\cos\theta_2)\)

    \(=\cos(\theta_1+\theta_2)+i\sin(\theta_1+\theta_2)\)

  • 將復平面的單位圓劃分成\(n\)份，定義第\(k\)份的向量對應的復數為\(w_n^k\)。特殊的\(w_n^0=w_n^n=1\)。

    為了方便和必要，以下的\(n\)均認為是\(2\)的正整次數冪。

  • 歐拉公式：

    \(w_n^k\)的坐標表示為\((\cos\frac{2k\pi}{n},i\sin\frac{2k\pi}{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];
    }
}

快速傅裏葉逆變換

顯然我們可以有
\(\begin{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}\times \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\a_{n-1}\\ \end{bmatrix}=\begin{bmatrix} y_0 \\ y_1 \\ \vdots \\y_{n-1}\\\end{bmatrix}\)

現在我們有\(y\)和\(w\)那些東西，要求\(a\)，按道理我們要找到一個\(w\)那個的逆矩陣然後乘在右邊，然後我們有結論

\(\begin{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}\times \begin{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\)優化

• 遞歸版一看就好慢的。

如上圖，

然後我們發現二進制表示反過來了。

於是我們可以得到一個轉換方法。

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]

快速傅裏葉變換 學習筆記