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