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

Posted on by hackdl

点值表示法

我们正常表示一个多项式的方式，形如 (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}timesfrac {c}{t_2}-frac {b}{t_1}timesfrac {d}{t_2}\
=frac {ac-bd}{t_1t_2}
]

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

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

把圆均分

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

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

那么，不难发现，(omega_n^k)表示的点为 ((cos(2pifrac k n),sin(2pifrac n k)))。

这是为什么呢？我们考虑它的辐角，由于其平分了一整个圆，所以其辐角为 (360frac k n°)，转换为弧度后则为 (2pifrac 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)^jsum_{i=0}^{n-1}(omega_n^{-k})^i\
=sum_{j=0}^{n-1}a_j(omega_n^j)^isum_{i=0}^{n-1}(omega_n^{-k})^i\
=sum_{j=0}^{n-1}a_jsum_{i=0}^{n-1}(omega_n^{j-k})^i\
]

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

因此：

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

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

快速傅里叶变换

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

我们可以使用分治的思想，使得时间复杂度降至(O(nlog 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

利用FFT求解多项式的乘积

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

求解出其傅里叶变换形式之后，对于每个点对应的复数，相乘即可。

为什么呢？

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

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

给一个简单的实现：

、#include
#define LL int
#define cp complex
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=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=0&&ans[i]==0;i--);//前导零
	if(i==-1)len=0;
	for(;i>=0;i--)printf("%d",ans[i]); 
}

非递归FFT

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

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

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

void FFT(cp *a,bool inv)
{
	LL lim=0;
	while((1>j)&1)t|=(1

蝴蝶操作

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

调调顺序即可。

void FFT(cp *a,bool inv)
{
	LL lim=0;
	while((1>j)&1)t|=(1

一些小优化

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

然后 (omega_n^k) 可以利用性质累乘，最后代码就长这样了：

#include
#define LL int
#define cp complex
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>j)&1)t|=(1=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=0&&ans[i]==0;i--);
	if(i==-1)len=0;
	for(;i>=0;i--)printf("%d",ans[i]);
}

参考

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

服务器托管，北京服务器托管，服务器租用 http://www.fwqtg.net
机房租用，北京机房租用，IDC机房托管， http://www.e1idc.net

服务器托管，北京服务器托管，服务器租用，机房机柜带宽租用

服务器托管

咨询：董先生

电话13051898268 QQ/微信93663045！

上一篇:
下一篇: