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

点值表示法

我们正常表示一个多项式的方式,形如 \(A(x)=a_0+a_1x+a_2x^2+...+a_nx^n\),这是正常人容易看懂的,但是,我们还有一种表示法。

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

于是,我们任意地取 \(n+1\) 个不同的值,表示 \(x\) ,求出的值与 \(x\) 对应,形成 \(n+1\)个点,这就可以表示。

复数

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

C++中有复数的模板,complex,可以直接作为变量类型使用。

运算规则,自然不用多说了,也就是直接拿式子算即可。

如果你不会,可以看看百度百科

我们将点至原点的距离称为模长,将其与原点相连之后与 \(x\) 轴形成的一个夹角称为辐角,不过呢,对于第三第四象限的点,自然要加上一个 \(180°\)了。

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

复数的乘法可以理解为,模长相乘,辐角相加。

Tips:

此部分的证明来自 cjx 犇犇。

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

\[(a+bi)(c+di)=ac-bd+i(bc+ad) \]

那么这个点就是 \((ac-bd,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)}\\ \]

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

接下来是辐角相加,我们设原来两个辐角为 \(\theta_1,\theta_2\),而模长为 \(t_1,t_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\) 值。

那么显然,辐角的值就是 \(\theta_1+\theta_2\)了。

把圆均分

img

如图是坐标轴上一个以原点为圆心的半径为 \(1\) 的圆。

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

那么,不难发现,\(\omega_n^k\)表示的点为 \((\cos(2\pi\frac k n),\sin(2\pi\frac n k))\)

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

简单推推式子不难发现 \((\omega_n^1)^k=\omega_n^k\),这是利用了模长相乘,辐角相加,因为模长是 \(1\) ,怎么乘都是 \(1\),于是辐角不断叠加,从定义上看是这样的。

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

这个东西有什么用呢?

离散傅里叶变换

离散傅里叶变换(Discrete Fourier Transform,简称DFT)的思想是利用 \(\omega_n^k\)将一个多项式转为点值表示法。

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

于是我们得到了 \(n-1\) 个点,使用复数形式表示,成为一个数组 \((y_0,y_1,y_2,...,y_{n-1})\)的。

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

傅里叶逆变换

我们再将其作为一个多项式 \(B(x)=y_0+y_1x+....+y_{n-1}x^{n-1}\)

对于这个多项式\(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 \]

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

有几点需要注意:

  • 我们提取 \(a_i\) 只提取实部,因为虚部是虚数,无法经过我们的转换后变成实数。
  • 数字有一定误差(毕竟你使用了三角函数等东西,小数是会有误差的),所以要四舍五入。
for(int i=0;i<=len-1;i++)
{
	ans[i]+=floor(a[i].real()/len+0.5);
}

快速傅里叶变换

快速傅里叶变换(Fast Fourier Transform,简称FFT),是在 DFT的基础上我们发现时间复杂度依然需要 \(O(n^2)\),没有含金量,所以我们要给他含金量!

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

对于\(A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\),我们可以:

\[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})\\ \]

不错!利用这两个式子,我们可以在 \(O(n \log n)\) 的时间复杂度求出 \(A(x)\)

不过注意这里一直有一个除二操作,为了方便,我们需要把多项式补成一个次数为\(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} \]

由于此处 \(a^2+b^2=1\) ,因此其共轭复数为其倒数。

void FFT(cp *a,LL n,bool inv)//inv 表示omega是否取倒
{
	if(n==1)return;
	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(inv)x=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\) 不是一个点,而是单纯因为 \(\omega_n^i\)是一个复数,所以 \(a_i\) 才变成了一个复数,\(a_i\)只是表示\(A(\omega_n^i)\)的值。

首先我们知道当前的 \(a_i\)表示的是 \(A(\omega_n^i)\)\(b_i\)表示的是 \(B(\omega_n^i)\),那么我们将两个值直接相乘。

因此多项式相乘以后,我们希望 \(a_i=A(\omega_n^i)B(\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*n)len*=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==-1)len=0;
	for(;i>=0;i--)printf("%d",ans[i]); 
}

非递归FFT

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

我们不妨想一下,对于一个数 \(x\),其位置可以根据其二进制确定,就是其二进制倒过来的数字。
img

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

void FFT(cp *a,bool inv)
{
	LL lim=0;
	while((1<<lim)<len)lim++;
	for(int i=0;i<=len-1;i++)
	{
		LL t=0;
		for(int j=0;j<lim;j++)
			if((i>>j)&1)t|=(1<<(lim-j-1));//处理其翻转后的值
		if(i<t)swap(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(inv)x=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)<len)lim++;
	for(int i=0;i<=len-1;i++)
	{
		LL t=0;
		for(int j=0;j<lim;j++)
			if((i>>j)&1)t|=(1<<(lim-j-1));
		if(i<t)swap(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(inv)x=conj(x);
				x*=a[i+j+m];
				a[i+j+m]=a[i+j]-x;
				a[i+j]=a[i+j]+x;
			}
		}
	}
}

一些小优化

对于 \(i\) 的二进制翻转可以先预处理出来,我们用 \(R_i\) 表示 \(i\) 的二进制翻转。

然后 \(\omega_n^k\) 可以利用性质累乘,会快很多,因为复数的各种操作常数巨大,直接累乘常数较小。

最后代码就长这样了:

#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<t)swap(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(inv)t=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*n)len*=2,lim++;
	for(int i=0;i<=len-1;i++)
	{
		LL t=0;
		for(int j=0;j<lim;j++)if((i>>j)&1)t|=(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==-1)len=0;
	for(;i>=0;i--)printf("%d",ans[i]);
}

参考

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

posted @ 2023-04-29 13:21  DengDuck  阅读(95)  评论(0编辑  收藏  举报