ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

FFT

2022-02-24 10:02:35  阅读:182  来源: 互联网

标签:stdin stdout int FFT ch freopen define


建议全文背诵

3->2优化

#include<bits/stdc++.h>
using namespace std;
#define LL long long 
#define file(a) freopen(#a".in","r",stdin),freopen(#a".out","w",stdout)
#define endless() fclose(stdin),fclose(stdout)
#define ReadOnly(a) freopen(#a,"r",stdin)
inline int read(){
	int s=0,f=1;char ch=getchar();
	while(ch<'0'||'9'<ch) {if(ch=='-') f=-1;ch=getchar();}
	while('0'<=ch&&ch<='9') {s=s*10+(ch^48);ch=getchar();}
	return s*f;
}
const int N=1e6+5e5+3;
const double pi=acos(-1);
struct Complex{
	double x,y;
	Complex(double _x=0,double _y=0){
		x=_x;y=_y;
	}
	inline friend Complex operator+(Complex a,Complex b){
		return Complex(a.x+b.x,a.y+b.y);
	} 
	inline friend Complex operator-(Complex a,Complex b){
		return Complex(a.x-b.x,a.y-b.y);
	}
	inline friend Complex operator*(Complex a,Complex b){
		return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
	}
}F[N<<1];
int rev[N<<1];int n,m;
inline void get_rev(){
	for(int i=0;i<n;++i){
		rev[i]=(rev[i>>1]>>1)|((i&1)?n>>1:0);
	}
}
inline void FFT(Complex *f,bool flag){
	for(int i=0;i<n;++i){
		if(i<rev[i]) swap(f[i],f[rev[i]]);
	}
	for(int p=2;p<=n;p<<=1){//枚举区间长度
		int len=p>>1;//获取分治长度
		Complex w1n(cos(2*pi/p),sin(2*pi/p) );
		if(!flag) w1n.y*=-1;
		for(int k=0;k<n;k+=p){//枚举初始化位置
			Complex wxn(1,0);
			for(int l=k;l<k+len;++l){//枚举detail
				Complex tmp=wxn*f[len+l];
				f[len+l]=f[l]-tmp;
				f[l]=f[l]+tmp;
				wxn=wxn*w1n;//获取下一个单位根
			}
		}
	}
}
int main(void){
	n=read();m=read();
	for(int i=0;i<=n;++i){
		F[i].x=read();
	}
	for(int i=0;i<=m;++i){
		F[i].y=read();
	}
	for(m+=n,n=1;n<=m;n<<=1);//补足2的幂次
	get_rev();//二进制翻转
	FFT(F,1);//DFT
	for(int i=0;i<n;++i) F[i]=F[i]*F[i];
	FFT(F,0);//IDFT
	for(int i=0;i<=m;++i){
		printf("%d ",(int)(F[i].y/n/2+0.49) );
	}
	return 0;
}

原版,更加灵活

#include<bits/stdc++.h>
using namespace std;
#define LL long long 
#define file(a) freopen(#a".in","r",stdin),freopen(#a".out","w",stdout)
#define endless() fclose(stdin),fclose(stdout)
#define ReadOnly(a) freopen(#a,"r",stdin)
inline int read(){
	int s=0,f=1;char ch=getchar();
	while(ch<'0'||'9'<ch) {if(ch=='-') f=-1;ch=getchar();}
	while('0'<=ch&&ch<='9') {s=s*10+(ch^48);ch=getchar();}
	return s*f;
}
const int N=1e6+5e5+3;
const double pi=acos(-1);
struct Complex{
	double x,y;
	Complex(double _x=0,double _y=0){
		x=_x;y=_y;
	}
	inline friend Complex operator+(Complex a,Complex b){
		return Complex(a.x+b.x,a.y+b.y);
	} 
	inline friend Complex operator-(Complex a,Complex b){
		return Complex(a.x-b.x,a.y-b.y);
	}
	inline friend Complex operator*(Complex a,Complex b){
		return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
	}
}F[N<<1],P[N<<1];
int rev[N<<1];int n,m;
inline void get_rev(){
	for(int i=0;i<n;++i){
		rev[i]=(rev[i>>1]>>1)|((i&1)?n>>1:0);
	}
}
inline void FFT(Complex *f,bool flag){
	for(int i=0;i<n;++i){
		if(i<rev[i]) swap(f[i],f[rev[i]]);
	}
	for(int p=2;p<=n;p<<=1){//枚举区间长度
		int len=p>>1;//获取分治长度
		Complex w1n(cos(2*pi/p),sin(2*pi/p) );
		if(!flag) w1n.y*=-1;
		for(int k=0;k<n;k+=p){//枚举初始化位置
			Complex wxn(1,0);
			for(int l=k;l<k+len;++l){//枚举detail
				Complex tmp=wxn*f[len+l];
				f[len+l]=f[l]-tmp;
				f[l]=f[l]+tmp;
				wxn=wxn*w1n;//获取下一个单位根
			}
		}
	}
}
int main(void){
	n=read();m=read();
	for(int i=0;i<=n;++i){
		F[i].x=read();
	}
	for(int i=0;i<=m;++i){
		P[i].x=read();
	}
	for(m+=n,n=1;n<=m;n<<=1);//补足2的幂次
	get_rev();//二进制翻转
	FFT(F,1);FFT(P,1);//DFT
	for(int i=0;i<n;++i) F[i]=F[i]*P[i];
	FFT(F,0);//IDFT
	for(int i=0;i<=m;++i){
		printf("%d ",(int)(F[i].x/n+0.49) );
	}
	return 0;
}

标签:stdin,stdout,int,FFT,ch,freopen,define
来源: https://www.cnblogs.com/cbdsopa/p/15930371.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有