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;
}