FFT 模板

FFT 板子

背板子的岁月又开始了

#include<cstdio>
#include<algorithm>
#include<queue>
#include<cstring>
#include<cmath>
#define r register
#define rep(i,x,y) for(r ll i=x;i<=y;++i)
#define per(i,x,y) for(r ll i=x;i>=y;--i)
using namespace std;
typedef long long ll;
const ll V=(1e7+10)/2;
const double pi=acos(-1.0);
ll in()
{
	ll res=0,f=1;
	char ch;
	while((ch=getchar())<'0'||ch>'9')
	 if(ch=='-') f=-1;
	res=res*10+ch-48;
	while((ch=getchar())>='0'&&ch<='9')
	 res=res*10+ch-48;
	return res*f;
}
ll n,m,f[V];
ll now,k;
struct complex 
{
	double x,y;
	complex (double xx=0,double yy=0) { x=xx,y=yy; }
}a[V],b[V];
complex operator + (complex a,complex b) { return complex(a.x+b.x,a.y+b.y); }
complex operator - (complex a,complex b) { return complex(a.x-b.x,a.y-b.y); }
complex operator * (complex a,complex b) { return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x); }
void FFT(complex *a,ll v)
{
	rep(i,0,now-1)
	 if(i<f[i]) swap(a[i],a[f[i]]);
	for(r ll mid=1;mid<now;mid<<=1)
	{
		complex Wn(cos(pi/mid),v*sin(pi/mid));
		ll R=mid<<1;
		for(r ll j=0;j<now;j+=R)
		{
			complex val(1,0);
			for(r ll k=0;k<mid;++k,val=val*Wn)
			{
				complex x=a[j+k],y=val*a[j+mid+k];
				a[j+k]=x+y;
				a[j+mid+k]=x-y;
			}
		}
	}
}
int main()
{
	scanf("%lld%lld",&n,&m);
	rep(i,0,n) a[i].x=in();
	rep(i,0,m) b[i].x=in();
	now=1;//单位根中的k值,默认为2的自然数次幂 
	while(now<=n+m) now<<=1,++k;
	rep(i,0,now-1)
	 f[i]=(f[i>>1]>>1)|((i&1)<<(k-1));
	FFT(a,1);
	FFT(b,1);
	rep(i,0,now) a[i]=a[i]*b[i];
	FFT(a,-1);
	rep(i,0,m+n)
	 printf("%lld ",(ll)(a[i].x/now+0.49));
	return 0;
}

上一篇:利用FFT 求频域功率 与时域平均功率 使用matlab验证


下一篇:分治fft