快速傅里叶变换(FFT)详解
本文只讨论FFT在信息学奥赛中的应用
文中内容均为个人理解,如有错误请指出,不胜感激
前言
先解释几个比较容易混淆的缩写吧
DFT:离散傅里叶变换—>O(n^2)计算多项式乘法
FFT:快速傅里叶变换—>O(n*log(n)计算多项式乘法
FNTT/NTT:快速傅里叶变换的优化版—>优化常数及误差
FWT:快速沃尔什变换—>利用类似FFT的东西解决一类卷积问题
MTT:毛爷爷的FFT—>非常nb
多项式
系数表示法
设A(x)表示一个n-1次多项式
则A(x)=sum_{i=0}^{n} a_i * x^i
例如:A(3)=2+3*x+x^2
利用这种方法计算多项式乘法复杂度为O(n^2)
(第一个多项式中每个系数都需要与第二个多项式的每个系数相乘)
点值表示法
将n互不相同的x带入多项式,会得到n个不同的取值y
则该多项式被这n个点(x_1,y_1),(x_2,y_2),dots,(x_n,y_n)唯一确定
其中y_i=sum_{j=0}^{n-1} a_j*x_i^j
例如:上面的例子用点值表示法可以为(0,2),(1,5),(2,12)
利用这种方法计算多项式乘法的时间复杂度仍然为O(n^2)
(选点O(n),每次计算O(n))
我们可以看到,两种方法的时间复杂度都为O(n^2),我们考虑对其进行优化
对于第一种方法,由于每个点的系数都是固定的,想要优化比较困难
对于第二种方法,貌似也没有什么好的优化方法,不过当你看完下面的知识,或许就不这么想了
复数
在介绍复数之前,首先介绍一些可能会用到的东西
向量
同时具有大小和方向的量
在几何中通常用带有箭头的线段表示
圆的弧度制
等于半径长的圆弧所对的圆心角叫做1弧度的角,用符号rad表示,读作弧度。用弧度作单位来度量角的制度叫做弧度制
公式:
1^{circ }=dfrac{pi}{180}rad
180^{circ }=pi rad
平行四边形定则
(好像画的不是很标准。。)
平行四边形定则:AB+AD=AC
复数
定义
设a,b为实数,i^2=-1,形如a+bi的数叫负数,其中i被称为虚数单位,复数域是目前已知最大的域
在复平面中,x代表实数,y轴(除原点外的点)代表虚数,从原点(0,0)到(a,b)的向量表示复数a+bi
模长:从原点(0,0)到点(a,b)的距离,即sqrt{a^2+b^2}
幅角:假设以逆时针为正方向,从$x$轴正半轴到已知向量的转角的有向角叫做幅角
运算法则
加法:
因为在复平面中,复数可以被表示为向量,因此复数的加法与向量的加法相同,都满足平行四边形定则(就是上面那个)
乘法:
几何定义:复数相乘,模长相乘,幅角相加
代数定义:(a+bi)*(c+di)
=ac+adi+bci+bdi^2
=ac+adi+bci-bd
=(ac-bd)+(bc+ad)i
单位根
下文中,默认n为2的正整数次幂
在复平面上,以原点为圆心,1为半径作圆,所得的圆叫单位圆。以圆点为起点,圆的n等分点为终点,做n个向量,设幅角为正且最小的向量对应的复数为omega_n,称为n次单位根。
根据复数乘法的运算法则,其余n-1个复数为omega_n^2,omega_n^3,ldots,omega_n^n
注意omega_n^0=omega_n^n=1(对应复平面上以x轴为正方向的向量)
那么如何计算它们的值呢?这个问题可以由欧拉公式解决omega_{n}^{k}=cos k *frac{2pi}{n}+isin k*frac{2pi}{n}
例如
图中向量AB表示的复数为4次单位根
单位根的幅角为周角的dfrac{1}{n}
在代数中,若z^n=1,我们把z称为n次单位根
单位根的性质
- omega _{n}^{k}=cos kdfrac{2pi}{n}+isin kdfrac {2pi }{n}(即上面的公式)
ω2k2n=ωkn
证明:
omega _{2n}^{2k}=cos 2k*frac{2pi}{2n}+isin2k*frac{2pi}{2n}
=omega _{n}^{k}
- omega _{n}^{k+frac{n}{2}}=-omega _{n}^{k}
omega _{n}^{frac{n}{2}}=cosfrac{n}{2}*frac{2pi}{n}+isinfrac{n}{2}*frac{2pi}{n}
=cos pi+isinpi
=-1
- omega _{n}^{0}=omega _{n}^{n}=1
讲了这么多,貌似跟我们的正题没啥关系啊。。
OK!各位坐稳了,前方高能!
快速傅里叶变换
我们前面提到过,一个n次多项式可以被n个点唯一确定。
那么我们可以把单位根的0到n-1次幂带入,这样也可以把这个多项式确定出来。但是这样仍然是O(n^2)的呀!
我们设多项式A(x)的系数为(a_o,a_1,a_2,ldots,a_{n-1})
那么A(x)=a_0+a_1*x+a_2*{x^2}+a_3*{x^3}+a_4*{x^4}+a_5*{x^5}+ dots+a_{n-2}*x^{n-2}+a_{n-1}*x^{n-1}
将其下标按照奇偶性分类
A(x)=(a_0+a_2*{x^2}+a_4*{x^4}+dots+a_{n-2}*x^{n-2})+(a_1*x+a_3*{x^3}+a_5*{x^5}+ dots+a_{n-1}*x^{n-1})
设
A_1(x)=a_0+a_2*{x}+a_4*{x^2}+dots+a_{n-2}*x^{frac{n}{2}-1}
A_2(x)=a_1+a_3*{x}+a_5*{x^2}+ dots+a_{n-1}*x^{frac{n}{2}-1}
那么不难得到
A(x)=A_1(x^2)+xA_2(x^2)
我们将omega_n^k (k<frac{n}{2}) 代入得
A(omega_n^k)=A_1(omega_n^{2k})+omega_n^kA_2(omega_n^{2k})
=A_1(omega_{frac{n}{2}}^{k})+omega_n^kA_2(omega_{frac{n}{2}}^{k})
同理,将omega_n^{k+frac{n}{2}}代入得
A(omega_n^{k+frac{n}{2}})=A_1(omega_n^{2k+n})+omega_n^{k+frac{n}{2}}(omega_n^{2k+n})
=A_1(omega_n^{2k}*omega_n^n)-omega_n^kA_2(omega_n^{2k}*omega_n^n)
=A_1(omega_n^{2k})-omega_n^kA_2(omega_n^{2k})
大家有没有发现什么规律?
没错!这两个式子只有一个常数项不同!
那么当我们在枚举第一个式子的时候,我们可以O(1)的得到第二个式子的值
又因为第一个式子的k在取遍[0,frac{n}{2}-1]时,k+frac{n}{2}取遍了[frac{n}{2},n-1]
所以我们将原来的问题缩小了一半!
而缩小后的问题仍然满足原问题的性质,所以我们可以递归的去搞这件事情!
直到多项式仅剩一个常数项,这时候我们直接返回就好啦
时间复杂度:
不难看出FFT是类似于线段树一样的分治算法。
因此它的时间复杂度为O(nlogn)
快速傅里叶逆变换
不要以为FFT到这里就结束了。
我们上面的讨论是基于点值表示法的。
但是在平常的学习和研究中很少用点值表示法来表示一个多项式。
所以我们要考虑如何把点值表示法转换为系数表示法,这个过程叫做傅里叶逆变换
(y_0,y_1,y_2,dots,y_{n-1})$为$(a_0,a_1,a_2,dots,a_{n-1})的傅里叶变换(即点值表示)
设有另一个向量(c_0,c_1,c_2,dots,c_{n-1})满足
c_k=sum_{i=0}^{n-1}y_i(omega_n^{-k})^i
即多项式B(x)=y_0,y_1x,y_2x^2,dots,y_{n-1}x^{n-1}在omega_n^{0},omega_n^{-1},omega_n^{-2},dots,omega_{n-1}^{-(n-1)}处的点值表示
emmmm又到推公式时间啦
(c_0,c_1,c_2,dots,c_{n-1})满足 c_k=sum_{i=0}^{n-1}y_i(omega_n^{-k})^i
=sum_{i=0}^{n-1}(sum_{j=0}^{n-1}a_j(omega_n^i)^j)(omega_n^{-k})^i
=sum_{i=0}^{n-1}(sum_{j=0}^{n-1}a_j(omega_n^j)^i)(omega_n^{-k})^i
=sum_{i=0}^{n-1}(sum_{j=0}^{n-1}a_j(omega_n^j)^i(omega_n^{-k})^i)
=sum_{i=0}^{n-1}sum_{j=0}^{n-1}a_j(omega_n^j)^i(omega_n^{-k})^i
=sum_{i=0}^{n-1}sum_{j=0}^{n-1}a_j(omega_n^{j-k})^i
=sum_{j=0}^{n-1}a_j(sum_{i=0}^{n-1}(omega_n^{j-k})^i)
设S(x)=sum_{i=0}^{n-1}x^i
将omega_n^k代入得
S(omega_n^k)=1+(omega_n^k)+(omega_n^k)^2+dots(omega_n^k)^{n-1}
当k!=0时
等式两边同乘omega_n^k得
omega_n^kS(omega_n^k)=omega_n^k+(omega_n^k)^2+(omega_n^k)^3+dots(omega_n^k)^{n}
两式相减得
omega_n^kS(omega_n^k)-S(omega_n^k)=(omega_n^k)^{n}-1
S(omega_n^k)=frac{(omega_n^k)^{n}-1}{omega_n^k-1}
S(omega_n^k)=frac{(omega_n^n)^{k}-1}{omega_n^k-1}
S(omega_n^k)=frac{1-1}{omega_n^k-1}
观察这个式子,不难看出它分母不为0,但是分子为0
因此,当K!=0时,S(omega^{k}_{n})=0
那当k=0时呢?
很显然,S(omega^{0}_{n})=n
继续考虑刚刚的式子
c_k=sum_{j=0}^{n-1}a_j(sum_{i=0}^{n-1}(omega_n^{j-k})^i) 当j neq k时,值为0 当j=k时,值为n 因此, c_k=na_k a_k=frac{c_k}{n}
这样我们就得到点值与系数之间的表示啦
理论总结
至此,FFT的基础理论部分就结束了。
我们来小结一下FFT是怎么成功实现的
首先,人们在用系数表示法研究多项式的时候遇阻
于是开始考虑能否用点值表示法优化这个东西。
然后根据复数的两条性质(这个思维跨度比较大)得到了一种分治算法。
最后又推了一波公式,找到了点值表示法与系数表示法之间转换关系。
emmmm
其实FFT的实现思路大概就是
系数表示法—>点值表示法—>系数表示法
引用一下远航之曲大佬的图
当然,再实现的过程中还有很多技巧
我们根据代码来理解一下
递归实现
递归实现的方法比较简单。
就是按找我们上面说的过程,不断把要求的序列分成两部分,再进行合并
在c++的STL中提供了现成的complex类,但是我不建议大家用,毕竟手写也就那么几行,而且万一某个毒瘤卡STL那岂不是很GG?
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int MAXN=2*1e6+10;
inline int read()
{
char c=getchar();int x=0,f=1;
while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
return x*f;
}
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 fast_fast_tle(int limit,complex *a,int type)
{
if(limit==1) return ;//只有一个常数项
complex a1[limit>>1],a2[limit>>1];
for(int i=0;i<=limit;i+=2)//根据下标的奇偶性分类
a1[i>>1]=a[i],a2[i>>1]=a[i+1];
fast_fast_tle(limit>>1,a1,type);
fast_fast_tle(limit>>1,a2,type);
complex Wn=complex(cos(2.0*Pi/limit) , type*sin(2.0*Pi/limit)),w=complex(1,0);
//Wn为单位根,w表示幂
for(int i=0;i<(limit>>1);i++,w=w*Wn)//这里的w相当于公式中的k
a[i]=a1[i]+w*a2[i],
a[i+(limit>>1)]=a1[i]-w*a2[i];//利用单位根的性质,O(1)得到另一部分
}
int main()
{
int N=read(),M=read();
for(int i=0;i<=N;i++) a[i].x=read();
for(int i=0;i<=M;i++) b[i].x=read();
int limit=1;while(limit<=N+M) limit<<=1;
fast_fast_tle(limit,a,1);
fast_fast_tle(limit,b,1);
//后面的1表示要进行的变换是什么类型
//1表示从系数变为点值
//-1表示从点值变为系数
//至于为什么这样是对的,可以参考一下c向量的推导过程,
for(int i=0;i<=limit;i++)
a[i]=a[i]*b[i];
fast_fast_tle(limit,a,-1);
for(int i=0;i<=N+M;i++) printf("%d ",(int)(a[i].x/limit+0.5));//按照我们推倒的公式,这里还要除以n
return 0;
}
这里还有一个听起来很装B的优化—蝴蝶操作
观察合并的过程,w*a2[i] 这一项计算了两次,因为理论上来说复数的乘法是比较慢的,所以我们可以把这一项记出来
for(int i=0;i<(limit>>1);i++,w=w*Wn)//这里的w相当于公式中的k
{
complex t=w*a2[i];//蝴蝶操作
a[i]=a1[i]+t,
a[i+(limit>>1)]=a1[i]-t;//利用单位根的性质,O(1)得到另一部分
}
woc? 脸好疼。。。。。。
咳咳。
速度什么的才不是关键呢?
关键是我们AC不了啊啊啊
表着急,AC不了不代表咱们的算法不对,只能说这种实现方法太low了
下面介绍一种更高效的方法
迭代实现
再盗一下那位大佬的图
观察一下原序列和反转后的序列?
聪明的你有没有看出什么显而易见的性质?
没错!
我们需要求的序列实际是原序列下标的二进制反转!
因此我们对序列按照下标的奇偶性分类的过程其实是没有必要的
这样我们可以$O(n)$的利用某种操作得到我们要求的序列,然后不断向上合并就好了
// luogu-judger-enable-o2
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int MAXN=1e7+10;
inline int read()
{
char c=getchar();int x=0,f=1;
while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
return x*f;
}
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);}//不懂的看复数的运算那部分
int N,M;
int l,r[MAXN];
int limit=1;
void fast_fast_tle(complex *A,int type)
{
for(int i=0;i<limit;i++)
if(i<r[i]) swap(A[i],A[r[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=A[j+k],y=w*A[j+mid+k];//蝴蝶效应
A[j+k]=x+y;
A[j+mid+k]=x-y;
}
}
}
}
int main()
{
int N=read(),M=read();
for(int i=0;i<=N;i++) a[i].x=read();
for(int i=0;i<=M;i++) b[i].x=read();
while(limit<=N+M) limit<<=1,l++;
for(int i=0;i<limit;i++)
r[i]= ( r[i>>1]>>1 )| ( (i&1)<<(l-1) ) ;
// 在原序列中 i 与 i/2 的关系是 : i可以看做是i/2的二进制上的每一位左移一位得来
// 那么在反转后的数组中就需要右移一位,同时特殊处理一下奇数
fast_fast_tle(a,1);
fast_fast_tle(b,1);
for(int i=0;i<=limit;i++) a[i]=a[i]*b[i];
fast_fast_tle(a,-1);
for(int i=0;i<=N+M;i++)
printf("%d ",(int)(a[i].x/limit+0.5));
return 0;
}
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- 微信小程序的Web API接口设计及常见接口实现
- @陈同学的专属Python教程之常见数据结构
- 用易语言写个简单的小爬虫其中的关键点
- python坐标获取经纬度或经纬度获取坐标免费模块--geopy
- 详解:小程序页面预加载优化,让你的小程序运行如飞
- RocketMQ学习六-消息存储
- swoole 实现 unixSocket 通信
- mybatis-plus一对多关联查询踩坑
- 深入Spring Security魔幻山谷-获取认证机制核心原理讲解
- 文本相似性的总结
- Java面试题总结之JDBC 和Hibernate
- Mac 下搭建 Clion + OpenCV4.x 的开发环境
- 超详细,Windows系统搭建Flink官方练习环境
- MySQL 覆盖索引与延迟关联
- Java面试题总结之数据结构、算法和计算机基础(刘小牛和丝音的爱情故事1)