[2019.2.16]快速傅里叶变换FFT

又是一篇咕咕咕了好久的文章。

什么是多项式

一个形如\(a_0+a_1x_1+a_2x_2+...+a_nx^n\)的东西叫做关于\(x\)\(n\)次多项式。

其中\(a_i\)叫做这个多项式的\(i\)次项系数。特别地,\(a_0\)是这个多项式的常数项。

多项式乘法

一个\(n\)次多项式和一个\(m\)次多项式的乘积是一个\(n+m\)次多项式。

设第一个多项式的系数构成的数列为\(a\),另一个为\(b\)

则这两个多项式的乘积的系数构成的数列\(c\)\(a\)\(b\)的卷积。

即:

\(c_x=\sum_{i=0}^xa_i\times b_{x-i}\)

显然\(c\)\(n+m+1\)项。

如何实现多项式乘法

首先考虑暴力。

直接按照上述式子求出\(c\)即可。时间复杂度\(O(nm)\)

能不能更快呢?

点值表示法

可以证明,对于一个\(n\)次多项式\(f(x)\),我们只需要它在\(n+1\)个不同的\(x\)的取值时对应的\(f(x)\)的之,就可以唯一确定这个多项式。

于是就有了点值表示法。

一个\(n\)次多项式的点值表示法形如一个二元组集合\(P\),\(|P|=n+1\)

\(P=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),...,(x_n,f(x_n))\}\)

有什么用?

如果我们知道\(n\)次多项式\(f(x)\)\(m\)次多项式\(g(x)\)在互不相同的\(x_0,x_1,x_2,...,x_{n+m}\)上的取值,我们就可以知道\(f(x)\times g(x)\)的点值表示

\(P_{f(x)\times g(x)}=\{(x_0,f(x_0)\times g(x_0)),(x_1,f(x_1)\times g(x_1)),(x_2,f(x_2)\times g(x_2)),...,(x_{n+m},f(x_{n+m})\times g(x_{n+m}))\}\)

也就是说,点值表示下的多项式乘法是\(O(n)\)的!

于是我们兴奋地尝试去把这两个多项式转成点值表示。

我们枚举\(1\)\(n+m\)作为选定的\(x\),对于每个\(x\)可以\(O(n+m)\)求得\(f(x),g(x)\)的值。

然后我们\(O(n+m)\)把它乘起来。

于是我们发现这样做的复杂度是\(O((n+m)^2)\)的。

比暴力还要慢QWQ

不过我们的直觉告诉我们,这种奇怪的方法虽然慢,但是有优化的余地。

于是有一个叫傅里叶的神仙,怀着这样的信念,于是就有了——

快速傅里叶变换(FFT)

我们之前的点值表示都局限于实数。

既然实数不行,我们就需要请出......

复数

复数被表示为\(a+bi\),其中\(a,b\)为实数,\(i=\sqrt{-1}\)

我们把复数表示在复平面上,就像我们把实数表示在数轴上。

复平面上有横轴和纵轴,即实轴和虚轴。复平面上的点\((x,y)\)表示复数\(x+yi\)

(个人喜欢用数对的形式表示复数)

然后介绍一下复数的四则运算。

类似实数,把\(i\)当成未知数,注意\(i^2=-1\)即可。

也就是

\((a,b)+(c,d)=(a+b,c+d)\)

\((a,b)-(c,d)=(a-c,b-d)\)

\((a,b)\times(c,d)=ac+adi+bci+bdi^2=(ac-bd,ad+bc)\)

\(\frac{(a,b)}{(c,d)}=\frac{(a,b)\times(c,-d)}{(c,d)\times(c,-d)}=\frac{(ac+bd,bc-ad)}{c^2+d^2}=(\frac{ac+bd}{c^2+d^2},\frac{bc-ad}{c^2+d^2})\)

哦对了,对于复数\((a,b)\),\((a,-b)\)是它的共轭复数。复数除法的过程相当于将分子分母同乘分母的共轭复数。

我们来重点看一看复数乘法。

复数乘法有一个很有用的性质。如果我们将复数\((a,b)\)对应向量\((a,b)\),那么复数的乘法可以表达为:

模长相乘,幅角相加。

模长:对于向量\((a,b)\),其模长为点\((a,b)\)与原点\((0,0)\)的距离;幅角:对于向量\((a,b)\),其幅角为它与横轴正半轴的夹角。

即对于复数乘法\((a,b)\times(c,d)\),其结果所对应的向量的模长为向量\((a,b)\)的模长乘向量\((c,d)\)的模长,幅角为向量\((a,b)\)的幅角加向量\((c,d)\)的幅角。

有啥用?

先别急,我们在复数下定义一个东西。

单位复根

满足\((a,b)^n=1\)的复数\((a,b)\)\(n\)次单位复根。

那么它的分布如何呢?

首先,单位复根的模长必然为1。因为当模长小于1,它乘\(n\)次以后模长还是小于1;大于1同理。

然后,我们发现如果复数\((a,b)\)是一个单位复根,那么它的幅角\(\theta\)一定满足

\(\frac{n\theta}{2\pi}\in Z\)

不难理解吧,也就是一个向量从\((1,0)\)开始,每次旋转\(\theta\)的角度,旋转\(n\)次以后还落在\((1,0)\)上,那么它一定旋转了整数圈。

于是我们发现\(n\)次单位根正好是模长为1,幅角为\(\frac{2k\pi}{n}\)的向量对应的复数。

我们记\(\omega_n^k\)表示模长为1,幅角为\(\frac{2k\pi}{n}\)的向量对应的\(n\)次单位复根,称为第\(k\)\(n\)次单位复根。

我们还能发现,\(\omega_n^k=\omega_n^{k\ mod\ n}\),所以一般情况下,我们认为\(n\)次单位复根有\(n\)个,即\(\omega_n^0,\omega_n^1,\omega_n^2,...,\omega_n^{n-1}\)

然后我们发现它们有一些美妙的性质。

1.\(n\)次单位复根对应的向量\(n\)等分单位圆。

显然吧。因为两个相邻\(n\)次单位复根(认为\(\omega_n^{n-1}\)\(\omega_n^0\)相邻)对应的向量的夹角是相等的。

2.\(\omega_n^k=\omega_n^{k\ mod\ n}\)

这个之前有提到,因为在弧度制下,任意弧度\(\theta\)\(\theta+2k\pi(k\in Z)\)表示一个相同的角。

3.\((\omega_n^k)^p=\omega_n^{kp}\)

\(k\)\(n\)次单位复根对应向量的幅角变为原来的\(p\)倍,自然是对应\(\omega_n^{kp}\)对应的向量啦。

4.\(\omega_n^{2k}=\omega_{\frac{n}{2}}^k\)

前提是\(n\)是偶数。

不难理解,\(n\)次单位复根只是在\(\frac{n}{2}\)次单位复根的基础上增加了\(k\)为奇数的\(\omega_n^{k}\),\(k\)为偶数的部分的集合刚好是\(\frac{n}{2}\)次单位复根。

5.\(\omega_n^{k+\frac{n}{2}}=-\omega_n^k\)

前提是\(n\)是偶数。

相当于一个复数对应的向量进行一次中心对称,复数\((a,b)\)变为\((-a,-b)\)

DFT

终于开始讲FFT了。

不知是不是收到了上天的旨意,傅里叶突发奇想将\(n\)次单位复根带入了\(n-1\)次多项式\(f(x)\),希望得到意想不到的效果。

假定\(n\)是偶数。

\(f(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\)

我们将这\(n\)个单项式按照次数的奇偶性分类。

\(f(x)=(a_0+a_2x^2+a_4x^4+...+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3+a_5x^5+...+a_{n-1}x^{n-1})\)

\(f(x)=(a_0+a_2x^2+a_4x^4+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+a_5x^4+...+a_{n-1}x^{n-2})\)

再将\(\omega_n^k\)带入

\(f(\omega_n^k)=(a_0+a_2(\omega_n^k)^2+a_4(\omega_n^k)^4+...+a_{n-2}(\omega_n^k)^{n-2})+\omega_n^k(a_1+a_3(\omega_n^k)^2+a_5(\omega_n^k)^4+...+a_{n-1}(\omega_n^k)^{n-2})\ \ \ \ (*)\)

\(f(\omega_n^k)=(a_0+a_2\omega_n^{2k}+a_4(\omega_n^{2k})^2+...+a_{n-2}(\omega_n^{2k})^{\frac{n-2}{2}})+\omega_n^k(a_1+a_3\omega_n^{2k}+a_5(\omega_n^{2k})^2+...+a_{n-1}(\omega_n^{2k})^{\frac{n-2}{2}})\)

\(f(\omega_n^k)=(a_0+a_2\omega_\frac{n}{2}^k+a_4(\omega_\frac{n}{2}^k)^2+...+a_{n-2}(\omega_\frac{n}{2}^k)^{\frac{n-2}{2}})+\omega_n^k(a_1+a_3\omega_\frac{n}{2}^k+a_5(\omega_\frac{n}{2}^k)^2+...+a_{n-1}(\omega_\frac{n}{2}^k)^{\frac{n-2}{2}})\)

发现左右两边明显是两个递归子问题。

但是我们发现\(\frac{n}{2}\)次单位复根只有\(\frac{n}{2}\)个,所以上面的式子只适用于\(k<\frac{n}{2}\)的情况。

但是我们需要\(n\)个值来确定这个多项式。

那么怎么求\(\omega_n^{k+\frac{n}{2}}(0\le k<\frac{n}{2})\)呢?

我们知道\(\omega_n^{k+\frac{n}{2}}=-\omega_n^k\)

然后\((*)\)式中括号内每一项的次数都是偶数,所以符号不变,只有括号外的\(\omega_n^k\)变为了\(-\omega_n^k\)

所以

\(f(\omega_n^{k+\frac{n}{2}})=-f(\omega_n^k)=(a_0+a_2\omega_\frac{n}{2}^k+a_4(\omega_\frac{n}{2}^k)^2+...+a_{n-2}(\omega_\frac{n}{2}^k)^{\frac{n-2}{2}})-\omega_n^k(a_1+a_3\omega_\frac{n}{2}^k+a_5(\omega_\frac{n}{2}^k)^2+...+a_{n-1}(\omega_\frac{n}{2}^k)^{\frac{n-2}{2}})\)

然后就是实现。

等一下!

我们之前一直要求\(n\)是偶数。

那万一不是呢?

其实很简单,我们把原来的\(n\)增加到大于等于\(n\)的2的整数次幂即可。

我们定义\(FFT(l,r)\)表示对区间\(a_{l,r}\)进行\(DFT\)。设\(r-l+1=len\),\(ln=\frac{len}{2}\)

\(len=1\),什么都不用做直接返回即可。

首先,奇偶分类。那么我们先将\(a_{l,r}\)复制到\(cpy_{l,r}\),然后从\(l\)\(l+ln-1\)枚举\(i\),计数器\(tmp\)一开始等于\(l\),每次令\(a_i=cpy_{tmp},a_{i+ln}=cpy_{tmp+1}\),并让\(tmp\)增加2。

然后,递归到\(FFT(l,l+ln-1),FFT(l+ln,r)\)

之后定义\(rt\)为第1个\(len\)次单位复根,即\(\omega_n^1\),由它所对应的向量幅角为\(\frac{2\pi}{n}\),模长为1得到\(\omega_n^1=(\cos(\frac{2\pi}{n}),\sin(\frac{2\pi}{n}))\)

定义当前\(n\)次单位复根为\(nw\),一开始为\(\omega_n^0\),即1。

然后将已经经过递归变换的\(a_{l,r}\)复制到\(cpy_{l,r}\),从\(l\)\(l+ln-1\)枚举\(i\),令\(a_i=cpy_i+nw\times cpy_{i+ln},a_{i+ln}=cpy_i-nw\times cpy_{i+ln}\)

结束。

时间复杂度\(O(n\log n)\)

代码会和IDFT一起放出。

IDFT

有了DFT,我们就可以轻松地将点值表示的多项式相乘。

可我们总不能输出点值表示的多项式吧。

我们需要将它转回系数表达的多项式。

其实,IDFT的过程,就是把DFT中的\(\omega_n^1\)换成了它的共轭复数,即\((\cos(\frac{2\pi}{n}),-\sin(\frac{2\pi}{n}))\),得到的系数再除以\(\frac{1}{n}​\)即可。

证明比较长,略。

其实我自己也是一知半解

那么终于到了放代码的时刻。

code:

const double pi=acos(-1);//pi
struct com{//复数结构体
    double re,in;
}f[300010],f1[300010],f2[300010],cpy[300010];
com operator+(const com &x,const com &y){
    return (com){x.re+y.re,x.in+y.in};
}
com operator-(const com &x,const com &y){
    return (com){x.re-y.re,x.in-y.in};
}
com operator*(const com &x,const com &y){
    return (com){x.re*y.re-x.in*y.in,x.re*y.in+x.in*y.re};
}
void FFT(com *f,int l,int len,int op){
    if(len<2)return;
    int nl=len>>1,tmp=l-1;
    for(int i=l;i<l+len;++i)cpy[i]=f[i];//复制原f数组到cpy
    for(int i=l;i<l+nl;++i)f[i]=cpy[++tmp],f[i+nl]=cpy[++tmp];//按照奇偶分类
    FFT(f,l,nl,op),FFT(f,l+nl,nl,op);//递归处理
    com w=(com){cos(pi/(1.0*nl)),sin(pi/(1.0*nl)*op)},nw=(com){1,0},t;//w是第一个单位复根,nw是当前单位复根
    for(int i=l;i<l+nl;++i,nw=nw*w)t=nw*f[i+nl],cpy[i]=f[i]+t,cpy[i+nl]=f[i]-t;
    for(int i=l;i<l+len;++i)f[i]=cpy[i];
}

例题

1.洛谷模板

code:

#include<bits/stdc++.h>
using namespace std;
const double pi=2*acos(-1);
struct comp{
    double re,in;
}tmp[2100010],rt,pw,t,a[2100010],b[2100010];
comp operator+(const comp &x,const comp &y){
    return (comp){x.re+y.re,x.in+y.in};
}
comp operator-(const comp &x,const comp &y){
    return (comp){x.re-y.re,x.in-y.in};
}
comp operator*(const comp &x,const comp &y){
    return (comp){x.re*y.re-x.in*y.in,x.re*y.in+x.in*y.re};
}
int n,m,w,tg,mid,nw,sz;
void Work(comp *p,const int &l,const int &r){
    for(int i=l;i<=r;i++)tmp[i]=p[i];
    nw=l,sz=l-1;
    while(nw<=r)p[++sz]=tmp[nw],nw+=2;
    nw=l+1;
    while(nw<=r)p[++sz]=tmp[nw],nw+=2;
}
void FFT(comp *p,const int &l,const int &nl,const int &op){
    if(nl<2)return;
    Work(p,l,l+nl-1);
    FFT(p,l,nl>>1,op),FFT(p,l+(nl>>1),nl>>1,op);
    rt=(comp){cos(pi/nl),sin(pi/nl)*op},pw=(comp){1,0};
    mid=(l+l+nl-1)/2;
    for(int i=l;i<l+nl;i++)tmp[i]=p[i];
    for(int i=l;i<=mid;i++,pw=pw*rt)t=pw*tmp[mid+i-l+1],p[i]=tmp[i]+t,p[mid+i-l+1]=tmp[i]-t;
}
void scan(double &x){
    x=0;
    char c=getchar();
    while('0'>c||c>'9')c=getchar();
    while('0'<=c&&c<='9')x=x*10+c-'0',c=getchar();
}
void print(int x){
    x>9?print(x/10),0:0,putchar(x%10+'0');
}
int main(){
    scanf("%d%d",&n,&m);
    n++,m++;
    for(int i=0;i<n;i++)scan(a[i].re);
    for(int i=0;i<m;i++)scan(b[i].re);
    n+=m-1,m=1;
    while(m<n)m<<=1;
    FFT(a,0,m,1),FFT(b,0,m,1);
    for(int i=0;i<m;i++)a[i]=a[i]*b[i];
    FFT(a,0,m,-1);
    for(int i=0;i<n;i++)print((int)(fabs(a[i].re)/m+0.1)),putchar(' ');
    return 0;
}

2.高精度乘法洛谷模板

code:

#include<bits/stdc++.h>
using namespace std;
const double pi=2*acos(-1);
struct comp{
    double re,in;
}tmp[150010],rt,pw,t,a[150010],b[150010];
comp operator+(const comp &x,const comp &y){
    return (comp){x.re+y.re,x.in+y.in};
}
comp operator-(const comp &x,const comp &y){
    return (comp){x.re-y.re,x.in-y.in};
}
comp operator*(const comp &x,const comp &y){
    return (comp){x.re*y.re-x.in*y.in,x.re*y.in+x.in*y.re};
}
int n,m,w,tg,mid,nw,sz;
void Work(comp *p,const int &l,const int &r){
    for(int i=l;i<=r;i++)tmp[i]=p[i];
    nw=l,sz=l-1;
    while(nw<=r)p[++sz]=tmp[nw],nw+=2;
    nw=l+1;
    while(nw<=r)p[++sz]=tmp[nw],nw+=2;
}
void FFT(comp *p,const int &l,const int &nl,const int &op){
    if(nl<2)return;
    Work(p,l,l+nl-1);
    FFT(p,l,nl>>1,op),FFT(p,l+(nl>>1),nl>>1,op);
    rt=(comp){cos(pi/nl),sin(pi/nl)*op},pw=(comp){1,0};
    mid=(l+l+nl-1)/2;
    for(int i=l;i<l+nl;i++)tmp[i]=p[i];
    for(int i=l;i<=mid;i++,pw=pw*rt)t=pw*tmp[mid+i-l+1],p[i]=tmp[i]+t,p[mid+i-l+1]=tmp[i]-t;
}
void scan(double &x){
    x=0;
    char c=getchar();
    while('0'>c||c>'9')c=getchar();
    x=c-'0';
}
int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++)scan(a[n-i].re);
    for(int i=1;i<=n;i++)scan(b[n-i].re);
    n=n*2-1,m=1;
    while(m<n)m<<=1;
    FFT(a,0,m,1),FFT(b,0,m,1);
    for(int i=0;i<m;i++)a[i]=a[i]*b[i];
    FFT(a,0,m,-1);
    for(int i=0;i<n;i++)a[i].re=fabs(a[i].re)/m,a[i].re=(int)(w+a[i].re+0.1),w=(int)a[i].re/10,a[i].re-=10*w,(w&&i==n-1)?n++:0;
    for(int i=n-1;i>=0;i--)tg|=a[i].re>0,tg?putchar((int)a[i].re+'0'):0;
    return 0;
}

3.[Hnoi2017]礼物 (洛谷,BZOJ)

题解

posted @ 2019-03-17 18:26 xryjr233 阅读(...) 评论(...) 编辑 收藏