文章目录

前言
FFT是什么？
多项式表示法
- 系数表示法
- 点值表示法
初步想法
问题
解决
问题2
- 复数
- - 定义
  - 运算
  - 复平面
  - 单位圆
  - 单位根的性质
- 解决
问题3
- IFFT（快速傅里叶逆变换）
- code
问题4
- 解决
- Code
后记

前言

最近学会了FFT，一个极其美妙的算法！
众所周知FFT英文全称是 Fast Fast TLE Fast Fourier Transform

FFT是什么？

FFT（快速傅里叶变换）是DFT（离散傅里叶变换）的优化，用于加速多项式乘法。
即：
设 F(x)=∑i=0nfixiF(x)=\sum_{i=0}^nf_ix^iF(x)=∑i=0nfixi , G(x)=∑i=0ngixiG(x)=\sum_{i=0}^ng_ix^iG(x)=∑i=0ngixi
求 F(x)G(x)F(x)G(x)F(x)G(x)

多项式表示法

系数表示法

系数表示法是大家最常见的一种表示多项式的方法。
比如 F(x)=∑i=0nfixiF(x)=\sum_{i=0}^nf_ix^iF(x)=∑i=0nfixi
那么这个多项式我们就可以用 {f0,f1,...,fn}\{f_0,f_1,...,f_n\}{f0,f1,...,fn} 这些系数表示。
普通的暴力多项式乘法也是基于这种系数表示法上的。

点值表示法

先抛出个定理：n+1n+1n+1 个不同的点能够确定一个次数为nnn的多项式
我也不会证
我们把多项式放到平面直角坐标系中就是个函数了。
那这个多项式 F(x)F(x)F(x) 就可以表示为 {(x0,f(x0)),(x1,f(x1)),...,(xn,f(xn))}\{(x_0,f(x_0)),(x_1,f(x_1)),...,(x_n,f(x_n))\}{(x0,f(x0)),(x1,f(x1)),...,(xn,f(xn))}

对于两个函数 F(x)={(x0,f(x0)),(x1,f(x1)),...,(xn,f(xn))}F(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),...,(x_n,f(x_n))\}F(x)={(x0,f(x0)),(x1,f(x1)),...,(xn,f(xn))} 和 G(x)={(x0,g(x0)),(x1,g(x1)),...,(xn,g(xn))}G(x)=\{(x_0,g(x_0)),(x_1,g(x_1)),...,(x_n,g(x_n))\}G(x)={(x0,g(x0)),(x1,g(x1)),...,(xn,g(xn))} 的乘积等于 {(x0,f(x0)g(x0)),(x1,f(x1)g(x1)),...,(xn,f(xn)g(xn))}\{(x_0,f(x_0)g(x_0)),(x_1,f(x_1)g(x_1)),...,(x_n,f(x_n)g(x_n))\}{(x0,f(x0)g(x0)),(x1,f(x1)g(x1)),...,(xn,f(xn)g(xn))}
那这样的时间复杂度就是 O(n)O(n)O(n) 的了！

初步想法

有了点值表示法，我们就有了一点初步的想法：
1.1.1. 读入两个系数表示法的多项式
2.2.2. 把读入进来的两个多项式转换为点值表示法
3.3.3. 用点值表示法做多项式乘法
4.4.4. 把两个多项式转换回系数表示法
5.5.5. 输出

问题

看上去十分美好，多项式乘法 O(n)O(n)O(n) ，不挺棒的吗？
然而，你仔细思考一下会发现系数表示法与点值表示法互相转换是 O(n2)O(n^2)O(n2) 的！
你想啊，我们要选出 n+1n+1n+1 个点，然后每个点代入求函数值，因此是 O(n2)O(n^2)O(n2) 的。

解决

我们观察一下偶函数或者奇函数，我们是不是可以找一对点使得它们是相反数，那么我们就只用求其中一个就好了！
那对于一般的多项式呢？
例如 F(x)=3x5+2x4+x3+7x2+5x+1F(x)=3x^5+2x^4+x^3+7x^2+5x+1F(x)=3x5+2x4+x3+7x2+5x+1 ,我们把它转换一下，按照奇偶性分类得
F(x)=(2x4+7x2+1)+(3x5+x3+5x)F(x)=(2x^4+7x^2+1)+(3x^5+x^3+5x)F(x)=(2x4+7x2+1)+(3x5+x3+5x) ,对于后半段提取一个 xxx 得
F(x)=(2x4+7x2+1)+x(3x4+x2+5)F(x)=(2x^4+7x^2+1)+x(3x^4+x^2+5)F(x)=(2x4+7x2+1)+x(3x4+x2+5) ,那么两部分都是偶函数了！
我们再看，设前面的多项式为 Pe(x2)P_e(x^2)Pe(x2) 后面的多项式为 Po(x2)P_o(x^2)Po(x2)。（为什么是x2x^2x2呢？观察到 Pe(x)=(2x2+7x+1)P_e(x)=(2x^2+7x+1)Pe(x)=(2x2+7x+1)，每一项次数刚好是两倍的关系，那么我们就需要 x2x^2x2）
我们再考虑求 Pe(x2)P_e(x^2)Pe(x2) 和 Po(x2)P_o(x^2)Po(x2)，显然，我们可以递归处理！
按照这样做的话，那FFTFFTFFT的时间复杂度就是 O(nlog2n)O(nlog_2n)O(nlog2n) 的了！

问题2

你以为就这么结束了吗？那有什么美妙的？
你是否观察到一个小问题：相反数平方后剩下的 n2\frac{n}{2}2n 个点都是非负数！这就不再存在相反数了！这个递归也就做不下去了！
那怎么办呢……这就是FFT最美妙的一个地方了，我们把这些点都取成复数！

复数

定义

在介绍取什么样的复数之前，我先讲讲什么是复数。

我们把形如z=a+bi（a,b均为实数）的数称为复数，其中a称为实部，b称为虚部，i称为虚数单位。当z的虚部等于零时，常称z为实数；当z的虚部不等于零时，实部等于零时，常称z为纯虚数。复数域是实数域的代数闭包，即任何复系数多项式在复数域中总有根。复数是由意大利米兰学者卡当在十六世纪首次引入，经过达朗贝尔、棣莫弗、欧拉、高斯等人的工作，此概念逐渐为数学家所接受。
——百度百科

iii为虚数单位，iii 的定义为 i2+1=0i^2+1=0i2+1=0 ，即 iii 是关于 xxx 的方程 x2=−1x^2=-1x2=−1 的一个解。

运算

复数的运算就有点类似多项式了，不过 i2i^2i2 需要约成 −1-1−1
假设有两个复数 z0=a+bi,z1=c+diz_0=a+bi,z_1=c+diz0=a+bi,z1=c+di
则：
z0+z1=(a+c)+(b+d)iz0×z1=(a+bi)(c+di)=ac+adi+bci+bdi2=(ac−bd)+(ad+bc)iz0z1=a+bic+di=(a+bi)(c−di)(c+di)(c−di)=(ac+bd)+(bc−ad)ic2+d2=ac+bdc2+d2+bc−adc2+d2iz_0+z_1=(a+c)+(b+d)i \\ z_0\times z_1=(a+bi)(c+di)=ac+adi+bci+bdi^2=(ac-bd)+(ad+bc)i \\ \frac{z_0}{z_1}=\frac{a+bi}{c+di}=\frac{(a+bi)(c-di)}{(c+di)(c-di)}=\frac{(ac+bd)+(bc-ad)i}{c^2+d^2}=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i z0+z1=(a+c)+(b+d)iz0×z1=(a+bi)(c+di)=ac+adi+bci+bdi2=(ac−bd)+(ad+bc)iz1z0=c+dia+bi=(c+di)(c−di)(a+bi)(c−di)=c2+d2(ac+bd)+(bc−ad)i=c2+d2ac+bd+c2+d2bc−adi

复平面

像实数与数轴上的点一一对应一样，复数能与“平面直角坐标系”中的点一一对应，而这个“平面直角坐标系”就是复平面。
在复平面中，x轴又叫做实轴，在实轴上的点都是实数；同理y轴又叫虚轴，在y轴上的点（原点除外）都是纯虚数。
比如一个复数 z0=a+biz_0=a+biz0=a+bi 则 z0z_0z0 在复平面中表示为点 (a,b)(a,b)(a,b)
特别的，在复平面中，复数也能表示一个向量。

单位圆

知道了什么是复数之后，我们再来考虑如何选数。
这就有点难了。
先讲讲什么是单位圆吧

如图。
然后我们把这个圆平均分成 nnn 份，第 kkk 个点记作 wnkw_n^kwnk，点 wn0w_n^0wn0 就是 (1,0)(1,0)(1,0) 然后接下来就是逆时针的剩下n-1个点。
现在我们思考如何求得 wnkw_n^kwnk，第kkk个点就是 (cosθ,sinθ)(cos\theta,sin\theta)(cosθ,sinθ)，即 wnk=cosθ+sinθiw_n^k=cos\theta+sin\theta iwnk=cosθ+sinθi

结合上图食用更加。
那这个 θ=2πkn\theta=\frac{2\pi k}{n}θ=n2πk

单位根的性质

以上的wnkw_n^kwnk就称为单位根，我们看看它有哪些有用的性质？
1.1.1. w2n2k=wnkw_{2n}^{2k}=w_n^kw2n2k=wnk
证明：wnk=cos2πkn+sin2πkni=cos4πk2n+sin4πk2ni=w2n2kw_n^k=cos\frac{2\pi k}{n}+sin\frac{2\pi k}{n}i=cos\frac{4\pi k}{2n}+sin\frac{4\pi k}{2n}i=w_{2n}^{2k}wnk=cosn2πk+sinn2πki=cos2n4πk+sin2n4πki=w2n2k
2.2.2. wnk=−wnk+n2w_n^k=-w_n^{k+\frac{n}{2}}wnk=−wnk+2n
证明：读者可自己尝试证明，也不难证。
3.3.3. wn0=wnn=0w_n^0=w_n^n=0wn0=wnn=0
证明：同上

解决

有了单位根，我们利用单位根的性质，然后把单位根塞到函数里面，进行分治，时间复杂度就有保证了。

问题3

我们把系数表示法转化成了点值表示法，但我们要输出，所以我们考虑如何转化回系数表示法

IFFT（快速傅里叶逆变换）

我们考虑用矩阵来看。

这是FFT，那我们现在相当于知道左边两个矩阵，要求最右边的那个矩阵——这不就是求矩阵的逆吗（虽然我不会）？
那这个矩阵的逆就是：

其实就是虚部取反，然后再乘上个 1n\frac{1}{n}n1

code

我们发现IFFT和FFT真的长得太像了，我们其实就可以考虑把他们两个的代码合在一起。

#include<cmath>
#include<math.h>
#include<cstdio>
#include<complex>
#include<iostream>
#define PI M_PI
using namespace std;
int n,m1,m2;
complex<double>A[4200005],B[4200005];
void FFT(complex<double> *A, int len,int type) {//type=1表示FFT，type=-1表示IFFT，具体见下if (len==1) return;int mid=len>>1;complex<double> A0[mid],A1[mid];for (int i=0;i<mid;++i){A0[i]=A[i<<1];A1[i]=A[i<<1|1];}FFT(A0,mid,type);FFT(A1,mid,type);complex<double>W(cos(2.0*PI/len),type*sin(2.0*PI/len));complex<double>w(1.0, 0.0);for (int i=0;i<mid;++i) {A[i]=A0[i]+w*A1[i];A[i+mid]=A0[i]-w*A1[i];w*=W;}return;
}
int main(){scanf("%d%d",&m1,&m2);for (int i=0;i<=m1;++i) cin>>A[i];for (int i=0;i<=m2;++i) cin>>B[i];n=m1+m2;//把n凑成2的幂for (int i=1;i<=25;++i) if ((1<<i)>n){n=1<<i;break;}FFT(A,n,1);FFT(B,n,1);for (int i=0;i<=n;++i)A[i]*=B[i];FFT(A,n,-1);for (int i=0;i<=m1+m2;++i)printf("%.0lf ",A[i].real()/n+0.00001);return 0;
}

问题4

真多问题啊
FFT常数真的极大，普通FFT连luogu上板题都过不了。
FFT过不了FFT板题

解决

发现FFT主要是在递归的地方耗费太多的时间了，我们考虑把递归改成迭代时间可能会更优秀一些。
我们看看递归是怎么递归的：
（网上偷了一张图）

发现每个数后来所处的位置其实就是二进制下反过来了而已。
那我们考虑怎样把二进制反过来呢？
很容易可以想到dp，具体看代码，这个dp感觉比较妙。
那现在我们就只需要枚举长度以及起点，就可以做FFT/IFFT了！

Code

终于没有问题了！

#include<cmath>
#include<math.h>
#include<cstdio>
#include<complex>
#include<iostream>
#include<algorithm>
#define PI M_PI
#define N 4000005
using namespace std;
int n,m1,m2;
int rev[N];
complex<double>A[N],B[N];//C++自带函数
void FFT(complex<double> *A, int type) {for (int i=0;i<n;++i)if (i<rev[i]) swap(A[i],A[rev[i]]);for (int j=1;j<n;j<<=1){complex<double>W(cos(PI/j),sin(PI/j)*type);for (int k=0;k<n;k+=(j<<1)){complex<double>w(1.0,0.0);for (int i=0;i<j;++i){complex<double>ye,yo;ye=A[i+k];yo=A[i+j+k]*w;A[i+k]=ye+yo;A[i+j+k]=ye-yo;w*=W;}}}return;
}
int main(){scanf("%d%d",&m1,&m2);for (int i=0;i<=m1;++i) cin>>A[i];for (int i=0;i<=m2;++i) cin>>B[i];n=m1+m2;for (int i=1;i<=25;++i) if ((1<<i)>n){n=1<<i;break;}//dpfor (int i=0;i<n;++i)rev[i]=(rev[i>>1]>>1)|((i&1)?(n>>1):0);FFT(A,1);FFT(B,1);for (int i=0;i<=n;++i)A[i]*=B[i];FFT(A,-1);for (int i=0;i<=m1+m2;++i)printf("%d ",(int)(A[i].real()/(double)n+0.5));return 0;
}

后记

曾经听hqx巨佬说过一件事，就是某次比赛，出题人出了一道多项式乘法的题，然后正解是 n2n^2n2 的，原因是FFT常数太大，暴力打打优化还能比FFT快（笑）。

【学习笔记】极其美妙的算法——FFT（快速傅里叶变换）相关推荐

【学习笔记】多项式相关算法
[学习笔记]多项式相关算法手动博客搬家: 本文发表于20181125 13:19:28, 原地址https://blog.csdn.net/suncongbo/article/details/844 ...
FFT快速傅里叶变换超详细的入门学习总结
FFT快速傅里叶变换说明本文创作的目的是为自己巩固该算法,加深印象并深入理解,同时也为FFT入门学者提供一份可鉴的学习总结. 原文链接:https://blog.csdn.net/qq_39565 ...
java fft 频谱算法_快速傅里叶变换（FFT）算法原理及代码解析
FFT与DFT关系: 快速傅里叶变换(Fast Fourier Transform)是离散傅里叶(DFT)变换的一种快速算法,简称FFT,通过FFT可以将一个信号从时域变换到频域:FFT(快速傅里叶变 ...
【经典算法实现 44】理解二维FFT快速傅里叶变换及 IFFT快速傅里叶逆变换（迭代法和递归法）
[经典算法实现 44]理解二维FFT快速傅里叶变换及 IFFT快速傅里叶逆变换(迭代法和递归法) 一.二维FFTFFTFFT快速傅里叶变换公式推导二.二维FFTFFTFFT 及 IFFTIF ...
Apollo星火计划学习笔记——Apollo开放空间规划算法原理与实践
文章目录前言 1. 开放空间规划算法总体介绍 1.1 Task: OPEN_SPACE_ROI_DECIDER 1.2 Task: OPEN_SPACE_TRAJECTORY_PROVIDER 1. ...
Boost库学习笔记（二）算法模块-C++11标准
Boost库学习笔记(二)算法模块-C++11标准一.综述 Boost.Algorithm是一系列人通用推荐算法的集合,虽然有用的通用算法很多,但是为了保证质量和体积,并不会将太多通用算法通过审查测 ...
【学习笔记】目标检测算法总结
[学习笔记]目标检测算法总结说明 MacOS操作系统. MindNote思维导图软件. B站学习视频+原论文学习. 初学者笔记如有问题请多多指教记录 Overfeat模型.R-CNN.Fast ...
C++ Primer 学习笔记第十章泛型算法
C++ Primer 学习笔记第十章泛型算法 336 find函数 #include <iostream> #include <vector> #include <s ...
FFT快速傅里叶变换C语言实现信号处理对振动信号进行实现时域到频域的转换
FFT快速傅里叶变换C语言实现信号处理对振动信号进行实现时域到频域的转换,可实现FFT8192个点或改成其他FFT1024.4096等等,可以直接运行,运行结果与matlab运行的一致,写好了注释, ...
快速傅里叶变换c语言函数,C语言实现FFT(快速傅里叶变换)
while(1); } #include #include /********************************************************************* ...

【学习笔记】极其美妙的算法——FFT（快速傅里叶变换）