FFT(傅里叶快速变换，详细讲解+推导) 每日一遍，算法再见！

2023-05-16

FFT详细推导

FFT(傅里叶快速变换)
- 一.前置知识
- - 1.复数和单位根
  - 2.单位根的三个引理
  - 3.多项式
- 二.FFT(快速傅里叶变换推导)
- 三.IFFT
- 四.FFT求解多项式乘积模板代码
- - 1.递归版
  - 2.非递归版(这个更快，省去了递归时间)
- 五.视频资源
- 六.FFT题目集

FFT(傅里叶快速变换)

FFT在实际工程中有着非常的广泛，尤其是在信号领域，在ACM算法竞赛领域主要可以用来快速计算多项式的乘积

一.前置知识

1.复数和单位根

有人会觉得FFT怎么会用到复数的知识，想必是有点古怪，后面在推导过程中会介绍到它的用处。
在这里插入图片描述

2.单位根的三个引理

在这里插入图片描述

3.多项式

在这里插入图片描述

前置知识学完了之后下面我们来推导FFT

二.FFT(快速傅里叶变换推导)

DFT的作用是把多项式的系数表示法，转化为点值表示法,复杂度 O ( n 2 ) O(n^2) O(n2)，而FFT则是快速版的DFT，作用是一样的，FFT的复杂度是 O ( n l o g n ) O(nlogn) O(nlogn).
在这里插入图片描述

在这里插入图片描述

上面我们说过n次多项式需要n+1个点来唯一确定，那么我们找点的时候为什么带入的x是复数(也即wi)而不是实数呢？这个视频讲的很详细FFT推导过程。

FFT模板：

/*这里的opt=1*/
void FFT(Complex *a,int lim,int opt)
{
    if(lim==1) return ;
	Complex a0[lim>>1],a1[lim>>1];
	/*我们把多项式的奇数项和偶数项拆开*/ 
	for(int i=0;i<lim;i+=2)
	{
	 a0[i>>1]=a[i],a1[i>>1]=a[i+1];
    }
    FFT(a0,lim>>1,opt);//然后递归拆解奇数项 
    FFT(a1,lim>>1,opt);//递归拆解偶数项 
    Complex wn = Complex(cos(2.0*PI/lim),opt*sin(2.0*PI/lim));//单位根
	Complex w = Complex(1,0);//第一个根w0 
    for(int k=0;k<(lim>>1);k++)
    {
    	Complex t=w*a1[k];//这样会少一次复数运算； 
    	a[k]=a0[k]+t;//最后我们把求得的值代回 
    	a[k+(lim>>1)]=a0[k]-t;
    	w=w*wn;//想当于次幂+1 
	}
	return ; 
}

三.IFFT

我们要解决多项式a和多项式b的乘积运算，FFT把多项式转化为点值表示之后，我们对两个多项式进行乘积运算，得到一个新的多项式C的点值表示，我们还需要把C的点值表示转化为系数表示，这样才能得到每一项的系数，这时就需要用到FFT的逆过程，IFFT(快速傅里叶逆变换)，也就是IDFT(离散傅里叶逆变换)的快速版。
在这里插入图片描述

IFFT模板：

/*这里的opt=-1*/
void FFT(Complex *a,int lim,int opt)
{
    if(lim==1) return ;
	Complex a0[lim>>1],a1[lim>>1];
	/*我们把多项式的奇数项和偶数项拆开*/ 
	for(int i=0;i<lim;i+=2)
	{
	 a0[i>>1]=a[i],a1[i>>1]=a[i+1];
    }
    FFT(a0,lim>>1,opt);//然后递归拆解奇数项 
    FFT(a1,lim>>1,opt);//递归拆解偶数项 
    Complex wn = Complex(cos(2.0*PI/lim),opt*sin(2.0*PI/lim));//单位根
	Complex w = Complex(1,0);//第一个根w0 
    for(int k=0;k<(lim>>1);k++)
    {
    	Complex t=w*a1[k];//这样会少一次复数运算； 
    	a[k]=a0[k]+t;//最后我们把求得的值代回 
    	a[k+(lim>>1)]=a0[k]-t;
    	w=w*wn;//想当于次幂+1 
	}
	return ; 
}

通过FFT和IFFT的对比我们发现只有opt参数不一样其余地方全是相同的，FFT的opt=1，IFFT的opt=-1，而opt就只作用于下面这一段代码块。

Complex wn = Complex(cos(2.0*PI/lim),opt*sin(2.0*PI/lim));

这就是我们上面这张图的解释
在这里插入图片描述
最后a/n，因为a是一个Complex，所以a/n指的是a的实部除以n，也即a.r/n，因为我们的系数是存在Complex的实部里面的，虚部只是用于FFT和IFFT的递归运算过程。

四.FFT求解多项式乘积模板代码

1.递归版

下面给出整个代码：

#include<bits/stdc++.h>
using namespace std;
const double PI=acos(-1);//3.1415926... 
const int Maxn=1e5;
struct Complex
{
	double r,i;
	Complex() {r=0.0,i=0.0;}//不带参数的构造函数 
	Complex(double real,double imag):r(real),i(imag){}//带参数的构造函数 
}F[Maxn],G[Maxn];
Complex operator + (Complex a,Complex b){
	return Complex(a.r+b.r,a.i+b.i);
}
Complex operator - (Complex a,Complex b)
{
	return Complex(a.r-b.r,a.i-b.i);
}
Complex operator * (Complex a,Complex b)
{
	return Complex(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);
}
/*快速傅里叶变换*/
void FFT(Complex *a,int lim,int opt)
{
    if(lim==1) return ;
	Complex a0[lim>>1],a1[lim>>1];
	/*我们把多项式的奇数项和偶数项拆开*/ 
	for(int i=0;i<lim;i+=2)
	{
	 a0[i>>1]=a[i],a1[i>>1]=a[i+1];
    }
    FFT(a0,lim>>1,opt);//然后递归拆解奇数项 
    FFT(a1,lim>>1,opt);//递归拆解偶数项 
    Complex wn = Complex(cos(2.0*PI/lim),opt*sin(2.0*PI/lim));//单位根
	Complex w = Complex(1,0);//第一个根w0 
    for(int k=0;k<(lim>>1);k++)
    {
    	Complex t=w*a1[k];//这样会少一次复数运算； 
    	a[k]=a0[k]+t;//最后我们把求得的值代回 
    	a[k+(lim>>1)]=a0[k]-t;
    	w=w*wn;//想当于次幂+1 
	}
	return ; 
}  
int main()
{
	int n,m,lim=1;
	scanf("%d %d",&n,&m);
	for(int i=0;i<=n;i++) scanf("%lf",&F[i].r);
	for(int i=0;i<=m;i++) scanf("%lf",&G[i].r);
	int len = n+m;
    while(lim<=len) lim<<=1;// 
    FFT(F,lim,1);//我们先通过FFT，把a多项式的系数表示转化为点对表示 
	FFT(G,lim,1);//我们先通过FFT，也把b多项式的系数表示转化为点对表示
	
	/*这时候我们对a，b进行卷积操作，就算F,G没有lim那么多项，
	但是我们给复数默认传参是0+0i，所以卷积对结果没有影响 */ 
	for(int i=0;i<=lim;i++) F[i]=F[i]*G[i]; 
	
	FFT(F,lim,-1);//然后我们通过IFFT，也就是快速傅里叶逆变换，把点对表示转化为系数表示 
	for(int i=0;i<=n+m;i++) printf("%d ",(int)(F[i].r/lim+0.5));//最后的每一个系数/n,然后四舍五入就是答案 

}

我们这里测试一下
第一行我们输入n，m分别是3，3表示多项式a，多项式b的最高次次数。
第二行我们输入多项式a各项的系数(默认从低次到高次项)
1 2 3 4
第三行我们输入多项式b各项的系数(默认从低次到高次项)
1 2 3 4

我们的结果是
1 4 10 20 25 24 16
在这里插入图片描述
上面想当于解决了多项式A×多项式B的问题，也即 A ⨁ B A\bigoplus B A⨁B
其中 A ( x ) = 1 + 2 x + 3 x 2 + 4 x 3 A(x) = 1+2x+3x^2+4x^3 A(x)=1+2x+3x2+4x3, B ( x ) = 1 + 2 x + 3 x 2 + 4 x 3 B(x) = 1+2x+3x^2+4x^3 B(x)=1+2x+3x2+4x3
计算得到 C ( x ) = 1 + 4 x + 10 x 2 + 20 x 3 + 25 x 4 + 24 x 5 + 16 x 6 C(x) = 1+4x+10x^2+20x^3+25x^4+24x^5+16x^6 C(x)=1+4x+10x2+20x3+25x4+24x5+16x6

2.非递归版(这个更快，省去了递归时间)

非递归版需要用到蝴蝶效应的知识，本蒟蒻还没有学习完，待更…

五.视频资源

如果看了笔记还是不明白的话，这里推荐一个非常详细的视频讲解，看不懂直接来打我！
建议两个都仔细看一遍，有必要做一下笔记尝试推导
1.FFT的推导细节说明
2.全网最详细FFT,DFT,IDFT，IFFT讲解。

六.FFT题目集

1.模板题，这道题虽然可以之前可以python，现在数据加强了不知道题目建议用FFT，注意高精度进位就可以了，还有就是，数组开到5e7不然会re
P1919 【模板】A*B Problem升级版（FFT快速傅里叶）
AC代码：

#include<bits/stdc++.h>
using namespace std;
const double PI=acos(-1);//3.1415926... 
const int Maxn=5e6+5;
struct Complex
{
	double r,i;
	Complex() {r=0.0,i=0.0;}//不带参数的构造函数 
	Complex(double real,double imag):r(real),i(imag){}//带参数的构造函数 
}F[Maxn],G[Maxn];
Complex operator + (Complex a,Complex b){
	return Complex(a.r+b.r,a.i+b.i);
}
Complex operator - (Complex a,Complex b)
{
	return Complex(a.r-b.r,a.i-b.i);
}
Complex operator * (Complex a,Complex b)
{
	return Complex(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);
}
/*快速傅里叶变换*/
void FFT(Complex *a,int lim,int opt)
{
    if(lim==1) return ;
	Complex a0[lim>>1],a1[lim>>1];
	/*我们把多项式的奇数项和偶数项拆开*/ 
	for(int i=0;i<lim;i+=2)
	{
	 a0[i>>1]=a[i],a1[i>>1]=a[i+1];
    }
    FFT(a0,lim>>1,opt);//然后递归拆解奇数项 
    FFT(a1,lim>>1,opt);//递归拆解偶数项 
    Complex wn = Complex(cos(2.0*PI/lim),opt*sin(2.0*PI/lim));//单位根
	Complex w = Complex(1,0);//第一个根w0 
    for(int k=0;k<(lim>>1);k++)
    {
    	Complex t=w*a1[k];//这样会少一次复数运算； 
    	a[k]=a0[k]+t;//最后我们把求得的值代回 
    	a[k+(lim>>1)]=a0[k]-t;
    	w=w*wn;//想当于次幂+1 
	}
	return ; 
} 
char a[Maxn],b[Maxn]; 
int ans[Maxn];
int main()
{
	int n=0,m=0,lim=1;
	scanf("%s",a);
	scanf("%s",b);
	n=strlen(a)-1;
	m=strlen(b)-1;
	for(int i=n;i>=0;i--) F[n-i].r = a[i]-'0';
	for(int i=m;i>=0;i--) G[m-i].r = b[i]-'0';
	
	int len = n+m;
    while(lim<=len) lim<<=1;
    FFT(F,lim,1);//我们先通过FFT，把a多项式的系数表示转化为点对表示 
	FFT(G,lim,1);//我们先通过FFT，也把b多项式的系数表示转化为点对表示
	for(int i=0;i<=lim;i++) F[i]=F[i]*G[i]; 
	
	FFT(F,lim,-1);//然后我们通过IFFT，也就是快速傅里叶逆变换，把点对表示转化为系数表示 
	for(int i=0;i<=n+m;i++) ans[i]=(int)(F[i].r/lim+0.5);//最后的每一个系数/n,然后四舍五入就是答案 
 	int  num=0;

	for(int i=0;i<=n+m;i++)
 	{
 	   ans[i]+=num;
 	   num=ans[i]/10;
 	   ans[i]%=10;
	}
	while(num)
	{
		ans[++len] = num%10;
		num/=10;
	}
	for(int i=len;i>=0;i--) cout<<ans[i];
	return 0;
}

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)