> 科技资讯 > 快速傅里叶变换FFT学习笔记

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

科技资讯 投稿 5600℃ 0 评论

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

点值表示法

我们知道，\(n+1\个点可以表示出一个 \(n\ 次的多项式。

复数

一种表示坐标的方法，对于坐标 \((x,y\,可写作 \(x+iy\，其中\(x\为实部，\(y\为虚部。

complex，可以直接作为变量类型使用。

如果你不会，可以看看百度百科。

我的表述自然不够专业，希望可以表述出这个意思吧。

Tips：

为啥是这样呢？证明如下:

\[(a+bi(c+di=ac-bd+i(bc+ad \]
\[\sqrt{(ac-bd^2+(bc+ad^2}\\ =\sqrt{(ac^2-2abcd+(bd^2+(bc^2+2abcd+(ad^2}\\ =\sqrt{(ac^2+(bd^2+(bc^2+(ad^2}\\ =\sqrt{(a^2+b^2(c^2+d^2}\\ \]

那么我们应该可以看出来这个模长相乘了。

\[\cos(\theta_1+\theta_2=\cos(\theta_1\cos(\theta_2-\sin(\theta_1\sin(\theta_2\\ =\frac {a}{t_1}\times\frac {c}{t_2}-\frac {b}{t_1}\times\frac {d}{t_2}\\ =\frac {ac-bd}{t_1t_2} \]

我们知道分母是乘积的模长，分子是横坐标，所以这个式子恰好就是乘积辐角对应的 \(\cos\ 值。

把圆均分

我们定义\(\omega_n^k\表示从\((1,0\开始，把圆均分为\(n\份的第\(k\个点的复数，其中\((1,0\为\(\omega_n^0\。

这是为什么呢？我们考虑它的辐角，由于其平分了一整个圆，所以其辐角为 \(360\frac k n°\，转换为弧度后则为 \(2\pi\frac k n\，且模长为 \(1\，利用三角函数易得其坐标了。

    性质1：\(\omega _{dn}^{dk}=\omega_n^k\，根据定义可证。
  • 性质2：\(\omega _{n}^{k+\frac n 2}=-\omega_n^k\，两点对称。

这个东西有什么用呢？

离散傅里叶变换

对于一个多项式\(A(x=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\,我们按照前文所云，将所有的 \(\omega_n^k\作为 \(x\ 代入。

这被称为 \(A(x\ 的傅里叶变换。

傅里叶逆变换

对于这个多项式\(B(x\，代入所有的 \(\omega_n^{-k}\，也就是\(\omega_n^k\的倒数，得到 \((z_0,z_1,z_2,...,z_{n-1}\。

\[z_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_{j=0}^{n-1}a_j(\omega_n^i^j\sum_{i=0}^{n-1}(\omega_n^{-k}^i\\ =\sum_{j=0}^{n-1}a_j(\omega_n^j^i\sum_{i=0}^{n-1}(\omega_n^{-k}^i\\ =\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k}^i\\ \]

关于后面那个等比数列，若 \(j=k\，可得 \(1\，否则用等比数列式子可知为 \(0\。

\[z_k=na_k\\ a_k=\frac {z_k} n \]

所以我们可以求出原来的多项式了。

快速傅里叶变换

我们可以使用分治的思想，使得时间复杂度降至\(O(n\log n\。

\[A_1(x=a_0+a_2x+...,A_2(x=a_1+a_3x+...\\ A(x=A_1(x^2+xA_2(x^2\\ A(\omega_n^k=A_1(\omega_n^{2k}+\omega_n^kA_2(\omega_n^{2k}\\ A(\omega_n^k=A_1(\omega_{\frac n 2}^{k}+\omega_n^kA_2(\omega_{\frac n 2}^{k}\\ \]

同理:

\[A(\omega_n^{k+\frac n 2}=A_1(\omega_n^{2k+n}+\omega_n^{k+\frac n 2}A_2(\omega_n^{2k+n}\\ A(\omega_n^{k+\frac n 2}=A_1(\omega_{\frac n 2}^{k}-\omega_n^{k}A_2(\omega_{\frac n 2}^{k}\\ \]

不过注意这里一直有一个除二操作，为了方便，我们需要把多项式补成一个次数为\(2^x-1\ 的多项式。

注意这个取倒可以利用共轭复数，对于 \(a+ib\，其共轭复数为 \(a-ib\。

\[\frac {1}{a+ib}=\frac {a-ib}{(a+ib(a-ib}=\frac {a-ib}{a^2+b^2} \]
void FFT(cp *a,LL n,bool inv//inv 表示omega是否取倒
{
	if(n==1return;
	static cp buf[N];
	LL m=n/2;
	for(int i=0;i<=m-1;i++buf[i]=a[2*i],buf[i+m]=a[2*i+1];//奇偶分开
	for(int i=0;i<=n-1;i++a[i]=buf[i];
	FFT(a,m,inv,FFT(a+m,m,inv;
	for(int i=0;i<=m-1;i++
	{
		cp x=omega(n,i;
		if(invx=conj(x;//conj(x求解共轭复数
		buf[i]=a[i]+x*a[i+m],buf[i+m]=a[i]-x*a[i+m];
	}
	for(int i=0;i<=n-1;i++a[i]=buf[i];
}

利用FFT求解多项式的乘积

这个还是十分简单的，直接把两个多项式转化为两个长度相同的次数为\(2^x-1\的多项式。

为什么呢？

因此多项式相乘以后，我们希望 \(a_i=A(\omega_n^iB(\omega_n^i\，那么就是直接相乘喽。

、#include<bits/stdc++.h>
#define LL int
#define cp complex<double>
using namespace std;
const double PI=acos(-1.0000;
const int N=5e6+5;
cp omega(LL n, LL k
{
    return cp(cos(2*PI*k/n,sin(2*PI*k/n;
}
LL n,x,len,ans[N];
cp a[N],b[N];
//上文的 FFT 实现省去
int main(
{
	scanf("%d",&n; 
	len=1;
	while(len<2*nlen*=2;
	for(int i=n-1;i>=0;i--scanf("%1d",&x,a[i].real(x;
	for(int i=n-1;i>=0;i--scanf("%1d",&x,b[i].real(x;
	FFT(a,len,0,FFT(b,len,0;
	for(int i=0;i<=len-1;i++a[i]*=b[i];
	FFT(a,len,1;
	for(int i=0;i<=len-1;i++//进位
	{
		ans[i]+=floor(a[i].real(/len+0.5;
		ans[i+1]+=ans[i]/10;
		ans[i]%=10;
	}
	int i=len;
	for(i;i>=0&&ans[i]==0;i--;//前导零
	if(i==-1len=0;
	for(;i>=0;i--printf("%d",ans[i]; 
}

非递归FFT

这里有一个优化，我们发现每次递归有一个把 \(a_i\ 奇偶分开的过程，本质来看，就是将二进制末尾为 \(0\ 的数字与二进制末尾为 \(1\ 的数字分开。

我们先将每个 \(a_i\ 放置在对应的位置，然后向上逐渐合并。

void FFT(cp *a,bool inv
{
	LL lim=0;
	while((1<<lim<lenlim++;
	for(int i=0;i<=len-1;i++
	{
		LL t=0;
		for(int j=0;j<lim;j++
			if((i>>j&1t|=(1<<(lim-j-1;//处理其翻转后的值
		if(i<tswap(a[i],a[t];
	}
	static cp buf[N];
	for(int l=2;l<=len;l*=2
	{
		LL m=l/2;
		for(LL j=0;j<=len-1;j+=l
		{
			for(LL i=0;i<=m-1;i++
			{
				cp x=omega(l,i+j;
				if(invx=conj(x;
				buf[i+j]=a[i+j]+x*a[i+j+m];
				buf[i+j+m]=a[i+j]-x*a[i+j+m];
			}
		}
		for(int j=0;j<=len-1;j++a[j]=buf[j];
	}
}

蝴蝶操作

这个东西其实就是想了个办法使得把工具人数组 buf 除掉了。

void FFT(cp *a,bool inv
{
	LL lim=0;
	while((1<<lim<lenlim++;
	for(int i=0;i<=len-1;i++
	{
		LL t=0;
		for(int j=0;j<lim;j++
			if((i>>j&1t|=(1<<(lim-j-1;
		if(i<tswap(a[i],a[t];
	}
	for(int l=2;l<=len;l*=2
	{
		LL m=l/2;
		for(LL j=0;j<=len-1;j+=l
		{
			for(LL i=0;i<=m-1;i++
			{
				cp x=omega(l,i+j;
				if(invx=conj(x;
				x*=a[i+j+m];
				a[i+j+m]=a[i+j]-x;
				a[i+j]=a[i+j]+x;
			}
		}
	}
}

一些小优化

对于 \(i\ 的二进制翻转可以先预处理出来。

#include<bits/stdc++.h>
#define LL int
#define cp complex<double>
using namespace std;
const double PI=acos(-1.0000;
const int N=5e6+5;
cp omega(LL n, LL k
{
    return cp(cos(2*PI*k/n,sin(2*PI*k/n;
}
LL n,len,lim,x,ans[N],r[N];
cp a[N],b[N];
void FFT(cp *a,bool inv
{
	for(int i=0;i<=len-1;i++
	{
		LL t=r[i];
		if(i<tswap(a[i],a[t];
	}
	for(int l=2;l<=len;l*=2
	{
		LL m=l/2;
		cp omg=omega(l,1;
		for(LL j=0;j<=len-1;j+=l
		{
			cp x(1,0; 
			for(LL i=0;i<=m-1;i++
			{
				cp t=x;
				if(invt=conj(t;
				t*=a[i+j+m];
				a[i+j+m]=a[i+j]-t,a[i+j]=a[i+j]+t;
				x*=omg;
			}
		}
	}
}
int main(
{
	scanf("%d",&n; 
	len=1;
	while(len<2*nlen*=2,lim++;
	for(int i=0;i<=len-1;i++
	{
		LL t=0;
		for(int j=0;j<lim;j++if((i>>j&1t|=(1<<(lim-j-1;
		r[i]=t;
	}
	for(int i=n-1;i>=0;i--scanf("%1d",&x,a[i].real(x;
	for(int i=n-1;i>=0;i--scanf("%1d",&x,b[i].real(x;
	FFT(a,0,FFT(b,0;
	for(int i=0;i<=len-1;i++a[i]*=b[i];
	FFT(a,1;
	for(int i=0;i<=len-1;i++
	{
		ans[i]+=floor(a[i].real(/len+0.5;
		ans[i+1]+=ans[i]/10;
		ans[i]%=10;
	}
	int i=len;
	for(i;i>=0&&ans[i]==0;i--;
	if(i==-1len=0;
	for(;i>=0;i--printf("%d",ans[i];
}

参考

胡小兔-小学生都能看懂的FFT！

编程笔记 » 快速傅里叶变换FFT学习笔记

赞同 (29) or 分享 (0)
游客 发表我的评论   换个身份
取消评论

表情
(0)个小伙伴在吐槽