FFT用途

很简单，多项式相乘。

FFT原理

如果暴力乘，复杂度是$O(n^{2})$的，显然不行

所以考虑点值法。点值表示下的多项式相乘是$O(n)$的。这就是DFT（离散傅里叶变换）。

但很显然，随机带入n个值的复杂度是$O(n^{2})$的。

因此我们考虑以下推导。

设$$F(n)=a_{0}+a_{1}x+a_{2}x^{2}+...+a_{n}x^{n}$$那么我们一定可以将其拆成两个函数，即$$G(n)=a_{0}+a_{2}x+a_{4}x^{2}+a_{6}x^{3}+...$$$$H(n)=a_{1}+a_{3}x+a_{5}x^{2}+a_{7}x^{3}+...$$所以$$F(n)=G(x^{2})+xH(x^{2})$$从这个式子我们可以看出，如果我们取这么一个点值集合，使得其中一半的值为所有值的平方，比如${-1,1}$，那么我们就可以使计算规模减半（对于$G(n)$和$H(n)$，我们只需带入其中一半的值即可）。看来这是个很完美的策略。

很不幸，如果我们只在实数中这么思考，那复杂度仍为$O(n^{2})$，可以说没什么用。但这给了我们一个启发，如果每次对函数进行拆分递归并使点值集合规模减半，我们就可以在$O(nlogn)$的复杂度内求出一个多项式的点值表示。毫无疑问，问题的关键是这个点值集合。那是否存在这样一个点值集合呢？实数域当然没有，但复数域有，也就是单位根！（用复数做自变量真的让我很长时间无法接受）至于单位根么，出门右转问百度

当然，既然用了FFT，自然要把多项式项数补成2的幂，具体做法为高次项补零。

递归式FFT的很好打，但常数大，且容易爆栈。因此我们考虑迭代。

易发现，每次拆分多项式的时候，我们总是按下标的奇偶性来拆分。这就导致了一个结果：原多项式经过递归拆分后的新多项式下标是原多项式下标的二进制反转。

 1 void unit() {
 2     lg[1]=0;
 3     fui(i,2,mxle,1)lg[i]=lg[i/2]+1;
 4     bz[0]=1;
 5     fui(i,1,mxlg,1)bz[i]=bz[i-1]<<1;
 6     n=max(n,m);
 7     n=bz[lg[n]+2]-1;
 8     root[0].x=1;
 9     root[1].x=cos((2*PI)/(n+1));
10     root[1].y=sin((2*PI)/(n+1));
11     fui(i,2,n,1)root[i]=root[i-1]*root[1];                                //求单位根
12     frot[0]=root[0];
13     fui(i,1,n,1)frot[i]=root[n-i+1];                                    //求逆单位根（后面用
14     for(int i=0; i<=n; i++)fz[i]=(fz[i>>1]>>1)|((i&1)<<(lg[n+1]-1)) ;    //二进制反转
15 }

预处理

然后就是FFT了

 1 void pre_FFT() {
 2     fui(i,0,n,1)if(i<fz[i]) {
 3         swp=fft[i];
 4         fft[i]=fft[fz[i]];
 5         fft[fz[i]]=swp;
 6     }
 7 }
 8 void FFT() {
 9     pre_FFT();                                                                        //对原多项式进行二进制反转
10     int zl=(n+1)>>1;                                                                //单位根增量
11     for(int i=1; i<=n; i=(i<<1)|1,zl>>=1) {                                            //每次处理的多项式长度
12         for(int j=0; j<n; j+=i+1) {                                                    //每次处理的多项式左端点
13             int mid=(j+j+i)>>1;
14             for(int k1=j,k2=mid+1,nw=0; k1<=mid; nw+=zl,k1++,k2++) {                //对该多项式的点值进行迭代处理
15                 x1=fft[k2]*root[nw]+fft[k1],x2=fft[k2]*root[nw+((n+1)>>1)]+fft[k1];
16                 fft[k1]=x1;
17                 fft[k2]=x2;
18             }
19         }
20     }
21 }

FFT

记住，千万不要用stl里的复数，会被坑死的，强烈建议手打

IFFT（快速傅里叶逆变换）

FFT求出了点值，然后相乘，求出了结果多项式的点值表示。但这不是结果，我们还要将其“翻译”回系数表示。所以就有了IFFT，可以叫其插值。

插值么，具体做法就是对结果求一遍类似FFT的东西，只不过把单位根改成{${{\omega}^{0},{\omega}^{-1},{\omega}^{-2},{\omega}^{-3},...,{\omega}^{-n}}$}，即上面预处理的时候求过的那个所谓逆单位根带入即可。证明么，我不会，反正作为一个OIer记结论就好了。具体证明参见自为风月马前卒大佬的博客。

 1 void IFFT() {
 2     pre_FFT();
 3     int zl=(n+1)>>1;
 4     for(int i=1; i<=n; i=(i<<1)|1,zl>>=1) {
 5         for(int j=0; j<n; j+=i+1) {
 6             int mid=(j+j+i)>>1;
 7             for(int k1=j,k2=mid+1,nw=0; k1<=mid; nw+=zl,k1++,k2++) {
 8                 x1=fft[k2]*frot[nw]+fft[k1],x2=fft[k2]*frot[nw+((n+1)>>1)]+fft[k1];    //改成逆单位根带入
 9                 fft[k1]=x1;
10                 fft[k2]=x2;
11             }
12         }
13     }
14 }

IFFT

模板题

洛谷P3803 【模板】多项式乘法（FFT）

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 #define INF 0x7fffffff
 4 #define ME 0x7f
 5 #define FO(s) freopen(s".in","r",stdin);freopen(s".out","w",stdout)
 6 #define fui(i,a,b,c) for(int i=(a);i<=(b);i+=(c))
 7 #define fdi(i,a,b,c) for(int i=(a);i>=(b);i-=(c))
 8 #define fel(i,a) for(register int i=h[a];i;i=b[i].ne)
 9 #define ll long long
10 #define ld double
11 #define MEM(a,b) memset(a,b,sizeof(a))
12 #define maxn 1000010
13 #define mxlg 21
14 #define mxle 4194310
15 const ld PI=3.1415926535897932384626;
16 int n,m,mn;
17 ld a0[mxle],b0[mxle];
18 ld ans[mxle];
19 int lg[mxle],bz[mxlg+3];
20 int fz[mxle];
21 struct Complex{
22     ld x,y;
23     Complex(){}
24     Complex(ld __x,ld __y){x=__x;y=__y;}
25     Complex operator *(Complex xx)const{return (Complex){x*xx.x-y*xx.y,x*xx.y+y*xx.x};}
26     Complex operator *(ld xx)const{return (Complex){x*xx,y*xx};}
27     Complex operator +(Complex xx)const{return (Complex){x+xx.x,y+xx.y};}
28 }root[mxle],frot[mxle],a[mxle],b[mxle],fft[mxle],c[mxle],swp,x1,x2;
29 template<class T>
30 inline T read(T &n){
31     n=0;int t=1;double x=10;char ch;
32     for(ch=getchar();!isdigit(ch)&&ch!='-';ch=getchar());(ch=='-')?t=-1:n=ch-'0';
33     for(ch=getchar();isdigit(ch);ch=getchar()) n=n*10+ch-'0';
34     if(ch=='.') for(ch=getchar();isdigit(ch);ch=getchar()) n+=(ch-'0')/x,x*=10;
35     return (n*=t);
36 }
37 template<class T>
38 T write(T n){
39     if(n<0) putchar('-'),n=-n;
40     if(n>=10) write(n/10);putchar(n%10+'0');return n;
41 }
42 template<class T>
43 T writeln(T n){write(n);putchar('\n');return n;}
44 int inc(int &x){int in=(n+1)>>1;for(;x&in;in>>=1)x^=in;x|=in;return x;}
45 void unit(){
46     lg[1]=0;fui(i,2,mxle,1)lg[i]=lg[i/2]+1;
47     bz[0]=1;fui(i,1,mxlg,1)bz[i]=bz[i-1]<<1;
48     n=max(n,m);n=bz[lg[n]+2]-1;
49     root[0].x=1;
50     root[1].x=cos((2*PI)/(n+1));root[1].y=sin((2*PI)/(n+1));
51     fui(i,2,n,1)root[i]=root[i-1]*root[1];
52     frot[0]=root[0];
53     fui(i,1,n,1)frot[i]=root[n-i+1];
54     fui(i,0,n,1)a[i]=root[0]*a0[i];
55     fui(i,0,n,1)b[i]=root[0]*b0[i];
56     for(int i=0;i<=n;i++)fz[i]=(fz[i>>1]>>1)|((i&1)<<(lg[n+1]-1)) ;
57 }
58 void pre_FFT(){
59     fui(i,0,n,1)if(i<fz[i]){
60         swp=fft[i];fft[i]=fft[fz[i]];fft[fz[i]]=swp;
61     }
62 }
63 void FFT(){
64     pre_FFT();int zl=(n+1)>>1;
65     for(int i=1;i<=n;i=(i<<1)|1,zl>>=1){
66         for(int j=0;j<n;j+=i+1){
67             int mid=(j+j+i)>>1;
68             for(int k1=j,k2=mid+1,nw=0;k1<=mid;nw+=zl,k1++,k2++){
69                 x1=fft[k2]*root[nw]+fft[k1],x2=fft[k2]*root[nw+((n+1)>>1)]+fft[k1];
70                 fft[k1]=x1;fft[k2]=x2;
71             }
72         }
73     }
74 //    pre_FFT();
75 }
76 void IFFT(){
77     pre_FFT();int zl=(n+1)>>1;
78     for(int i=1;i<=n;i=(i<<1)|1,zl>>=1){
79         for(int j=0;j<n;j+=i+1){
80             int mid=(j+j+i)>>1;
81             for(int k1=j,k2=mid+1,nw=0;k1<=mid;nw+=zl,k1++,k2++){
82                 x1=fft[k2]*frot[nw]+fft[k1],x2=fft[k2]*frot[nw+((n+1)>>1)]+fft[k1];
83                 fft[k1]=x1;fft[k2]=x2;
84             }
85         }
86     }
87 }
88 int main(){
89     read(n);read(m);mn=m+n;
90     fui(i,0,n,1)read(a0[i]);fui(i,0,m,1)read(b0[i]);unit();
91     fui(i,0,n,1)fft[i]=a[i];FFT();fui(i,0,n,1)a[i]=fft[i];
92     fui(i,0,n,1)fft[i]=b[i];FFT();fui(i,0,n,1)b[i]=fft[i];
93     fui(i,0,n,1)c[i]=a[i]*b[i];
94     fui(i,0,n,1)fft[i]=c[i];IFFT();fui(i,0,n,1)c[i]=fft[i];
95     fui(i,0,mn,1)ans[i]=(c[i].x+0.5)/((n+1));
96     fui(i,0,mn,1)printf("%.0f ",ans[i]);
97     return 0;
98 }

AC代码

常数太大，导致最后两个点2s多。。。

转载于:https://www.cnblogs.com/Axjcy/p/9800440.html

FFT（快速傅里叶变换）摘要相关推荐

【经典算法实现 44】理解二维FFT快速傅里叶变换及 IFFT快速傅里叶逆变换（迭代法和递归法）
[经典算法实现 44]理解二维FFT快速傅里叶变换及 IFFT快速傅里叶逆变换(迭代法和递归法) 一.二维FFTFFTFFT快速傅里叶变换公式推导二.二维FFTFFTFFT 及 IFFTIF ...
FFT快速傅里叶变换超详细的入门学习总结
FFT快速傅里叶变换说明本文创作的目的是为自己巩固该算法,加深印象并深入理解,同时也为FFT入门学者提供一份可鉴的学习总结. 原文链接:https://blog.csdn.net/qq_39565 ...
FFT快速傅里叶变换C语言实现信号处理对振动信号进行实现时域到频域的转换
FFT快速傅里叶变换C语言实现信号处理对振动信号进行实现时域到频域的转换,可实现FFT8192个点或改成其他FFT1024.4096等等,可以直接运行,运行结果与matlab运行的一致,写好了注释, ...
c语言fft乘法步骤,C语言实现FFT(快速傅里叶变换).doc
C语言实现FFT(快速傅里叶变换) 择蚁牙幸帆揣邓淌港烬粹甩滋整维含兔忿茂慨渔下餐随扼哇房坏鹅穆礼围引介害芝共茨恿把喜恤寇杖除冕嗓停揍猫调锚遭傀个碱晓频斌硕宾撕坪莱哩腊养掘蹄轴国繁蔬虞靡砖焙倍勾呸怀怒 ...
快速傅里叶变换c语言函数,C语言实现FFT(快速傅里叶变换)
while(1); } #include #include /********************************************************************* ...
如何 FFT(快速傅里叶变换) 求幅度、频率（超详细含推导过程）
目录如何 FFT(快速傅里叶变换) 求幅度.频率(超详细含推导过程) 一. 打颗栗子二. 求幅度 1. 快速傅里叶变换 2. 求出复数的绝对值 3. 归一化小结三. 求频率 1. 频率公式 ...
FFT快速傅里叶变换的应用——画单边频谱图matlab
FFT快速傅里叶变换的应用--画单边频谱图matlab 快速傅里叶变换在数字信号处理里用的十分广泛,在matlab仿真中,处理信号的时频域变换十分有效,这里结合两个做过的仿真,来说一说fft的应用:画 ...
FFT 快速傅里叶变换初探
一直认为很高深的东西其实也并不很难. 以下内容部分来自qy大神的ppt,同时结合了自己的理解.但理解还不是很深,需要继续研究. 开头首先什么是傅里叶变换:傅立叶变换能将满足一定条件的某个函数表示成三 ...
FFT快速傅里叶变换详解
介绍简而言之,这个东西用来做多项式乘法. 比如说,有两个多项式: 3x2+2x+1,2x2+x+43x^2+2x+1~,~2x^2+x+4 3x2+2x+1 , 2x2+x+4 那么他们乘起来就是 ...
MATLAB——FFT(快速傅里叶变换)
基础知识 FFT即快速傅里叶变换,利用周期性和可约性,减少了DFT的运算量.常见的有按时间抽取的基2算法(DIT-FFT)按频率抽取的基2算法(DIF-FFT). 1.利用自带函数fft进行快速傅里叶 ...

FFT（快速傅里叶变换）摘要

FFT用途

FFT原理

IFFT（快速傅里叶逆变换）

模板题

FFT（快速傅里叶变换）摘要相关推荐

最新文章

热门文章