快速傅里叶变换（FFT）（学习笔记）

学习了一波FFTFFTFFT，只是浅浅的入门。还有很多前置知识，有一些还不是太了解，完了深入学习之后再补博客qwqqwqqwq
以下内容大部分参考秦岳学长的课件

多项式

形如A(x)=∑k=0n−1ak×xkA(x)=\sum_{k=0}^{n-1}a_k\times x^kA(x)=∑k=0n−1ak×xk，其中aka_kak为多项式系数
C(x)=A(x)+B(x)C(x)=A(x)+B(x)C(x)=A(x)+B(x) ： ck=ak+bkc_k=a_k+b_kck=ak+bk
C(x)/A(x)=B(x)C(x)/A(x)=B(x)C(x)/A(x)=B(x) ： cn=∑ak×bn−kc_n=\sum a_k\times b_{n-k}cn=∑ak×bn−k

多项式的表示法：
系数表达式：上面说的
点值表达式：给出N+1N+1N+1个不同的xxx代入A(x)A(x)A(x)的点值，这样的N+1N+1N+1个元素构成的集合{(x0,y0),(x1,y1),…,(xn,yn)}\{(x0,y0),(x1,y1),…,(xn,yn)\}{(x0,y0),(x1,y1),…,(xn,yn)}，其中yiy_iyi=A(xi)A(x_i)A(xi)，称为点值表示法。

唯一性定理：证明可用范德蒙矩阵行列式

系数与点值：
系数表示法的nnn次多项式A(x),B(x)A(x),B(x)A(x),B(x)可以O(n)O(n)O(n)快速求得A(x)+B(x)A(x)+B(x)A(x)+B(x)，但是卷积通常需要O(n2)O(n^2)O(n2)计算

点值表示法的nnn次多项式卷积也可以O(n)O(n)O(n)计算，只需将值域yyy对应相乘即可{(xi,yi)}×{(xi,zi)}={(xi,yi∗zi)}\{(xi,yi)\}×\{(xi,zi)\}=\{(xi,yi*zi)\}{(xi,yi)}×{(xi,zi)}={(xi,yi∗zi)} (由于nnn项确定次数界为nnn的多项式，故计算卷积是至少保留2n2n2n项)

实际上点值表示法的卷积计算的是圆周卷积

点值与系数的转化：
1.1.1.系数表示法→→→点值表示法
朴素计算O(n3)O(n^3)O(n3)
秦九韶算法/霍纳法则O(n2)O(n^2)O(n2)
2.2.2.点值表示法→→→系数表示法
高斯消元O(n3)O(n^3)O(n3)
拉格朗日插值法O(n2)O(n^2)O(n2)：A(x)=∑0<=k<=nyk∗Π(j≠k)(x−xj)/(xk−xj)A(x)=∑_{0<=k<=n}y_k*Π(j≠k)(x-x_j)/(x_k-x_j)A(x)=∑0<=k<=nyk∗Π(j̸=k)(x−xj)/(xk−xj)
(原理:构造求和式当x=xkx=x_kx=xk时只有一项为yky_kyk，其余全为000)

(霍纳法则不会qwqqwqqwq，等以后补

而FFTFFTFFT可以在O(nlogn)O(nlogn)O(nlogn)时间内求出多项式卷积

复数

单位复数i=sqrt(−1)i=sqrt(-1)i=sqrt(−1)，复数表示为(a+b×i)(a+b\times i)(a+b×i)的形式，运算和实数类似。
其实复数就相当于向量，复数的计算可以用向量计算来解决

struct complex{double x,y;complex(double xx=0,double yy=0) {x=xx,y=yy;}
}a[maxn],b[maxn];
complex operator +(complex a,complex b) {return complex(a.x+b.x,a.y+b.y);}
complex operator -(complex a,complex b) {return complex(a.x-b.x,a.y-b.y);}
complex operator *(complex a,complex b) {return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

证明用向量推一推

欧拉公式：eix=cos(x)+isin(x)e^{ix}=cos(x)+isin(x)eix=cos(x)+isin(x) (泰勒级数)
NNN次单位根：11/n1^{1/n}11/n，在复数域中有nnn个解：e2πik/ne^{2πik/n}e2πik/n
主NNN次单位根：ωn=e2πi/nω_n=e^{2πi/n}ωn=e2πi/n

由于NNN次单位根的循环性质，若将其作为点值的xix_ixi将会简化计算至O(nlogn)O(nlogn)O(nlogn)

欧拉公式的证明可以把公式化一化很容易推
主NNN次单位根可以几何意义来理解，就是复平面上转nnn次可以回到xxx轴

complex Wn(cos(2.0*Pi/lim),1.0*type*sin(2.0*Pi/lim));

快速傅里叶变换

（公式不好打直接贴图）

有点值表达式变回系数表达式是快速傅里叶逆变换，又叫IFFTIFFTIFFT

图解快速傅里叶变换：

代码实现的时候FFTFFTFFT分递归版和非递归版，递归版就是上面介绍的最基础的，而非递归的要用到蝴蝶变换。

递归版：


void FFT(complex *F,int lim,int type){if(lim==1) return;complex a1[(lim>>1)+5],a2[(lim>>1)+5];for(int i=0;i<=lim;i+=2) a1[i>>1]=F[i],a2[i>>1]=F[i+1];//拆成奇数项和偶数项 FFT(a1,lim>>1,type); FFT(a2,lim>>1,type);complex Wn(cos(2.0*Pi/lim),1.0*type*sin(2.0*Pi/lim)),w(1,0);for(int i=0;i<(lim>>1);i++,w=w*Wn){F[i]=a1[i]+w*a2[i];F[i+(lim>>1)]=a1[i]-w*a2[i];}
}

蝴蝶变换

（不知为什么像素如此渣，将就着看qwq）

实现的时候要预处理一个revrevrev数组（说实话只要记住就好原理什么的不重要）

while(limit<=n+m) limit<<=1,l++;for(int i=0;i<limit;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));

非递归版FFTFFTFFT：

void fft(complex *F,int type){for(int i=0;i<limit;i++)if(i<rev[i]) swap(F[i],F[rev[i]]);for(int mid=1;mid<limit;mid<<=1){complex Wn(cos(Pi/mid),type*sin(Pi/mid));//单位根for(int r=mid<<1,j=0;j<limit;j+=r){//r是区间长度，j是位置complex w(1,0);//幂for(int k=0;k<mid;k++,w=w*Wn){complex x=F[j+k],y=w*F[j+mid+k];F[j+k]=x+y; F[j+mid+k]=x-y;}}}
}

因为非递归版比递归版常数小所以现在常用非递归版的

例题

luogu3803模板题

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int maxn=4e6+10;
const double Pi=acos(-1.0);inline int rd(){int x=0,f=1;char c=' ';while(c<'0' || c>'9') f=c=='-'?-1:1,c=getchar();while(c<='9' && c>='0') x=x*10+c-'0',c=getchar();return x*f;
}struct complex{double x,y;complex(double xx=0,double yy=0){x=xx,y=yy;}
}a[maxn],b[maxn];
complex operator +(complex a,complex b) {return complex(a.x+b.x,a.y+b.y);}
complex operator -(complex a,complex b) {return complex(a.x-b.x,a.y-b.y);}
complex operator *(complex a,complex b) {return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}int n,m,l,rev[maxn],limit=1;void fft(complex *F,int type){for(int i=0;i<limit;i++)if(i<rev[i]) swap(F[i],F[rev[i]]);for(int mid=1;mid<limit;mid<<=1){complex Wn(cos(Pi/mid),type*sin(Pi/mid));//单位根for(int r=mid<<1,j=0;j<limit;j+=r){//r是区间长度，j是位置complex w(1,0);//幂for(int k=0;k<mid;k++,w=w*Wn){complex x=F[j+k],y=w*F[j+mid+k];F[j+k]=x+y; F[j+mid+k]=x-y;}}}
}int main(){n=rd(); m=rd();for(int i=0;i<=n;i++) a[i].x=rd();for(int i=0;i<=m;i++) b[i].x=rd();while(limit<=n+m) limit<<=1,l++;for(int i=0;i<limit;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-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));return 0;
}

bzoj2179
其实就是A∗BproblemA*BproblemA∗Bproblem

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define maxn 240005
using namespace std;
int n,lena,lenb,limit=1,l,rev[maxn],ans[maxn];
char s1[maxn],s2[maxn];
const double Pi=acos(-1.0);struct complex{double x,y;complex(double xx=0,double yy=0){x=xx,y=yy;}
}a[maxn],b[maxn];
complex operator +(complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
complex operator -(complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
complex operator *(complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}void FFT(complex *F,int type){for(int i=0;i<limit;i++)if(i<rev[i]) swap(F[i],F[rev[i]]);for(int mid=1;mid<limit;mid<<=1){complex Wn(cos(Pi/mid),1.0*type*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=F[j+k],y=w*F[j+mid+k];F[j+k]=x+y,F[j+mid+k]=x-y;}}}
}int main(){scanf("%d%s%s",&n,s1,s2); lena=lenb=n-1;for(int i=0;i<n;i++) a[n-i-1].x=s1[i]-'0';for(int i=0;i<n;i++) b[n-i-1].x=s2[i]-'0';while(a[lena].x==0) lena--;while(b[lenb].x==0) lenb--;while(limit<=lena+lenb) limit<<=1,++l;for(int i=0;i<limit;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-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<=lena+lenb;i++){ans[i]+=(int)(a[i].x/limit+0.5);if(ans[i]>9) ans[i+1]+=ans[i]/10,ans[i]%=10;}while(ans[lena+lenb+1]) lena++;int now;for(now=lena+lenb;now;now--)if(ans[now]) break;for(;now>=0;now--) printf("%d",ans[now]);return 0;
}

bzoj2194
也是非常裸的模板题了，了解性质就能知道把数组翻转一下就好了

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#define maxn 400005
using namespace std;
const double Pi=acos(-1.0);inline int rd(){int x=0,f=1;char c=' ';while(c<'0' || c>'9') f=c=='-'?-1:1,c=getchar();while(c<='9' && c>='0') x=x*10+c-'0',c=getchar();return x*f;
}struct complex{double x,y;complex(double xx=0,double yy=0){x=xx,y=yy;}
}a[maxn],b[maxn];
complex operator +(complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
complex operator -(complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
complex operator *(complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}int n,limit=1,l,rev[maxn];void FFT(complex *F,int type){for(int i=0;i<limit;i++)if(i<rev[i]) swap(F[i],F[rev[i]]);for(int mid=1;mid<limit;mid<<=1){complex Wn(cos(Pi/mid),1.0*type*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=F[j+k],y=w*F[j+mid+k];F[j+k]=x+y,F[j+mid+k]=x-y;}}}
}int main(){n=rd();for(int i=0;i<n;i++) a[n-i-1].x=rd(),b[i].x=rd();while(limit<2*n) limit<<=1,l++;for(int i=0;i<limit;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-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=n-1;i>=0;i--) printf("%d\n",(int)(a[i].x/limit+0.5));return 0;
}

哦对了多说一句，FFTFFTFFT通常可以用来解决长这个样子的式子：
c(n)=∑i=0n−1a(i)×b(n−i)c(n)=\sum_{i=0}^{n-1}a(i)\times b(n-i)c(n)=∑i=0n−1a(i)×b(n−i)（这其实就是卷积）
也就是下标和为定值

然后是一些不那么裸的题（但其实还是很裸
ZJOI2014力
把分子和分母分别看成两个多项式然后前面后面分开算，只要翻转一个就好了

#include<iostream>
#include<cstdio>
#include<cmath>
#define maxn 400005
#define LL long long
using namespace std;
const double Pi=acos(-1.0);struct complex{double x,y;complex(double xx=0,double yy=0){x=xx,y=yy;}
}a1[maxn],a2[maxn],b[maxn];
complex operator +(complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
complex operator -(complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
complex operator *(complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}int n,rev[maxn],l,limit=1;inline void FFT(complex *F,int type){for(int i=0;i<limit;i++)if(i<rev[i]) swap(F[i],F[rev[i]]);for(int mid=1;mid<limit;mid<<=1){complex Wn(cos(Pi/mid),type*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=F[j+k],y=w*F[j+k+mid];F[j+k]=x+y,F[j+mid+k]=x-y;}}}
}int main(){scanf("%d",&n);for(int i=1;i<=n;i++) scanf("%lf",&a1[i].x),a2[n-i+1].x=a1[i].x;for(int i=1;i<=n;i++) b[i].x=1.0/i/i;while(limit<=2*n) limit<<=1,l++;for(int i=0;i<limit;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));FFT(a1,1); FFT(b,1);for(int i=0;i<=limit;i++) a1[i]=a1[i]*b[i];FFT(a1,-1); FFT(a2,1);for(int i=0;i<=limit;i++) a2[i]=a2[i]*b[i];FFT(a2,-1);for(int i=1;i<=n;i++)printf("%.5lf\n",(a1[i].x-a2[n-i+1].x)/(double)limit);return 0;
}

HNOI2017礼物
bzoj3160fft+manacher
bzoj4503

再说一句这玩意卡精度会死人的。。掉精度掉的很厉害
所以看起来现在出题都用NTTNTTNTT？

快速傅里叶变换（FFT）（学习笔记）相关推荐

快速傅里叶变换 FFT 学习笔记
文章目录 FFT ( 快速傅里叶变换 ) 学习笔记参考文章: 多项式系数表示法点值表示法复数前置芝士向量弧度制定义运算法则单位根快速傅里叶变换: 快速傅里叶逆变换 (IFFT): ...
【文文殿下】快速傅里叶变换(FFT)学习笔记
多项式定义形如\(A(x)=\sum_{i=0}^{n-1} a_i x^i\)的式子称为多项式. 我们把\(n\)称为该多项式的次数界. 显然,一个\(n-1\)次多项式的次数界为\(n\). ...
FFT（快速傅里叶变换）学习笔记
简介 FFTFFTFFT (法法塔)是个什么玩意?他的全名叫快速傅里叶变换(然而貌似和傅里叶并没有太大关系),用来快速求出多项式的点值表示,这个东西一般用来解决多项式相乘的问题. 一般的高精度乘法,我 ...
快速傅里叶变换FFT学习小记
FFT学得还是有点模糊,原理那些基本还是算有所理解了吧,不过自己推这个推不动. 看的资料主要有这两个: http://blog.miskcoo.com/2015/04/polynomial-multi ...
【近万字】分数傅里叶变换课程学习笔记
学习自"课堂在线"平台,北京理工大学陶然教授的课程视频,讲解的非常详细全面,数学公式推导都有,以下为学习笔记,仅记录要点部分. 注:学习此课程,按重要程度排序,需要有信号与系统.数 ...
OpenCV快速傅里叶变换(FFT)用于图像和视讯流的模糊检测
OpenCV快速傅里叶变换(FFT)用于图像和视频流的模糊检测翻译自[OpenCV Fast Fourier Transform (FFT) for blur detection in images ...
图像傅里叶变换（快速傅里叶变换FFT）
学习DIP第7天,图像傅里叶变换转载请标明出处:http://blog.csdn.net/tonyshengtan,欢迎大家转载,发现博客被某些论坛转载后,图像无法正常显示,无法正常表达本人观点,对 ...
基于python的快速傅里叶变换FFT（二）
基于python的快速傅里叶变换FFT(二) 本文在上一篇博客的基础上进一步探究正弦函数及其FFT变换. 知识点 FFT变换,其实就是快速离散傅里叶变换,傅立叶变换是数字信号处理领域一种很重要的算 ...
基于python的快速傅里叶变换FFT（一）
基于python的快速傅里叶变换FFT(一) FFT可以将一个信号变换到频域.有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了.这就是很多信号分析采用FFT变换的原因. ...
MIT 线性代数 Linear Algebra 26:复矩阵，傅里叶矩阵, 快速傅里叶变换 FFT
这一讲我们来讲一下复矩阵.线性代数中,复矩阵是避免不了的话题,因为一个简单实矩阵都有可能有复数特征值. 复矩阵我们着重看一下复矩阵和实矩阵在运算上的区别. 距离首先,一个复数向量的的距离求法发生了 ...

快速傅里叶变换（FFT）（学习笔记）

多项式

复数

快速傅里叶变换

蝴蝶变换

例题

快速傅里叶变换（FFT）（学习笔记）相关推荐

最新文章

热门文章