fft求调
  • 板块学术版
  • 楼主zzyaba
  • 当前回复0
  • 已保存回复0
  • 发布时间2024/7/23 09:57
  • 上次更新2024/7/23 10:50:16
查看原帖
fft求调
360025
zzyaba楼主2024/7/23 09:57
#include<bits/stdc++.h>
#include<windows.h>
#define MAXN 16
using namespace std;
const int sz=1<<MAXN;
struct z{
	double a,b;
	z(double _a=0,double _b=0):a(_a),b(_b){
	}
	z operator=(double _a){
		a=_a;
		b=0;
	}
	z operator+(z o){
		return z(a+o.a,b+o.b);
	}
	z operator-(z o){
		return z(a-o.a,b-o.b);
	}
	z operator*(z o){
		return z(a*o.a-b*o.b,a*o.b+b*o.a);
	}
	z operator/(z o){
		return z((a*o.b-b*o.a)/(o.a*o.a+o.b*o.b),(a*o.a+b*o.b)/(o.a*o.a+o.b*o.b));
	}
	z operator-(){
		return z(-a,-b);
	}
	z operator!(){
		return z(a,-b);
	}
}o,qpow[sz],t[4];
ostream&operator<<(ostream&x,z b){
	if(abs(b.a)>=1e-6){
		if(abs(b.b)>=1e-6){
			if(b.b>0){
				x<<b.a<<"+"<<b.b<<"i";
			}else{
				x<<b.a<<"-"<<-b.b<<"i";
			}
		}else{
			x<<b.a;
		}
	}else{
		if(abs(b.b)>=1e-6){
			x<<b.b<<"i";
		}else{
			x<<0;
		}
	}
	return x;
}
int n,m,l;
double a[sz],b[sz],pre[sz];
z x[sz],y[sz],l1,l2;
void FFT(double*a,z*b,int len){
	if(len==1){
		*b=*a;
		return;
	}
	FFT(a,b,len/2);
	FFT(a+len/2,b+len/2,len/2);
	for(int i=0;i<len/2;i++){
		l1=b[i];
		l2=b[len/2+i];
		b[i]=l1+l2*qpow[i*sz/len%sz];
		b[len/2+i]=l1-l2*qpow[i*sz/len%sz];
	}
}
void fft(double*a,z*b,int len){
	for(int i=0;i<len;i++){
		l=0;
		for(int j=0;j<MAXN;j++){
			if(i&(1<<j)){
				l+=(1<<MAXN-j-1);
			}
		}
		pre[l]=a[i];
	}
	FFT(pre,b,len);
}
void DFFT(double*a,z*b,int len){
	if(len==1){
		*a=(*b*t[MAXN%4]).a;
		return;
	}
	for(int i=0;i<len/2;i++){
		l1=b[i];
		l2=b[len/2+i];
		b[i]=(l1+l2)/2;
		b[len/2+i]=(l1-l2)/2*qpow[(sz-i*sz/len)%sz];
	}
	DFFT(a,b,len/2);
	DFFT(a+len/2,b+len/2,len/2);
}
void dfft(double*a,z*b,int len){
	DFFT(pre,b,len);
	for(int i=0;i<len;i++){
		l=0;
		for(int j=0;j<MAXN;j++){
			if(i&(1<<j)){
				l+=(1<<MAXN-j-1);
			}
		}
		a[i]=pre[l];
	}
}
int main(){
	o=z(cos(acos(-1.)*2/sz),sin(acos(-1.)*2/sz));
	qpow[0]=1;
	for(int i=1;i<sz;i++){
		qpow[i]=qpow[i-1]*o;
	}
	t[0]=1;
	for(int i=1;i<4;i++){
		t[i]=t[i-1]*z(0,-1);
	}
	cin>>n>>m;
	for(int i=0;i<=n;i++){
		cin>>a[i];
	}
	for(int i=0;i<=m;i++){
		cin>>b[i];
	}
	fft(a,x,sz);
	fft(b,y,sz);
	for(int i=0;i<sz;i++){
		x[i]=x[i]*y[i];
	}
	dfft(a,x,sz);
	for(int i=0;i<=n+m;i++){
		cout<<a[i]<<" ";
	}
}

当MAXN>16时输出错误

2024/7/23 09:57
加载中...