在学习一项算法以前,咱们总该关心这个算法到底是为了干什么。c++
(如下应用只针对OI)git
一句话:求多项式乘法(固然它的实际用处不少)算法
设多项式函数
\(A(x)=a_0+a_1x+a_2x^2+\ldots+a_nx^n\)学习
\(B(x)=b_0+b_1x+b_2x^2+\ldots+b_mx^m\)优化
咱们想求 \(F(x)=A(x)B(x)=\sum\limits_{i=0}^n\sum\limits_{j=0}^ma_ib_jx^{i+j}\)spa
肉眼可见的复杂度\(O(nm)\).net
而使用通常的FFT则复杂度可以在\(O(nlogn)\) (实则带比较大的常数,但数据范围大时确定比\(O(nm)\)好)code
写这篇博客是为了让本身和看到的人可以梳理一遍FFT这个算法。orm
如下FFT的实现使用比较常见的复数单位根作法。
首先要知道复数是什么,高中选修即有,网上资料也不少,不赘述,这里只讲FFT直接相关。
定义: \(z^n=1\) 那么 \(z\) 就是n次的单位根
对于一个\(n\)次单位根,就是将复数单位圆等分红\(n\)份,每一个数 \(w_n^k (k\in[0,n-1) )\) 就被称做n次的第k个单位根
如图
图中表示了全部8次单位根,其中\(w_8^0=1,w_8^1=1+i,w_8^2=i,w_8^3=-1+i\ldots\) 所谓单位根就是这些复数,它们的平方均为1。
虽然单位根的次数\(n\)在定义上能够取全体正整数,可是咱们在FFT中只须要考虑n为2的整数次幂,以方便咱们利用性质。
单位根有一些性质,咱们先列出来看看 (下边的n均为2的整次幂!)
\(w_n^i\times w_n^j=w_n^{i+j}\),这个性质实际上不能叫性质,从复数乘法,辐角的概念出发即可以得证。
\(w_{2n}^{2k}=w_n^k\)
理解,能够看到 \(2k\) 和 \(k\) 在 \(2n\) 和 \(n\) 中占的比例是相同的,联想单位圆便可(\(w_8^2\)显然等于\(w_4^1\),均为1)
我看了好几篇博客,包括知乎,对折半引理的说法是: 对于 \(n\) 等于 \(2\) 的整次幂, \(w_n^k\)通过平方能够获得全部\(w_{\frac{n}{2}}^k\)(单纯来讲只须要\(n\)是二的倍数就好)
两个小性质均可以证实
\(w_n^{2k}=w_{\frac{n}{2}}^k\)
由消去引理,这个小性质是显然的。从这个性质出发,若是 \(k\) 取遍\([0,n-1]\),那么首先右边是全了,对于左边,\(2k\in[0,2n-2]\),经历了一次循环, (按道理这里须要指出一个循环性质,但我以为不必) \(w_n^{k+n}=w_n^{k}\times w_n^{n}\), 而\(w_n^n\)即\(w_n^0=1\),因此 \(2k\) 过了\(n-1\)以后会通过\(\frac{n}{2}\)个相同的点
\(w_n^k=-w_n^{k+\frac{n}{2}}\)
旋转\(180^\circ\) 事情变得显然;从复数乘法的角度来讲也能够:\(w_n^{k+\frac{n}{2}}=w_n^k\times w_n^{\frac{n}{2}}\),而\(w_n^{\frac{n}{2}}\)刚好为\(-1\)。
那么这个性质告诉咱们前一半\(k\in[0,\frac{n}{2}-1]\)和后一半 \(k\in[\frac{n}{2},n-1]\)的平方是同样的,因此同理,两半的平方值是同样的,也就是两部分平方只肯定了\(\frac{n}{2}\)个值,那么利用前边第一个小性质,咱们不须要循环性质就可以证实折半引理。
折半引理告诉咱们,若是知道了n的全部单位根,那么\(\frac{n}{2}\)的单位根咱们能够知道,反之,咱们利用第二个性质也能够从\(\frac{n}{2}\)的全部单位根,加一个负号,就能推出全部的 \(n\) 次单位根。
\(\sum\limits_{i=0}^{n-1}w_n^k=0\)
从折半引理看,咱们把 \(w_n^k\) 分红两部分\(w_n^a (a\in[0,\frac{n}{2}-1] ),w_n^b(b\in[\frac{n}{2},n-1])\),也就是圆的上半部分和下半部分,对应有每一个\(w_n^a+w_n^b=0\),其实是\(w_n^a+w_n^{a+\frac{n}{2}}=0\),因此两部分的累加和就是0(咱们上边已经假设n是2的整次幂了,因此必定是相同大小的两部分)
首先简单说一下辐角,在复数平面上,从实数轴正半轴逆时针旋转到复数\(z\)到原点的连线,中间通过的角度就叫作辐角(大小在\(2\pi\)以内)
从圆的知识类比一下,咱们能够知道
在叙述 \(n\) 次单位根这件事情的时候,咱们是将复数单位圆 \(n\) 等分以后获得单位根,这样来看的话,显然\(w_n^1\)的辐角是 \(\theta=\frac{2\pi}{n}\), 那么\(w_n^k=(w_n^1)^k\) ,那么\(w_n^k\)的辐角就是 \(k\theta\),咱们知道在普通二维平面上单位圆上的点表示为\((cos\theta,sin\theta)\),类比到复平面,就会获得 \(w_n^k=cos\theta+sin\theta\ i \ \ \ \ (\theta=\frac{2\pi}n)\)
至此,咱们已经了解单位根的性质及其表示,然而咱们还不知道为何要介绍这么一个东西。下面咱们会了解的。
一开始,咱们作介绍时设了两个多项式
\(A(x)=a_0+a_1x+a_2x^2+\ldots+a_nx^n\)
\(B(x)=b_0+b_1x+b_2x^2+\ldots+b_mx^m\)
仔细看的话,咱们发现这两个多项式像什么?没错,像函数,实际上它彻底能够当成一个函数表达式。
即\(y=A(x)\)
咱们知道,求解一次函数只须要知道两个点的坐标,求解二次函数只须要知道三个点的坐标,那么咱们能够换个说法,两个点肯定了一个一次函数,三个点肯定了一个二次函数,一样的,咱们能够说 \(n+1\) 个点肯定了一个 \(n\) 次函数,肯定的方法很简单,咱们带入联当即可(想一想咱们是怎么经过两个点肯定一条直线的)。
如今咱们把名词换一下,”函数“变成“多项式”,就能够改为说 \(n+1\) 个不一样点肯定一个 \(n\) 次多项式。
即点集\(((x_0,y_0),(x_1,y_1),\ldots ,(x_n,y_n))\) 可以肯定一个 \(n\) 次多项式。
代入\(y=A(x)\)看的形象一点 ,咱们解下边的方程组(点是咱们已知且在函数上的),就可以求出来全部的系数 \(a\) (就是和求函数解析式同样的概念),若是用高斯消元来解这样一个方程组,复杂度是 \(O(n^3)\)的,这个过程被称做离散傅里叶逆/反变换(IDFT) 。
同时,若是咱们知道一个多项式的系数表达式,咱们只须要随便找\(n+1\)个点 \(x\) 代入到\(A(x)\)中,就能够获得全部的\(y\)。
这样咱们就能够将一个多项式从系数表达式 \(O(n^2)\) 的转换成点值表达式。
把多项式的系数表示换成点值表示的变换被称做离散傅里叶变换(DFT)。(不当心先介绍了逆变换)
点值表达式之间想要肯定新的点值表达式,只用将y相乘便可(\(x\)对应相等),好比说
\(A=((x_0,y_0),(x_1,y_1)\ldots,(x_n,y_n))\)
\(B=((x_0,y_1'),(x_1,y_2')\ldots,(x_n,y_n') )\)
\(F=AB=((x_0,y_1y_1'),(x_1,y_2y_2'),(x_2,y_2y_2')\ldots,(x_n,y_ny_n') )\)
由于咱们的每一个\(y_iy_i'\)都是一次多项式乘法的结果,按照定义,\(F(x_k)=A(x_k)B(x_k)=y_ky_k'=\sum\limits_{i=0}^n\sum\limits_{j=0}^ma_ib_jx_k^{i+j}\)咱们的多项式乘法自己是只关心系数的,和 \(x\) 无关(就像函数表达式\(y=f(x)\)只和各项系数有关而与具体的点无关同样),而根据咱们点值表示法的定义 \(y_k=F(x_k)\),因此\((x_0,y_1y_1')\)就是\(F(x)\)上的一个点 \((x_k,F(x_k))\)
一个小问题 一开始咱们知道,多项式\(A(x)B(x)\)的结果理应获得最高次项为\(n+m\)的\(F(x)\),而A,B的直接转化的点值表示法分别只有\(n\)个点和\(m\)个点,最后要真正肯定\(F(x)\),须要 \(n+m+1\)个点怎么办?好办,咱们一开始把A和B多算一些y,补若干个点补到\(n+m+1\)个点,多一些点并不会影响什么,而这样就可以创造出\(n+m+1\)个可以肯定\(F(x)\)的点集了。
AB的运算获得的新点集,就可以肯定新的 多项式 \(F(x)\) 。而在已知 \(y_i\)和 \(y_i'\) 的状况下,运算AB须要的复杂度是 \(O(n)\)的!
从前置知识点值表示法的内容中,咱们能够获得一种新的求多项式 \(F(x)=A(x)B(x)\) 的方法
DFT,将\(A(x)\) \(B(x)\) 经过离散傅里叶变换 (DFT)转化为点值表示法,为了可以真正肯定\(F(x)\),A和B都要转化为\(n+m+1\)个点以肯定最高次项为\(n+m\)的\(F(x)\),复杂度 \(O(n^2)\)
将\(A(x)\)的全部点值和\(B(x)\)的全部点值对应相乘,获得\(F(x)\)的点值表示法。复杂度 \(O(n)\)
IDFT,以全部 \(x\) 的各个次项作系数矩阵(以\(x_i^1,x_i^2\ldots\)作系数矩阵),进行高斯消元,获得\(F(x)\)的全部系数。复杂度 \(O(n^3)\)
很是漂亮啊!!咱们一通操做,直接获得一个O(n³)的作法可以进行多项式乘法!!秒!
咱们进行这么多操做,研究了点值表示法可以对应多项式的系数表示法的事实,固然不是为了花里胡哨获得一个\(O(n^3)\)的作法
这个时候,咱们就要考虑优化了,优化什么呢?咱们不忍心点值表示运算明明已经达到了极优秀的\(O(n)\),而DFT和IDFT却严重拖慢运行。
这个时候,咱们就要用到复数单位根了。
平时咱们作函数之类的事情,一般规定 \(x\in R\) 实际上咱们能够扩展数域, 使得 \(x\in C\) 即让\(x\)从实数扩展到整个复数域,做为一个横坐标来使用。而后咱们考虑将DFT的过程当中任选的\(x\),变成有规律的复数单位根,为了使得这些 \(x\) 有规律咱们还规定,对于一个\(n-1\)次多项式,咱们向多项式 \(A(x)\) 和\(B(x)\) 代入全部的单位根,即代入全部的 \(w_n^k\) ,而全部的单位根是不一样的复数,因此点集\((w_n^k,A(w_n^k))\) 就可以做为\(A(x)\)的点值表示法,\(B(x)\)同理。
接下来咱们将 \(n\) 次单位根代入到咱们的多项式中,为了使用咱们上边指出过的复数单位根的性质,以及为了可以惟一肯定咱们的目标多项式,咱们选择的 \(n\) 应当是2的整次幂,而且它要大于等于n+m+1。
代入以后 \(A(x)=((w_n^0,A(w_n^0),),(w_n^1,A(w_n^1))\ldots,(w_n^{n-1},A(w_n^{n-1}))\),考虑如何快速求出每一个\(A(w_n^k)\),咱们从消去引理知道,只要知道偶数项的\(w_n^k=w_\frac{n}2^{\frac{k}2}\),这提示咱们须要将原来的表达式进行一波奇偶分组,从系数表达式入手,如今\(A(x)\)是一个\(n-1\)次多项式。
设\(A_1(x)=(a_0+a_2x+a_4x^2+\ldots+a_{n-2}x^{\frac{n}2-1}),A_2(x)=(a_1+a_3x+\ldots+a_{n-1}x^{\frac{n}2-1})\)
咱们会获得\(A(x)=A_1(x^2)+xA_2(x^2)\)
将全部\(x\) 写成单位根的样子
对于\(k\in[0,\frac{n}2-1]\),\(A(w_n^k)=A_1(w_\frac{n}2^k)+w_n^kA_2(w_\frac{n}2^k)\)
对于\(k\in[\frac{n}2,n-1]\), 从\(w_n^k=-w_n^{k+\frac{n}{2}},w_n^k=w_\frac{n}2^\frac{k}2\)因为\(k\)已经大于\(\frac{n}2\),因此这里全部的\(A_1(w_\frac{n}2^k)=A_1(w_\frac{n}2^{k-\frac{n}2})\),\(A_2\)同理
因此\(A(w_n^k)=A_1(w_{\frac{n}2}^{k-\frac{n}2})-w_n^{k-\frac{n}2}A_2(w_{\frac{n}2}^{k-\frac{n}2})\) 也就是咱们能够直接引用在前一半已经算出来过的\(A_1A_2\),来计算出后一半的全部\(A(w_n^k)\) ,这样咱们把计算的规模缩小了一半,同理,咱们要知晓前一半\(A_1\)和\(A_2\)的值,也能够从一半的前一半推之,即递归即可以求解,最后须要进行log层计算,而后只须要在递归的底层对一个点进行 \(O(n)\) 的多项式运算而后每层\(O(n)\)的向上合并,最终的复杂度是 \(T(n)=2T(\frac{n}2)+O(n)\) 其最终复杂度是\(O(nlogn)\)的。这样咱们经过将多项式拆分红两部分而后进行类似性的递归,使得咱们的DFT可以在\(O(nlogn)\) 的复杂度下获得多项式的点值表达式,这个过程就被称做FFT。
如何快速的把获得的点值表示法的多项式变成系数表示法,先放一个结论吧
也就是说,若是咱们知道了一个多项式,能够用FFT化做点值表示,而若是咱们知道了全部点值反过来求全部系数,咱们发现求每一个系数的形式和求每一个点值的形式是相同的不一样的是多了一个系数\(\frac{1}n\)以及在\(w\)上标多了一个负号,而这是不影响单位根的三项引理的,也就说,这个反向的过程依然能够用FFT解决。如此一来,咱们就能够一样的在\(O(nlogn)\) 的复杂将点值表示转化为系数表示。
关于IFFT的证实:
咱们把每一个点展开成多项式,即把\((w_n^k,A(w_n^k) )\)还原成多项式的模样而后以多项式的系数做为列向量相乘,就可以获得全部的\(y\)
咱们的点值表示法所获得的矩阵是一个范德蒙矩阵,而咱们进行IFFT的过程就是已知向量 \(y\) 以及点矩阵\(W\) 求向量\(a\)
那么咱们只须要计算 \(yW^{-1}\) 便可,对于这个范德蒙矩阵求逆。
直接构造矩阵S
根据求和引理,\(c_{i,j}=w_{i,k}s_{k,j}\)两个矩阵相乘就会获得矩阵C
这个矩阵乘上一个\(\frac{1}n\)就是单位矩阵,因此以\(w_n^{-k}\)做为点的横坐标依次代入矩阵再乘一个\(\frac{1}n\)就是咱们须要的\(W^{-1}\),因此只须要咱们再用\(w_n^{-k}\)作一遍FFT,将结果乘\(\frac{1}n\)便可。
这样一来,咱们用点值表示法就可以\(O(nlogn)\)快速求出系数表示了。
蝴蝶变换是FFT的一个优化,之因此要优化,是由于纯递归求解FFT效率并不高,因此考虑如何迭代求解。
每次奇偶分组的时候,咱们都会挑一个当前所在位置编号\(i\) 为偶数的分到左边,为奇数的分到右边。
偶数的偶数的偶数次项到最后,就会是0,多一个奇数就会在头多1...(感性理解)
总之从变换的规律来看,咱们递归到最后,系数的分组是原来系数的二进制反转
好比咱们有一个7次多项式,有8项系数\(a_0-a_7\)
分组获得的就是
你看最后一组的分配,就是每一个数二进制反转的结果。
换言之,咱们一开始若是将系数二进制反转,而后从下往上推,每次都能推出一个长度更长的多项式,这样子问题就能够一步步向上推出全局的答案。
终于,咱们获得了整个FFT的算法(typroa显示共4862词(不包括下边的代码),码了一天,淦)!下面会在例题里贴上一个模板,而后随着作题会更新一些题目总结吧。
【学习笔记】超简单的快速傅里叶变换(FFT)(含全套证实)
单位根-百度百科
OIwiki(里边还有较详细的二进制反转证实,我懒了不想写了)
百度上的傅里叶变换(主要看严谨概念来的,然鹅看不大懂,而后题目就多了个【OI向】)
[学校的非公开资料]这个就莫得链接
代码以下:
#include<bits/stdc++.h> #define reg register typedef long long ll; using namespace std; inline int qr(){ int x=0,f=0;char ch=0; while(!isd igit(ch)){f|=ch=='-';ch=getchar();} while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();} return f?-x:x; } const int maxn=3e6+199; const double PI=3.14159265358979; int n,m; struct Complex{ double x,y; Complex operator+(const Complex &t)const { return (Complex){x+t.x,y+t.y}; } Complex operator-(const Complex &t)const{ return (Complex){x-t.x,y-t.y}; } Complex operator* (const Complex &t)const { return (Complex){x*t.x-y*t.y,x*t.y+y*t.x}; } }a[maxn],b[maxn]; int rev[maxn],bit,tot;//反转的二进制数 void fft(Complex c[],int inv){//inv表示的是单位根转动的方向 for(reg int i=0;i<tot;i++){//蝴蝶变换! if(i<rev[i]){//小于rev[i]时交换一次就好 swap(c[i],c[rev[i]]); } } for(reg int mid=1;mid<tot;mid<<=1){//从底层开始枚举,即从最小的子问题开始 Complex w1=(Complex){cos(PI/mid),inv*sin(PI/mid)};//求得单位根 for(reg int i=0;i<tot;i+=mid<<1){//枚举每个子问题 Complex wk=(Complex){1,0}; for(reg int j=0;j<mid;j++,wk=wk*w1){//扫左半边,此时大项目A的第i+j项就是原来处理出的左半边和右半边的结合,即以前的第i+j项和i+j+mid的组合,模拟递归 Complex x=c[i+j],y=wk*c[i+j+mid];//A1和A2 c[i+j]=x+y,c[i+j+mid]=x-y;//左半边和右半边的处理 }//c存储点值A(wk)=A1*(wk)+wk*A2(wk) } } } int main(){ // freopen(".in","r",stdin); // freopen(".out","w",stdout); n=qr();m=qr(); for(reg int i=0;i<=n;i++) a[i].x=qr(); for(reg int j=0;j<=m;j++) b[j].x=qr(); while((1<<bit)<n+m+1) bit++;//找到第一个大于n+m+1的2的整次幂 tot=1<<bit;//扩展多项式长度 for(reg int i=0;i<tot;i++){ rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));//记录每一个数的二进制反转 } //rev记录的是每一个二进制数的二进制逆转,咱们这个递推是求出来除去第0位的逆再把第0位放到最高位,就可以获得这个数的逆 fft(a,1),fft(b,1);//转化为点 for(reg int i=0;i<tot;i++) a[i]=a[i]*b[i];//求出目标多项式的点值表示 fft(a,-1);//求逆转化 for(reg int i=0;i<=n+m;i++){ printf("%d ",(int)(a[i].x/tot+0.5));//除去tot,就像咱们以前分析的 ,+0.5来以防精度出问题 } return 0; }