6
17
2016
0

【新姿势】再谈快速傅里叶变换(FFT)

之前学过但一直没有用它来写过题(太弱),这次再来复习一下~

快速傅里叶变换能够快速计算卷积。

离散傅里叶变换(DFT)可以求循环卷积,N,M为两类的个数,当循环卷积长度L>=N+M-1,就可以做线性卷积了。

快速傅里叶变换就是加速这一过程的。

FFT理解:

http://blog.csdn.net/acdreamers/article/details/39005227

http://blog.csdn.net/acdreamers/article/details/39026505

目的是加速多项式乘法(计算系数),而多项式乘积是卷积的形式,因为相同次项的系数为原来下标之和为该值的两个系数之积。所以该算法可以推广至快速求所有卷积形式的式子(系数)。

注意,最高次项为n,m的两个多项式,其乘积的最高次项为m+n;


代码理解与注意事项:

1、

#include<complex>中的赋值操作:

typedef complex<double> Ex;//之后会出现double数

Ex a(real,imag);

int x; scanf("%d",&x);//只要变量类型与%d一致即可,不一定要是double

Ex a=x;//单参数赋值给实部,虚部为0

2、 m=m+n;for(n=1;n<=m;n<<=1) len++;//找比m+n大的最小2幂次,len记录该指数。

3、 F(i,0,n-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));//由i>>1的逆位序置换递推i的你逆位序置换的过程,注意前导0/后缀0的处理.

4、(int)(a[i].real()+0.5)//四舍五入而不是默认的去尾,减小精度误差。理论上讲IDFT结果是整数。

5、

const double pi=acos(-1.0);//因为cos(pi)=-1;注意acos不要写成cos!

const int N=(int)1e6;//数组开2*N往往是不够的,会RE!

6、FFT代码是用迭代方式模拟递归归并的过程。

具体见代码注释:

void FFT(Ex *a,int n,int on)
{
	F(i,0,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);//防止重复交换
	for(int i=1;i<n;i<<=1){//枚举当前操作组(归并后的组)的个数(1-->2,n/2-->n) 
		Ex wn(cos(pi/i),on*sin(pi/i));
/*旋转因子 wn(cos(on*2*pi/i),sin(on*2*pi/i))
同时约掉2,并根据三角函数奇偶性将on提出*/ 
		for(int j=0;j<n;j+=i*2){//枚举起始元素标号(每个组2i个元素) 
			Ex w(1.0,0);//n次单位复数根 
			for(int k=0;k<i;k++){//枚举每次操作的每个元素 
				Ex x=a[j+k],y=a[j+k+i]*w;//两小组合并为当前组的过程 
				a[j+k]=x+y,a[j+k+i]=x-y;
				w=w*wn;//折半引理&&相消引理 更新旋转因子
			}
		}
	}
	if(on==-1) for(int i=0;i<n;i++) a[i].real()/=n;//归一系数,将点值转化为系数 
}

比如说:

1、大整数乘法,高精度相乘复杂度为O(M*N),但由于满足c[i+j]=c[i]+c[j](mod 10),可知乘法就是做线性卷积的过程!可以用FFT来加速,达到O(NlogN).

bzoj2179: FFT快速傅 立叶

【假装有代码】

2、多项式乘法(得到系数),满足卷积的基本形式:c[k]=∑(a[i]*b[k-i])下标和为定值

UOJ#34: 多项式乘法

#include<bits/stdc++.h>
#include<complex>
#define F(i,s,t) for(int i=(s);i<=(t);i++)
using namespace std;
int n,m,len;
typedef complex<double> Ex;
const double pi=acos(-1.0);
const int N=(int)1e6;
int rev[N]; 
Ex a[N],b[N],c[N];
void FFT(Ex *a,int n,int on)
{
	F(i,0,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int i=1;i<n;i<<=1){
		Ex wn(cos(pi/i),on*sin(pi/i));
		for(int j=0;j<n;j+=i*2){
			Ex w(1.0,0);
			for(int k=0;k<i;k++){
				Ex x=a[j+k],y=a[j+k+i]*w;
				a[j+k]=x+y,a[j+k+i]=x-y;
				w=w*wn;
			}
		}
	}
	if(on==-1) for(int i=0;i<n;i++) a[i].real()/=n;  
}
int main()
{
//	freopen("fft2.txt","r",stdin); 
	scanf("%d%d",&n,&m);int x;
	F(i,0,n) scanf("%d",&x),a[i]=x;
	F(i,0,m) scanf("%d",&x),b[i]=x;
	m=m+n;for(n=1;n<=m;n<<=1) len++;
	F(i,0,n-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	FFT(a,n,1);FFT(b,n,1);
	F(i,0,n) a[i]=a[i]*b[i];
	FFT(a,n,-1);
	F(i,0,m) printf("%d ",(int)(a[i].real()+0.5)); 
	return 0;
}
当然,形如c[k]=∑(a[i]*b[i-k])的形式可以将一个数组翻转,再用FFT.
 
bzoj2194: 快速傅立叶之二
 
Category: FFT | Tags: 学习笔记 卷积 | Read Count: 889

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter

Host by is-Programmer.com | Power by Chito 1.3.3 beta | Theme: Aeros 2.0 by TheBuckmaker.com