快速傅立叶变换(FFT)

多项式

系数表示法

设\(f(x)\)为一个\(n-1\)次多项式,则\(f(x)=\sum\limits_{i=0}^{n-1}a_i*x_i\)

其中\(a_i\)为\(f(x)\)的系数,用这种方法计算两个多项式相乘(逐位相乘)复杂度为\(O(n^2)\)

点值表示法

根据小学知识,一个\(n-1\)次多项式可以唯一地被\(n\)个点确定

即,如果我们知道了对于一个多项式的\(n\)个点\((x_1,y_1),(x_2,y_2)……(x_n,y_n)\)

那么这个多项式唯一满足,对任意\(1\le i \le n\),满足\(y_i=\sum\limits_{j=0}^{n-1}a_j*x_i^j\)

那么用点值实现多项式相乘是什么复杂度呢?

首先我们需要选\(n\)个点,每个点需要求出其在多项式中的值,复杂度为\(O(n^2)\)

然后把两个点值表示的多项式相乘,由于\(c(x_i)=a(x_i)*b(x_i)\),复杂度为\(O(n)\)

最后插值法用点值求出系数,复杂度为\(O(n^2)\)(我还不会插值)

考虑如果可以快速实现点值转插值和插值转点值,岂不是可以快速计算多项式乘法(说的简单,你倒是告诉我怎么快速转化啊)

前置芝士

复数

定义虚数单位\(i=\sqrt{-1}\),\(a,b\)为实数,则形如\(a+bi\)的数叫复数

其中\(a\)为复数的实部,\(b\)为复数的虚部

在复平面中,复数可以被表示为向量,所以和向量具有很多相似的性质,我们也可以用向量来理解复数,但是复数具有更多性质,比如作为一个数代入多项式

其中模长定义为\(\sqrt{a^2+b^2}\),幅度定义为\(x\)轴正半轴到向量转角的有向角

复数运算法则:

加减法与向量相同,重点是乘法:

几何定义为:模长相乘,幅角相加

代数定义为:
\[(a+bi)*(c+di)\]
\[=ac+adi+bci+bdi^2\]
\[=(ac-bd)+(ad+bc)i\]

单位根

我们首先定义圆心为坐标原点,半径为\(1\)的圆叫做单位圆

我们将这个圆\(n\)等分,得到\(n\)个圆上的点,以这\(n\)个圆上点的横坐标作为实部,纵坐标作为虚部,就得到了\(n\)个复数

快速傅立叶变换(FFT)

网上扒的图片,侵删\(QwQ\)

首先我们不自找麻烦,以\((1,0)\)作为这\(n\)个点的起点,记作\({\omega _n}^0\),逆时针方向第\(k\)个点记作\({\omega _n}^k\)

根据模长相乘幅角相加,我们可以看出\(\omega _n^k\)是\(\omega _n^0\)的\(k\)次方,其中\(\omega _n^1\)被称为\(n\)次单位根

根据幅角,我们可以计算出\(\omega _n^k\)表示的复数为\(cos{\frac{k}{n}2\pi}+i*sin{\frac{k}{n}2\pi}\)

单位根具有一些性质:

\(1.\omega _n^k=cos{\frac{k}{n}2\pi}+i*sin{\frac{k}{n}2\pi}=e^{\frac{2\pi ki}{n}}\)

证明:这个第一步到第二步由定义得出,第二步到第三步有欧拉公式得出

\(2.\omega _{2n}^{2k}=\omega_{n}^{k}\)

证明:\(\omega _{2n}^{2k}=e^{\frac{2\pi 2ki}{2n}}=e^{\frac{2\pi ki}{n}}=\omega_{n}^{k}\)

\(3.\omega _{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}\)

证明:\(\omega _{n}^{k+\frac{n}{2}}=\omega _{2n}^{2k+n}=\omega _{2n}^{2k}*\omega _{2n}^{n}=\omega_{n}^{k}*(cos\pi+i*sin\pi)=-\omega_{n}^{k}\)

\(4.\omega _{n}^{0}=\omega _{n}^{n}=1\)

证明:不用证了吧……

离散傅立叶变换(DFT)

假设\(f(x)=\sum\limits_{i=0}^{n-1}a_i*x_i\)

\(DFT(a)=(f(1),f(\omega _{n}),f(\omega _{n}^{2}),……,f(\omega _{n}^{n-1}))\)

通俗点说,就是对于一个系数表示法的多项式\(f(x)\),将\((1,\omega _{n},\omega _{n}^{2},……,\omega _{n}^{n-1})\)带入求出该多项式的点值表示法

离散傅立叶逆变换(IDFT)

将\(f(x)\)在\(n\)个\(n\)次单位根处的点值表示转化为系数表示

这里就可以回答,为什么我们要让\(n\)次单位根作为\(x\)代入多项式

假设\((y_0,y_1,y_2,……,y_{n-1})\)是多项式\(A(x)=\sum\limits_{i=0}^{n-1}a_i*x_i\)的离散傅立叶变换

我们另有一个多项式\(B(x)=\sum\limits_{i=0}^{n-1}=y_i*x_i\)

将上述\(n\)次单位根的倒数\((1,\omega _{n}^{-1},\omega _{n}^{-2},……,\omega _{n}^{-(n-1)})\)代入\(B(x)\)得到新的离散傅立叶变换\((z_0,z_1,z_2,……,z_{n-1})\)

则我们发现
\[z_k=\sum\limits_{i=0}^{n-1}y_i*(\omega _n^{-k})^i\]
\[=\sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}a_j*(\omega _n^{i})^j)*(\omega _n^{-k})^i\]
\[=\sum\limits_{j=0}^{n-1}a_j(\sum\limits_{i=0}^{n-1}(\omega _n^{j-k})^i)\]

对于\(\sum\limits_{i=0}^{n-1}(\omega _n^{j-k})^i\)我们单独考虑:

当\(j-k=0\)时,答案为\(n\),\(j\ne k\)时

等比数列求和得到\(\frac{(\omega _n^{j-k})^n-1}{\omega _n^{j-k}-1}=\frac{(\omega _n^n)^{j-k}-1}{\omega _n^{j-k}-1}=\frac{1^{j-k}-1}{\omega _n^{j-k}-1}=0\)

所以
\[\sum\limits_{j=0}^{n-1}a_j(\sum\limits_{i=0}^{n-1}(\omega _n^{j-k})^i)=n*a_k\]

\[a_k=\frac{z_k}{n}\]
得出结论:对于以\(A(x)\)的离散傅立叶变换作为系数的多项式\(B(x)\),取单位根的倒数\((1,\omega _{n}^{-1},\omega _{n}^{-2},……,\omega _{n}^{-(n-1)})\)作为\(x\)代入,再将结果除以\(n\)即为\(A(x)\)的系数

这个结论实现了将多项式点值转化为系数

快速傅立叶变换

我们一顿分析最后发现复杂度……仍然是\(O(n^2)\)

我们学这破玩意的意义不就是降低复杂度嘛,所以我们接下来讲怎么降复杂度

我们先设\(A(x)=\sum\limits_{i=0}^{n-1}a_i*x_i\)

我们将\(A(x)\)的下标按奇偶分类,得到
\[A(x)=(a_0+a_2*x^2+……+a_{n-2}*x^{x-2})+(a_1*x+a_3*x^3+……+a_{n-1}*x^{n-1})\]
设两个多项式
\[A_1(x)=a_0+a_2*x+……+a_{n-2}*x^{\frac{n}{2}-1}\]
\[A_2(x)=a_1+a_3*x+……+a_{n-1}*x^{\frac{n}{2}-1}\]
那么就可以发现\(A(x)=A_1(x^2)+xA_2(x^2)\)

将\(x=\omega _n^{k}(k<\frac{n}{2})\)代入
\[A(\omega _n^k)=A_1(\omega _n^{2k})+\omega _n^kA_2(\omega _n^{2k})\]
\[=A_1(\omega _{\frac{n}{2}}^{k})+\omega _n^kA_2(\omega _{\frac{n}{2}}^{k})\]

将将\(x=\omega _n^{k+\frac{n}{2}}(k<\frac{n}{2})\)代入
\[A(\omega _n^{k+\frac{n}{2}})=A_1(\omega _n^{2k+n})+\omega _n^{k+\frac{n}{2}}A_2(\omega _n^{2k+n})\]
\[=A_1(\omega _n^{2k}*\omega _n^n)-\omega _n^kA_2(\omega _n^{2k}*\omega _n^n)\]
\[=A_1(\omega _{\frac{n}{2}}^{k})-\omega _n^kA_2(\omega _{\frac{n}{2}}^{k})\]

我们一点也不惊喜地发现,只要求出\(A_1(x)\)和\(A_2(x)\)在\((\omega _{\frac{n}{2}}^0,\omega _{\frac{n}{2}}^1,……,\omega _{\frac{n}{2}}^{\frac{n}{2}-1})\)的点值表示,就可以\(O(n)\)地求出\(A(x)\)在\((1,\omega _{n},\omega _{n}^{2},……,\omega _{n}^{n-1})\)

所以我们可以递归实现\(O(nlogn)\)求解多项式乘法了

#include<bits/stdc++.h>
using namespace std;
namespace red{
#define eps (1e-8)
    inline int read()
    {
        int x=0;char ch,f=1;
        for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
        if(ch=='-') f=0,ch=getchar();
        while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
        return f?x:-x;
    }
    const int N=5e6+10;
    const double pi=acos(-1.0);
    int n,m;
    int limit=1,len;
    int p[N];
    struct complex
    {
        double x,y;
        complex(double tx=0,double ty=0){x=tx,y=ty;}
        inline complex operator + (const complex t) const
        {
            return complex(x+t.x,y+t.y);    
        }
        inline complex operator - (const complex t) const
        {
            return complex(x-t.x,y-t.y);
        }
        inline complex operator * (const complex t) const
        {
            return complex(x*t.x-y*t.y,x*t.y+y*t.x);
        }
    }a[N],b[N];
    inline void fft(complex *a,int inv)
    {
        for(int i=0;i<limit;++i)
            if(i<p[i]) swap(a[i],a[p[i]]);
        for(int mid=1;mid<limit;mid<<=1)
        {
            complex Wn(cos(pi/mid),inv*sin(pi/mid));
            for(int r=mid<<1,j=0;j<limit;j+=r)
            {
                complex w(1,0);
                for(int k=0;k<mid;++k,w=w*Wn)
                {
                    complex x=a[j+k],y=w*a[j+k+mid];
                    a[j+k]=x+y;
                    a[j+k+mid]=x-y;
                }
            }
        }
    }
    inline void main()
    {
        n=read(),m=read();
        for(int i=0;i<=n;++i) a[i].x=read();
        for(int i=0;i<=m;++i) b[i].x=read();
        while(limit<=n+m) limit<<=1,++len;
        for(int i=0;i<limit;++i)
            p[i]=(p[i>>1]>>1)|((i&1)<<(len-1));
        fft(a,1);
        fft(b,1);
        for(int i=0;i<=limit;++i) a[i]=a[i]*b[i];
        fft(a,-1);
        for(int i=0;i<=n+m;++i) printf("%d ",(int)(a[i].x/limit+0.5));
    }
}
signed main()
{
    red::main();
    return 0;
}

下面先咕咕咕一下

上一篇:快速傅里叶变换(FFT)学习笔记


下一篇:XGBoost 完整推导过程