学习:多项式算法----FFT

FFT,即快速傅里叶变换,是离散傅里叶变换的快速方法,可以在很低复杂度内解决多项式乘积的问题(两个序列的卷积

 

1|0卷积



 

卷积通俗来说就一个公式(本人觉得卷积不重要)

Ck=i+j=kAiBi

 那么这个表达式是啥意思了:

  有两个序列AB,其中A={A1,A2,...}B={B1,B2,...}

  A序列有a个元素,B序列有b个元素。于是,由这两个序列可以推出另一个序列C,C序列如何确定了?确定方法就按照卷积的公式得来的,即C={i+k=0AiBj,i+k=1AiBj,...,i+k=2a+2bAiBj}

 

这里只介绍一下卷积,下面来探究卷积和多项式之间的关系。


 

 

 

2|0多项式(预备知识)



 

2|1多项式的定义


形如

f(x)=a0+a1x+a2x2+...+anxn

的函数关系式叫做关于x的多项式,其中多项式系数可以构成一个序列,即

A={a0,a1,a2,...,an}

 

 

2|2多项式乘法与序列的卷积


假如有两个多项式,其中

f(x)=a0+a1x+a2x2+...+anxn

g(x)=b0+b1x+b2x2+...+bmxm

 

现在要求f(x)*g(x)的序列,很明显

f(x)g(x)=i+j=0aibj+i+j=1aibjx+i+j=2aibjx2+...+i+j=m+naibjxm+n

 

现在把f(x),g(x)以及f(x)g(x)三个多项式的系数拿出来写成三个序列,可得:

A={a0,a1,a2,...an}

B={b0,b1,b2,...bm}

C={i+j=0aibj,i+j=1aibj,i+j=2aibj,...i+j=m+naibj}

 

于是惊讶的发现,两个序列的卷积相当于两个多项式乘积得到的多项式系数序列。而FFT算法的目的,就是求两个多项式乘积得到的多项式系数序列

 

 

2|3多项式的表示方法


多项式有两种表示方法,即系数表示法点值表示法

 

1.系数表示法是我们常用的表示方法,即

f(x)=a0+a1x+a2x2+...+anxn

的形式

 

2.点值表示法,顾名思义,就是在这个多项式上任意选取不重复的n+1个点,即

f(x)={(x0,y0),(x1,y1),...,(x2,y2)}

可以证明:任何n+1个点可以确定唯一一个最高次项为n次方的多项式(下面是证明,可以看看,也可以忽略)

将一个多项式的系数写成一个系数矩阵(当然,这些系数我们是不知道的)

[a0a1a2an]

然后将刚刚选取的n+1个点写成下面两个矩阵

[1x0x02x03x0n1x1x12x13x1n1x2x22x23x2n1xnxn2xn3xnn][y0y1y2yn]

可以得出以下三个矩阵的关系

[1x0x02x03x0n1x1x12x13x1n1x2x22x23x2n1xnxn2xn3xnn][a0a1a2an]=[y0y1y2yn]

解出系数矩阵

[a0a1a2an]=[1x0x02x03x0n1x1x12x13x1n1x2x22x23x2n1xnxn2xn3xnn]1[y0y1y2yn]

只能唯一确定一个系数矩阵,即只能唯一确定一组系数,即只能唯一确定一个多项式,证明完毕

 

 

2|4多项式乘法与FFT算法


为了实现多项式乘法并得到系数序列,我们可以考虑实现的方法,如果直接暴力(通过定义直接算系数)是O(n2)复杂度,肯定会超时,于是我们这样考虑

 

首先可以选取n+1个点,把两个多项式转化为点值表示法

 

然后将两个点值表示法的多项式相乘(复杂度为O(n)),然后将新得到的多项式的点值表示法转化为系数表示法

 

注意:假设f(x)最高有n次方g(x)最高有m次方,所以f(x)*g(x)最高有m+n次方,但是m+n可能不是2的幂次方,如果不是,则需要选取点的数量应该是大于m+n2的幂次方个,假设这个值为1<<k,所以从系数表示法点值表示法的转化过程中,我们要在f(x)和g(x)内选择1<<k的点才能保证点值表示法的f(x)和g(x)之后,至少仍有m+n个点才能确定唯一一个f(x)*g(x)的多项式。所以在选取点的数量的时候,要保证点的个数能确定f(x)*g(x)后的多项式。

 

如果任意的选取n+m个点然后转化,复杂度还是太高,所以我们需要巧妙的选取1<<k个点,从而使得系数表示法与点值表示法之间转化的复杂度降为O(nlogn),这就是FFT算法。


 

 

 

3|0虚数的乘积及其表示(预备知识)



 

一个虚数,是由实部和虚部组成,表示为(a+bi),虚数的几何表示可以在坐标系中体现,如下图所示

 

如图所示,r为此虚数的模长,θ为此虚数的幅角,于是虚数还能表示为r(cosθ,isinθ)的形式

 

z1=r1(cosθ1,isinθ1),z2=r2(cosθ2,isinθ2)

于是z1z2=r1(cosθ1,isinθ1)r2(cosθ2,isinθ2)=r1r2(cos(θ1+θ2),isin(θ1+θ2))

 

可以得到虚数相乘的几何意义:模长相乘,幅角相加


 

 

 

4|0n次单位根(预备知识)



 

4|1定义


n次单位根wn,即xn=1在虚数范围内的解,即wnn=1,故n次单位根wn也是一个复数,且模长为1所以n次单位根在乘方的时候模长不变,幅角等差增大

 

如图所示,将一个单位圆分成n块幅角为正最小的那个虚数即为n次单位根wn1红色点),下图是n=8的情况

还要注意,wn0=wnn=1

由图及n次单位根的性质(n次单位根模长为1,幅角为2πn)可得

wnk=cos2kπn+isin2kπn

 

 

4|2性质


比较简单的性质:

1.wnm+n=wnmwnn=wnm

2.(wnm)n=(wnn)m=1,m[0,n1]

 

比较复杂(重要)的性质:

1.w2n2k=wnk,如图所示

2.wnk=wnk+n2

通过上面那个图也可以看出来,注意一点:虚数乘-1,那么虚数的实部与虚部都要乘上-1,所以wnk在单位圆上的点和wnk在单位圆上的点关于原点对称

3.(wnk)2=wn2k

推导:(wnk)2=(wnk+n2)2=wn2k=wn2k

 

更加复杂(重要)的性质

k=0n1wnik={0imodn0nimodn=0

证明:

  当imodn0时,相当如等比数列前n项和求和,故k=0n1wnik=wn0(1(wni)n)1wni=1(wnn)i1wni=0

  当imodn=0时,wnik=1,故k=0n1wnik=n1=n

 

至此,你应该知道,多项式从系数表示法转化为点值表示法时,选取的n个点的x值即为n次单位根序列{wn0,wn1,....,wnn1}

 

注意,此处的n是f(x)*g(x)的多项式最高次项的次数,不是f(x)或者g(x)最高此项的系数,前面(多项式乘法与FFT)讲过要选取1<<k(再次强调这个值比f(x)*g(x)的多项式最高次的次数m+n还要大,且是2的幂次方)个点


 

 

 

5|0FFT



 

假设现在多项式为

f(x)=a0+a1x+a2x2+...+ak1xk1,注意,多项式最高次项为k-1次方,共k项,k如果不等于2的n次幂,就凑成2的n次幂(加系数为0的项)

g(x)=b0+b1x+b2x2+...+bh1xh1,注意,多项式最高次项为h-1次方,共h项,h如果不等于2的n次幂,就凑成2的n次幂(加系数为0的项)

 

那么现在,我们从多项式中选取了这样一些点

(wn0,f(wn0)),(wn1,f(wn1)),(wn2,f(wn2)),...,(wnn1,f(wnn1)),共n个点

注意n是大于k+h的最小2次幂,而不是等于k+h!!!!!

 

时选取得点数n大于k+h,故可以确定f(x)*g(x)这个多项式,现在,只需确定f(wn0),f(wn1),...,f(wnn1)的值即可。

 

如何算值?我们可以观察一下f(wni)的值(0i<n)

f(wni)=a0+a1wni+a2(wni)2+a3(wni)3+...+ak1(wni)k1

=a0+a2(wni)2+a4(wni)4+...+ak2(wni)k2+wni(a1+a3(wni)2+...+ak1(wni)k2)

根据(wnk)2=wn2k可得

=a0+a2wn2i+a4(wn2i)2+...+ak2(wn2i)k22+wni(a1+a3wn2i+...+ak1(wn2i)k22)

=f1(wni)+wnif2(wni)

 

再算f1(wni)的值和f2(wni)的值的时候,我们又可以按照奇偶分开,像算f(wni)一样变成两个问题

 

同时我们发现

f(wni+n2)=f(wni)

=a0+a2(wni)2+a4(wni)4+...+ak2(wni)k2wni(a1+a3(wni)2+...+ak1(wni)k2)

=a0+a2wn2i+a4(wn2i)2+...+ak2(wn2i)k22wni(a1+a3wn2i+...+ak1(wn2i)k22)

=f1(wni)wnif2(wni)

 

又发现,如果我们要算f(wni)的值,我们会先算f1(wni)f2(wni)的值,但是算出f1(wni)f2(wni)的值之后,不仅能算出f(wni)的值,还能同时算出f(wni+n2)的值

 

总的来说,如果n=8算出了f(wn0)的值,就算出了f(wn4)的值;算出了f(wn1),就算出来f(wn5)的值;...;算出了f(wn3)的值,就算出了f(wn7)的1值

 

同理,算f1(wni)f2(wni)的值的时候,也是算出一半,得到另一半的值------每次解决问题,都会变成解决一半问题然后直接得到另一半的答案,在解决这一半问题的时候,也变成了解决一半的一半的问题,另一半的一半的问题又被直接得到答案

 

这样,就可以把f(wn0),f(wn1),f(wn2),...,f(wnn1)的值全部算出来,得到了n个点的值;同理算出g(wn0),g(wn1),f(wn2),...,g(wnn1)的值。得到了f(x)和g(x)的点值表示法,点值表示法相乘,O(n)复杂度就得到了f(x)*g(x)的点值表示法

 

代码如下:

#include <complex>//c++自带复数模板 typedef complex<double> Complex;//重定义数据类型,记得要用double类型 void FFT(Complex *a,int n,int inv){ //a是你要进行变换的数组(多项式系数数组),n是数组大小,inv暂时不管,你就当做它等于1 if(n==1) return;//如果分治后的序列大小为1,就直接返回,不用继续分治了 int mid=n>>1;//如果n不等于1,那就继续分治,mid是分治后变成两个序列的长度 static Complex b[1000];//创建一个辅助用的b数组,后面体现用处 for(int i=0;i<mid;i++) b[i]=a[2*i],b[i+mid]=a[2*i+1];//将a数组奇数位置的值和偶数位置的值分开,b数组前mid个位置存a数组偶数位置值,后mid个位置存a数组奇数位置值 for(int i=0;i<n;i++) a[i]=b[i];//将b数组值重新赋给a数组 //对分治后的两个序列(一个是a数组前mid个元素组成的序列,另一个数a数组后mid元素组成的序列) 进行FFT变换 FFT(a,mid,inv); FFT(a+mid,mid,inv); //算出一半,得到另一半的值,a数组前mid元素相当于之前说的f1,a数组后mid元素相当于之前说的f2 for(int i=0;i<mid;i++){ Complex com(cos(2*pi*i/n),inv*sin(2*pi*i/n)); b[i]=a[i]+com*a[i+mid];b[i+mid]=a[i]-com*a[i+mid]; } for(int i=0;i<n;i++) a[i]=b[i];//重新将b数组赋值给a数组 return; }
FFT变换(递归版)

 

这里再来详细解释一下,上面代码中的n值

加入开始两个序列是n次多项式和m次多项式,那么开始n值就等于大于等于(m+n)的最小2次幂,如下代码所示

cin>>n>>m; int k=n+m,lim=1; while(k>>=1) lim<<=1; if(lim<n+m) lim<<=1; FFT(a,lim,1);//对a序列进行变换 FFT(b,lim,1);//对b序列进行变换

此时lim值就是n的初值

 

强调一遍wnk=cos2kπn+isin2kπn


 

 

 

6|0FFT逆变换



 

FFT逆变换则是将点值表示法转化为系数表示法的步骤。在逆变换之前,我们已经对两个多项式系数序列进行了正变换,即

FFT(a,lim,1);//正变换 FFT(b,lim,1); for(int i=0;i<lim;i++) a[i]=a[i]*b[i];//将正变换的两个序列乘在一起,得到新的多项式的点值表示法序列。

 

 

如上所示,只要将a序列进行逆变换,就可以得到新多项式的系数序列。在FFT正变换之中,我们一直在求f(wni)的值。

假设所求多项式系数序列为k0,k1,...kn1

 

现在对于新的多项式a(新得到的序列,不是原来的a序列)来求卷积ch=i=0n1ai(wnh)i

相当于对a序列又进行FFT变换,只不过,我们现在以新序列a为系数的多项式在x=wni的一系列值

相当于求C(x)=a0+a1x+a2x2+...+an1xn1x=wn0,wn1,wn2,...wn(n1)的点值表示法

而我们知道ax=k0+k1wnx+k2(wnk)2+...=j=0n1kj(wnx)j

 

i=0n1ai(wnh)i

=i=0n1(j=0n1kj(wni)j)(wnh)i

=j=0n1kj(i=0n1(wnjh)i)

=j=0n1kj(i=0n1wn(jh)i)

 

由于当且只当j=h时,(jh)才为i的倍数,其余时候为0,根据n次单位根更加复杂的性质可得(不知道的看看n次单位根预备性质

j=0n1kj(i=0n1wn(jh)i)

=nkh

所以

kh=i=0n1ai(wnh)ih=chn

 

说明:wniwni实部相同,虚部相反。现在知道代码中inv的作用了吧

当inv=1,表示FFT变换;当inv=-1,表示FFT逆变换

 

FFT(a,lim,1); FFT(b,lim,1); for(int i=0;i<lim;i++) a[i]=a[i]*b[i]; FFT(a,lim,-1); //经过FFT逆变换之后,a序列虚部都已经为0,只有实部有值,且为整数 for(int i=0;i<n+m-1;i++) printf("%d ",(int)(a[i].real()/lim+0.5));//实部根据前面推导要除以n(这里除以lim)

为啥要加0.5捏:首先我们确定答案肯定是整数,但是在代码实现过程中,由于有π介入计算,导致答案可能变小了(误差很小),所以加一个0.5,然后强制类型转换切掉小数部

比如说本来答案是1,但是代码实现的结果可能是0.999999999999,此时加上了0.5,然后转换为int,答案就是1了。

 

强调一遍wnk=cos2kπn+isin2kπn

FFT代码

#include <complex>//c++自带复数模板 typedef complex<double> Complex;//重定义数据类型,记得要用double类型 void FFT(Complex *a,int n,int inv){ //a是你要进行变换的数组(多项式系数数组),n是数组大小,inv=1表示正变换,inv=-1表示逆变换 if(n==1) return;//如果分治后的序列大小为1,就直接返回,不用继续分治了 int mid=n>>1;//如果n不等于1,那就继续分治,mid是分治后变成两个序列的长度 static Complex b[1000];//创建一个辅助用的b数组,后面体现用处 for(int i=0;i<mid;i++) b[i]=a[2*i],b[i+mid]=a[2*i+1];//将a数组奇数位置的值和偶数位置的值分开,b数组前mid个位置存a数组偶数位置值,后mid个位置存a数组奇数位置值 for(int i=0;i<n;i++) a[i]=b[i];//将b数组值重新赋给a数组 //对分治后的两个序列(一个是a数组前mid个元素组成的序列,另一个数a数组后mid元素组成的序列) 进行FFT变换 FFT(a,mid,inv); FFT(a+mid,mid,inv); //算出一半,得到另一半的值,a数组前mid元素相当于之前说的f1,a数组后mid元素相当于之前说的f2 for(int i=0;i<mid;i++){ Complex com(cos(2*pi*i/n),inv*sin(2*pi*i/n)); b[i]=a[i]+com*a[i+mid];b[i+mid]=a[i]-com*a[i+mid]; } for(int i=0;i<n;i++) a[i]=b[i];//重新将b数组赋值给a数组 return; }
FFT变换

 


 

 

 

 

7|0FFT系数储存注意



 

在进行FFT之前,我们有两个多项式系数序列,用数组存系数的不要用普通的数组存,而要用复数数组,因为在计算的时候系数涉及复数计算,所以这样

#include <complex>//c++自带复数模板 Complex a[1000],b[1000]; for(int i=0;i<n;i++) scanf("%d",&t),a[i]=Complex(t,0);//构造函数 for(int i=0;i<m;i++) scanf("%d",&t),b[i]=Complex(t,0);

这样,复数的实部存系数的值,虚部为0

 

一个算多项式相乘后系数序列的程序

//开头说明一下,最下面有输入输出样例 #include <iostream> #include <cmath> #define pi acos(-1) #include <complex>//c++自带复数模板 using namespace std; typedef complex<double> Complex;//重定义数据类型,记得要用double类型 Complex a[1000],b[1000]; int n,m; void FFT(Complex *a,int n,int inv){ //a是你要进行变换的数组(多项式系数数组),n是数组大小,inv=1表示正变换,inv=-1表示逆变换 if(n==1) return;//如果分治后的序列大小为1,就直接返回,不用继续分治了 int mid=n>>1;//如果n不等于1,那就继续分治,mid是分治后变成两个序列的长度 static Complex b[1000];//创建一个辅助用的b数组,后面体现用处 for(int i=0;i<mid;i++) b[i]=a[2*i],b[i+mid]=a[2*i+1];//将a数组奇数位置的值和偶数位置的值分开,b数组前mid个位置存a数组偶数位置值,后mid个位置存a数组奇数位置值 for(int i=0;i<n;i++) a[i]=b[i];//将b数组值重新赋给a数组 //对分治后的两个序列(一个是a数组前mid个元素组成的序列,另一个数a数组后mid元素组成的序列) 进行FFT变换 FFT(a,mid,inv); FFT(a+mid,mid,inv); //算出一半,得到另一半的值,a数组前mid元素相当于之前说的f1,a数组后mid元素相当于之前说的f2 for(int i=0;i<mid;i++){ Complex com(cos(2*pi*i/n),inv*sin(2*pi*i/n)); b[i]=a[i]+com*a[i+mid];b[i+mid]=a[i]-com*a[i+mid]; } for(int i=0;i<n;i++) a[i]=b[i];//重新将b数组赋值给a数组 return; } int main(){ cin>>n>>m;//输入两个多项式的项数 int k=n+m,lim=1,t; for(int i=0;i<n;i++) scanf("%d",&t),a[i]=Complex(t,0);//第一个多项式系数 for(int i=0;i<m;i++) scanf("%d",&t),b[i]=Complex(t,0);//第二个多项式系数 while(k>>=1) lim<<=1; if(lim<n+m) lim<<=1; FFT(a,lim,1); FFT(b,lim,1); for(int i=0;i<lim;i++) a[i]=a[i]*b[i]; FFT(a,lim,-1); for(int i=0;i<n+m-1;i++) printf("%d ",(int)(a[i].real()/lim+0.5)); } /* //输入 2 3//第一个多项式有两项 ,第二个多项式有三项 1 2//第一个多项式为(1+2x) 2 1 2//第二个多项式为(2+x+2x^2) //输出 2 5 4 4//(1+2x)*(2+x+2x^2)=2+5x+4x^2+4x^3,系数为2 5 4 4 */
多项式相乘

测试用例

//输入 2 2//两个多项式都有两项 1 1//第一个多项式为(1+x) 1 1//第二个多项式为(1+x) //输出 1 2 1//(1+x)*(1+x)=1+2*x+x^2,系数为1 2 1
测试

 


 

 

 

8|0FFT变换的示例演示



 

在讲FFT迭代版本之前,我打算先用一个例子来展示FFT变换,加深认识

假设

  多项式f(x)=3+2x+x2

  多项式g(x)=2+x+2x2

在草稿纸上算一下,这两个多项式相乘的结果为h(x)=6+7x+10x2+5x3+2x4

 

先得到f(x)g(x)多项式的系数序列,注意应该写成复数形式的序列(复数形式为(real,image))

  多项式f(x)的系数序列为A={(3,0),(2,0),(1,0)}

  多项式g(x)的系数序列为B={(2,0),(1,0),(2,0)}

 

然后判断两个多项式相乘后有多少项,应该不超过3+3=6项。而大于6的最小2次幂应该为8,所以要把两个多项式序列补成8项。

 

(用(0,0)来补项)

  多项式f(x)的系数序列为A={(3,0),(2,0),(1,0),(0,0),(0,0),(0,0),(0,0),(0,0)}

  多项式g(x)的系数序列为B={(2,0),(1,0),(2,0),(0,0),(0,0),(0,0),(0,0),(0,0)}

 

对于每一个f(w8i),可以寻得规律

f(w8i)=A0+A1w8i+A2(w8i)2+A3(w8i)3+A4(w8i)4+A5(w8i)5+A6(w8i)6+A7(w8i)7

=A0+A2(w8i)2+A4(w8i)4+A6(w8i)6+w8i(A1+A3(w8i)2+A5(w8i)4+A7(w8i)6)

=A0+A2w4i+A4(w4i)2+A6(w4i)3+w8i(A1+A3w4i+A5(w4i)2+A7(w4i)3)

=A0+A4(w4i)2+w4i(A2+A6(w4i)2)+w8i(A1+A5(w4i)2+w4i(A3+A7(w4i)2))

=A0+A4w2i+w4i(A2+A6w2i)+w8i(A1+A5w2i+w4i(A3+A7w2i))

 

然后

f(w8i+4)=f(w8i)=A0+A1(w8i)+A2(w8i)2+A3(w8i)3+A4(w8i)4+A5(w8i)5+A6(w8i)6+A7(w8i)7

=A0+A2(w8i)2+A4(w8i)4+A6(w8i)6w8i(A1+A3(w8i)2+A5(w8i)4+A7(w8i)6)

 

只要算出来A0+A2w4i+A4(w4i)2+A6(w4i)3A1+A3w4i+A5(w4i)2+A7(w4i)3的值,就可以算出f(w8i)f(w8i+4)的值

 

同理,算出A0+A4(w4i)2A2+A6w2i的值,不仅可以得到A0+A2w4i+A4(w4i)2+A6(w4i)3的值,还可以得到A0+A2w4i+2+A4(w4i+2)2+A6(w4i+2)3

 

然后将每一个序列分成长度为2的小序列(下面是用1~8来编号的)

 

后面的是重点,对于每一个长度为2的小序列(设为(x,y)):

  第一个元素x重新赋值为x+w10y

  第二个元素y重新赋值为xw10y

(这个w11这样来的:$w_{当前序列一半长度l}^{一半长度l的第i个元素,i从0到l-1}$)

重新赋值完之后就可以合并,每两个长度为2的相邻序列合并为一个长度为4的序列

 

对于每一个长度为4的序列(设为(a,b,c,d)):

  第一个元素a重新赋值为a+w20c

  第三个元素c重新赋值为aw20c

  第二个元素b重新赋值为b+w21d

  第四个元素d重新赋值为bw21d

重新赋值完之后就可以合并,每两个长度为4的相邻序列合并为一个长度为8的序列(此时全部小序列已经合并为一个序列了)

 

对于整个序列(长度为8),按照上面的赋值方法重新赋值(可以参考代码),于是FFT变换就ok了

两个序列就变成了

  A={A0,A1,A2,A3,A4,A5,A6,A7}

  B={B0,B1,B2,B3,B4,B5,B6,B7}

将两个序列对应元素相乘,得到新序列

  C={A0B0,A1B1,A2B2,A3B3,A4B4,A5B5,A6B6,A7B7}

然后将C序列进行FFT逆变换,就是多项式相乘后的系数序列了。

 

ps:FFT逆变换就是把当前序列再来一次FFT变换,只不过在重新赋值环节的时候是乘上wmid1,而不是wmidi,最后得到的序列虚部为0,将实部全部除以8即可(这里不懂参考代码或前面逆变换的讲解)

 


 

 

 

9|0FFT迭代版本



 

9|1FFT的序列更新规律


综上所述,FFT的序列更新是有规律的,先把序列分解为序列长度为2的小段每一段更新完毕后就合并小段,继续一组一组更新,再合并,最后合并成一个序列,然后再最后一次更新完毕。

 

如图所示:

最后序列{A0,A1,...,An1}为原序列进过FFT变换后的序列。

 

不像递归版本,迭代版本是直接从每一个小区间开始更新,然后合并,然后更新...,但是我们需要找到原系数序列经过奇偶分离的一系列操作后每个值的位置发生了什么变化。

 

假设原系数序列为{A0,A1,A2,A3,A4,A5,A6,A7},下图是奇偶分离的示意图

 

没错,这不是巧合,序列经过奇偶分离之后,前后位置数的二进制是对称的,根据这一点我们就不需要用递归来实现,可以直接来得到奇偶分离最终序列,然后再进行FFT

 

 

9|2FFT代码


FFT迭代版比递归版快的不是一点。。

关于序列奇偶分离可以用下面方法

int rev[(n+2)>>1]; for(i=0;i<n;i++){ rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); if (i<rev[i])swap(a[i],a[rev[i]]); }

 

于是就可以实现迭代版本FFT

void FFT(Complex *a,int n,int inv){ int bit=0,i,rev[(n+2)>>1]; while((1<<bit)<n) bit++; //交换序列 for(i=0;i<n;i++){ rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); if (i<rev[i])swap(a[i],a[rev[i]]);//if来保证只存在小的和大的交换,不会发生大的和小的交换 } for (int mid=1;mid<n;mid*=2){//mid枚举区间长度的一半 Complex com(cos(pi/mid),inv*sin(pi/mid));//单位根 for (int i=0;i<n;i+=mid*2){//i来枚举每个区间的开头位置 Complex omega(1,0); for (int j=0;j<mid;j++,omega*=com){//j枚举每个区间的前一半元素 Complex x=a[i+j],y=omega*a[i+j+mid]; a[i+j]=x+y,a[i+j+mid]=x-y;//蝴蝶效应 } } } }
FFT变换(迭代版)

 


 

 

 

10|0总接代码



 

FFT递归版本

void FFT(Complex *a,int n,int inv){ //a是你要进行变换的数组(多项式系数数组),n是数组大小,inv=1表示正变换,inv=-1表示逆变换 if(n==1) return;//如果分治后的序列大小为1,就直接返回,不用继续分治了 int mid=n>>1;//如果n不等于1,那就继续分治,mid是分治后变成两个序列的长度 static Complex b[1000];//创建一个辅助用的b数组,后面体现用处 for(int i=0;i<mid;i++) b[i]=a[2*i],b[i+mid]=a[2*i+1];//将a数组奇数位置的值和偶数位置的值分开,b数组前mid个位置存a数组偶数位置值,后mid个位置存a数组奇数位置值 for(int i=0;i<n;i++) a[i]=b[i];//将b数组值重新赋给a数组 //对分治后的两个序列(一个是a数组前mid个元素组成的序列,另一个数a数组后mid元素组成的序列) 进行FFT变换 FFT(a,mid,inv); FFT(a+mid,mid,inv); //算出一半,得到另一半的值,a数组前mid元素相当于之前说的f1,a数组后mid元素相当于之前说的f2 for(int i=0;i<mid;i++){ Complex com(cos(2*pi*i/n),inv*sin(2*pi*i/n));// b[i]=a[i]+com*a[i+mid];b[i+mid]=a[i]-com*a[i+mid]; // cout<<n<<" "<<i<<" "<<b[i].real()<<" "<<b[i].imag()<<endl; } for(int i=0;i<n;i++) a[i]=b[i];//重新将b数组赋值给a数组 return; }
FFT变换(迭代版)

 

FFT迭代版本

void FFT(Complex *a,int n,int inv){ //a是你要进行变换的数组(多项式系数数组),n是数组大小,inv=1表示正变换,inv=-1表示逆变换 int bit=0,i,rev[(n+2)>>1]; while((1<<bit)<n) bit++; //交换序列 for(i=0;i<n;i++){ rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); if (i<rev[i])swap(a[i],a[rev[i]]);//if来保证只存在小的和大的交换,不会发生大的和小的交换 } for (int mid=1;mid<n;mid*=2){//mid枚举区间长度的一半 Complex com(cos(pi/mid),inv*sin(pi/mid));//单位根 for (int i=0;i<n;i+=mid*2){//i来枚举每个区间的开头位置 Complex omega(1,0); for (int j=0;j<mid;j++,omega*=com){//j枚举每个区间的前一半元素 Complex x=a[i+j],y=omega*a[i+j+mid]; a[i+j]=x+y,a[i+j+mid]=x-y;//蝴蝶效应 } } } }
FFT变换(迭代版)

 

 

 

 

 

 

 

 

 


__EOF__

本文作者七月之流
本文链接https://www.cnblogs.com/qiyueliu/p/11237318.html
关于博主:评论和私信会在第一时间回复。或者直接私信我。
版权声明:本博客所有文章除特别声明外,均采用 BY-NC-SA 许可协议。转载请注明出处!
声援博主:如果您觉得文章对您有帮助,可以点击文章右下角推荐一下。您的鼓励是博主的最大动力!
posted @   七月流  阅读(3427)  评论(2编辑  收藏  举报
编辑推荐:
· DeepSeek 解答了困扰我五年的技术问题
· 为什么说在企业级应用开发中,后端往往是效率杀手?
· 用 C# 插值字符串处理器写一个 sscanf
· Java 中堆内存和栈内存上的数据分布和特点
· 开发中对象命名的一点思考
阅读排行:
· PPT革命!DeepSeek+Kimi=N小时工作5分钟完成?
· What?废柴, 还在本地部署DeepSeek吗?Are you kidding?
· DeepSeek企业级部署实战指南:从服务器选型到Dify私有化落地
· 程序员转型AI:行业分析
· 重磅发布!DeepSeek 微调秘籍揭秘,一键解锁升级版全家桶,AI 玩家必备神器!
点击右上角即可分享
微信分享提示