[POJ 2821]TN's Kindom III(任意长度循环卷积的Bluestein算法)

[POJ 2821]TN's Kindom III(任意长度循环卷积的Bluestein算法)

题面

给出两个长度为(n)的序列(B,C),已知(A)和(B)的循环卷积为(C),求(A).

(n<2^{17})

分析

Bluestein算法的模板题，可以参考这篇博客

再探快速傅里叶变换(FFT)学习笔记(其三)(循环卷积的Bluestein算法+分治FFT+FFT的优化+任意模数NTT)

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#define maxn (1<<17)
const double pi=acos(-1.0);
using namespace std; 
struct com{
	double real;
	double imag;
	com(){
		
	} 
	com(double _real,double _imag){
		real=_real;
		imag=_imag;
	}
	com(double x){
		real=x;
		imag=0;
	}
	void operator = (const com x){
		this->real=x.real;
		this->imag=x.imag;
	}
	void operator = (const double x){
		this->real=x;
		this->imag=0;
	}
	friend com operator + (com p,com q){
		return com(p.real+q.real,p.imag+q.imag);
	}
	friend com operator + (com p,double q){
		return com(p.real+q,p.imag);
	}
	void operator += (com q){
		*this=*this+q;
	}
	void operator += (double q){
		*this=*this+q;
	}
	friend com operator - (com p,com q){
		return com(p.real-q.real,p.imag-q.imag);
	}
	friend com operator - (com p,double q){
		return com(p.real-q,p.imag);
	}
	void operator -= (com q){
		*this=*this-q;
	}
	void operator -= (double q){
		*this=*this-q;
	}
	friend com operator * (com p,com q){
		return com(p.real*q.real-p.imag*q.imag,p.real*q.imag+p.imag*q.real);
	}
	friend com operator * (com p,double q){
		return com(p.real*q,p.imag*q);
	} 
	void operator *= (com q){
		*this=(*this)*q;
	}
	void operator *= (double q){
		*this=(*this)*q;
	}
	friend com operator / (com p,double q){
		return com(p.real/q,p.imag/q);
	} 
	void operator /= (double q){
		*this=(*this)/q;
	} 
	friend com operator / (com p,com q){//复数的除法,类似解二元一次方程,代入复数乘法公式解出答案
		return com((p.real*q.real+p.imag*q.imag)/(q.real*q.real+q.imag*q.imag),(p.imag*q.real-p.real*q.imag)/(q.real*q.real+q.imag*q.imag));
	}
	void print(){
		printf("%lf + %lf i ",real,imag);
	}
};


void fft(com *x,int *rev,int n,int type){
	for(int i=0;i<n;i++) if(i<rev[i]) swap(x[i],x[rev[i]]);
	for(int len=1;len<n;len*=2){
		int sz=len*2;
		com wn1=com(cos(2*pi/sz),type*sin(2*pi/sz));
		for(int l=0;l<n;l+=sz){
			int r=l+len-1;
			com wnk=1;
			for(int i=l;i<=r;i++){
				com tmp=x[i+len];
				x[i+len]=x[i]-wnk*tmp;
				x[i]=x[i]+wnk*tmp;
				wnk=wnk*wn1;
			}
		}
	}
	if(type==-1) for(int i=0;i<n;i++) x[i]/=n;
} 
void bluestein(com *a,int n,int type){ 
	static com x[maxn*4+5],y[maxn*4+5];
	static int rev[maxn*4+5];
	memset(x,0,sizeof(x));
	memset(y,0,sizeof(y));
	int N=1,L=0;
	while(N<n*4){
		L++;
		N*=2;
	}
	for(int i=0;i<N;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
	for(int i=0;i<n;i++) x[i]=com(cos(pi*i*i/n),type*sin(pi*i*i/n))*a[i];
	for(int i=0;i<n*2;i++) y[i]=com(cos(pi*(i-n)*(i-n)/n),-type*sin(pi*(i-n)*(i-n)/n));
	fft(x,rev,N,1);
	fft(y,rev,N,1);
	for(int i=0;i<N;i++) x[i]*=y[i];
	fft(x,rev,N,-1);
	for(int i=0;i<n;i++){
		a[i]=x[i+n]*com(cos(pi*i*i/n),type*sin(pi*i*i/n));
		if(type==-1) a[i]/=n;//一定记得除以n,因为做一次Bluestein相当于一次FFT,IFFT最后要除n,这里也要除n 
	} 
}
void div(com *a,com *b,com *c,int n){//求解A*B=C 
	bluestein(b,n,1);
	bluestein(c,n,1);
	for(int i=0;i<n;i++) a[i]=c[i]/b[i];
	bluestein(a,n,-1);
}

int n;
com a[maxn+5],b[maxn+5],c[maxn+5];
int main(){
	scanf("%d",&n);
	for(int i=0;i<n;i++) scanf("%lf",&b[i].real);
	for(int i=0;i<n;i++) scanf("%lf",&c[i].real);
	div(a,b,c,n);
	for(int i=0;i<n;i++) printf("%.4f
",a[i].real);
}