FFT(快速傅里葉變換)

用途: 計算多項式乘積。

\(A(x)=a_{0}+a_{1}x+a_{2}x^{2}+···+a_{n-1}x^{n-1}\)

\(B(x)=b_{0}+b_{1}x+b_{2}x^{2}+···+b_{n-1}x^{n-1}\)

求 : \(A(x)B(x)\)

前置芝士

多項式性質: 用任意 \(n\) 個多項式函數圖像上的點，可惟一肯定一個 \(n-1\) 次多項式。

複數,不會的話右轉數學選修2-2.

作法

將 \(A(x)B(x)\) 的係數表示法轉化爲點表示法，

即 \((x_{1},A(x_{1}))\ \ \ (x_{2},A(x_{2}))\ ···\ (x_{n},A(x_{n}))\)
和 \((x_{1},B(x_{1}))\ \ \ (x_{2},B(x_{2}))\ ···\ (x_{n},B(x_{n}))\)

那麼 \(A(x)B(x)\) 就能夠表示爲

\[(x_{1},A(x_{1})B(x_{1}))\ \ \ (x_{2},A(x_{2})B(x_{2}))\ ···\ (x_{n},A(x_{n})B(x_{n})) \]

再將其轉化爲係數表示法。

如何轉化？

點的橫座標選擇複數域上的單位根。

將複平面中以原點爲圓心的單位元等分爲 \(n\) 份，以 \(\omega _{n}^{k}\) 表示從 \(x\) 正半軸一次逆時針取 \(n\) 份中的 \(k\) 份， \(\omega _{n}^{k}\) 就被稱爲 \(n\) 次單位根。

\(n\) 次單位根的性質:

• \(\omega ^{i}_{n} \ne \omega^{j}_{n} \ \ \ \forall i\ne j\)
• \(\omega^{k}_{n}=\cos \frac{2k\pi}{n}+\sin \frac{2k\pi}{n}i\)
• \(\omega^{0}_{n}=\omega^{n}_{n}=1\)
• \(\omega^{2k}_{2n}=\omega^{k}_{n}\)
• \(\omega^{k+\frac{n}{2}}_{n}=-\omega^{k}_{n}\)

那麼點表示法選取的橫座標爲: \(\omega^{0}_{n}\ \ \omega^{1}_{n}\ \ \omega^{2}_{n}···\omega^{n-1}_{n}\)

如今，須要快速計算:

\((\omega^{0}_{n},A(\omega^{0}_{n})\ \ \ (\omega^{1}_{n},A(\omega^{1}_{n}))\ ···\ (\omega^{n-1}_{n},A(\omega^{n-1}_{n}))\)

快速傅里葉正變換(係數表示法轉化爲點表示法)

\[\begin{aligned} A(x)&=a_{0}+a_{1}x+a_{2}x^{2}+···+a_{n-1}x^{n-1}\ \ \ n爲2的整次冪\\ &=(a_{0}+a_{2}x^{2}+a_{4}x^{4}+···+a_{n-2}x^{n-2})+(a_{1}x+a_{3}x^{3}+a_{5}x^{5}+···a_{n-1}x^{n-1})\ \ \ 奇偶次項分開 \end{aligned}\]

設:

\[A_{1}(x)=a_{0}+a_{2}x+a_{4}x^{2}+···+a_{n-2}x^{\frac{n}{2}-1} \]

\[A_{2}(x)=a_{1}+a_{3}x+a_{5}x^{2}+···+a_{n-1}x^{\frac{n}{2}-1} \]

那麼有:

\[A(x)=A_{1}(x^2)+xA_{2}(x^2) \]

\(k\in [0,n-1]\) , 將值域分爲兩段 \([0,\frac{n}{2}-1]\ ,\ [\frac{n}{2},n-1]\)

若是 \(k\in[0,\frac{n}{2}-1]\) , 那麼有 \((\omega^{k}_{n})^{2}=\omega^{2k}_{n}\)

\[\begin{aligned} A(\omega^{k}_{n})&=A_{1}(\omega^{2k}_{n})+\omega^{k}_{n}A_{2}(\omega^{2k}_{n})\\ &=A_{1}(\omega^{k}_{\frac{n}{2}})+\omega^{k}_{n}A_{2}(\omega^{k}_{\frac{n}{2}})\\ \end{aligned}\]

若是 \([\frac{n}{2},n-1]\) , 將這一段所有減去 \(\frac{n}{2}\) , \(k\) 仍是遍歷 \([0,\frac{n}{2}-1]\) 將 \(k\) 加上 \(\frac{n}{2}\) 就好了。

\[\begin{aligned} A(\omega^{k+\frac{n}{2}}_{n})&=A_{1}(\omega^{2k+n}_{n})+\omega^{k+\frac{n}{2}}_{n}A_{2}(\omega^{2k+n}_{n})\\ &=A_{1}(\omega^{2k}_{n})-\omega^{k}_{n}A_{2}(\omega^{2k}_{n})\\ &=A_{1}(\omega^{k}_{\frac{n}{2}})-\omega^{k}_{n}A_{2}(\omega^{k}_{\frac{n}{2}}) \end{aligned}\]

發現 \(A(\omega^{k}_{n})\ ,\ A(\omega^{k+\frac{n}{2}}_{n})\) 計算很是類似，只要求出 \(A_{1}(\omega^{k}_{\frac{n}{2}})\ ,\ A_{2}(\omega^{k}_{\frac{n}{2}})\) 便可。

這個過程能夠迭代 ，每次將區間分爲兩段 ，最多會遞歸 \(log\ n\) 層 ，總複雜度 \(O(nlog\ n)\) .

快速傅里葉逆變換(點表示法轉化爲係數表示法)

有點 \((\omega^{1}_{n},A(\omega^{1}_{n}))\ \ (\omega^{2}_{n},A(\omega^{2}_{n}))···(\omega^{n-1}_{n},A(\omega^{n-1}_{n}))\).

求 \(A(x)=a_{0}+a_{1}x+a_{2}x^2+···+a_{n-1}x^{n-1}\)

設 \(y_{k}=\omega^{k}_{n}\) 則有

\[c_{k}=\sum_{i=0}^{n-1}y_{i}(\omega^{-k}_{n})^{i}=na_{k} \]

證實:

\[\begin{aligned} c_{i}&=\sum_{j=0}^{n-1}y_{j}(\omega^{-i}_{n})^{j}\\ &=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_{j}(\omega^{i}_{n})^j))(\omega^{-k}_{n})^i\\ &=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_{j}(\omega^{j}_{n})^i(\omega^{-k}_{n})^i\\ &=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_{j}(\omega^{j-k}_{n})^i\\ &=\sum_{j=0}^{n-1}a_{j}(\sum_{i=0}^{n-1}(\omega^{j-k}_{n})^i) \end{aligned}\]

設 \(S(x)=1+x+x^2+···+x^{n-1}\) .

\(k\ne 0\) 時 , 有：

\[S(\omega^{k}_{n})=\omega^{0}_{n}+\omega^{k}_{n}+\omega^{2k}_{n}+···+\omega^{(n-1)k}_{n} \]

\[\omega^{k}_{n}S(\omega^{k}_{n})=\omega^{k}_{n}+\omega^{2k}_{n}+\omega^{3k}_{n}+···+\omega^{0}_{n} \]

兩柿相減,有：

\[(1-\omega^{k}_{n})S(\omega^{k}_{n})=0 \]

故 \(S(\omega^{k}_{n})=0\)

\(k=0\) 時 , \(S(\omega^{k}_{n})=n\) ,代入顯然

綜上 : 把 \(S\) 代入柿子可知：

\[c_{i}=na_{i} \]

證畢.

設 \(S(x)=1+x+x^2+···+x^{n-1}\) .

設 \(B(x)=y_{0}+y_{1}x^1+y_{2}x^2+···+y_{n-1}x^{n-1}\).

則有 \(c_{i}=B(\omega^{-i}_{n})\) ,而後用上面說的正變換求一下就行。

FFT板子

#include<bits/stdc++.h>
#define N 300005
#define LL long long 
#define LB long double
using namespace std;

int n,m;
struct node
{
    LB x,y;
    node operator+ (const node &t) const
    {
        return (node){x+t.x,y+t.y};
    }
    node operator- (const node &t) const
    {
        return (node){x-t.x,y-t.y};
    }
    node operator* (const node &t) const
    {
        return (node){x*t.x-y*t.y,x*t.y+y*t.x};
    }
}a[N],b[N];
int rev[N],bit,tot;
const LB pi=3.141592653589793638462643383279502;

inline int qr()
{
    int x=0,w=1;char ch=0;
    while(ch<'0'||ch>'9'){if(ch=='-')w=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+(ch^48);ch=getchar();}
    return x*w;
}

void FFT(node a[],int opt)
{
    for(register int i=0;i<tot;i++)
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    for(register int mid=1;mid<tot;mid<<=1)
    {
        node w1=(node){cos(pi/mid),opt*sin(pi/mid)};
        for(register int i=0;i<tot;i+=(mid<<1))
        {
            node wk=(node){1,0};
            for(register int j=0;j<mid;j++,wk=wk*w1)
            {
                node x=a[i+j];
                node y=wk*a[i+j+mid];
                a[i+j]=x+y;
                a[i+j+mid]=x-y;
            }
        }
    }
}

int main()
{
    //freopen("data.in","r",stdin);
    //freopen("data.out","w",stdout);
    n=qr();m=qr();
    for(register int i=0;i<=n;i++)
        scanf("%Lf",&a[i].x);
    for(register int i=0;i<=m;i++)
        scanf("%Lf",&b[i].x);
    while((1<<bit)<n+m+1)
        bit++;
    tot=1<<bit;
    for(register int i=0;i<tot;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    FFT(a,1);
    FFT(b,1);
    for(register int i=0;i<tot;i++)
        a[i]=a[i]*b[i];
    FFT(a,-1);
    for(register int i=0;i<=m+n;i++)
        printf("%d ",(int)(a[i].x/tot+0.5));
    system("pause");
    return 0;
}