【模板】Bluestein's algorithm 求循环卷积

例题：POJ2821

---

解析：

首先我们知道一般的，最常用的FFT求的就是在$\%2^n$意义下的循环卷积。

换句话说，长度为$n$的$DFT$求的就是在长度$\%n$下的循环卷积。

现在考虑长度不为$2$的整数次幂的时候我们怎么办。

令长度为$n$，还是参考$DFT$的式子：$$A_k=\sum _{i=0}^{n-1}{a}_i{\omega }_n^{ik}$$

由于$2ik=i^2+k^2-(i-k{)}^2$，于是我们可以把单位根换成$2n$次再来考虑。

于是$$A_k={\omega }_{2n}^{k^2}\sum _{i=0}^{n-1}{a}_i{\omega }_{2n}^{i^2}\cdot {\omega }_{2n}^{-(i-k{)}^2}$$

发现$\Sigma$后面是一个卷积，于是我们考虑用卷积实现$DFT$就行了。

构造两个序列，长度分别为$n$和$2n$，具体实现可以参考代码里面。

---

代码：

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
#define re register
#define cs const

cs double PI=acos(-1);
struct Complex{
	double x,y;
	Complex(){}
	Complex(cs double &_x,cs double &_y):x(_x),y(_y){}
	Complex conj()cs{return Complex(x,-y);}
	friend Complex operator+(cs Complex &a,cs Complex &b){return Complex(a.x+b.x,a.y+b.y);}
	friend Complex operator-(cs Complex &a,cs Complex &b){return Complex(a.x-b.x,a.y-b.y);}
	friend Complex operator*(cs Complex &a,cs Complex &b){return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
	friend Complex operator/(cs Complex &a,cs double &b){return Complex(a.x/b,a.y/b);}
	friend Complex operator/(cs Complex &a,cs Complex &b){return a*b.conj()/(b.x*b.x+b.y*b.y);}
};

cs int N=1<<17|1;

int r[N<<2|1];
Complex w[N<<2|1];
inline void FFT(Complex *A,int len,int typ){
	if(typ==-1)std::reverse(A+1,A+len);
	for(int re i=0;i<len;++i)if(i<r[i])std::swap(A[i],A[r[i]]);
	for(int re i=1;i<len;i<<=1)
	for(int re j=0;j<len;j+=i<<1)
	for(int re k=0;k<i;++k){
		static Complex x,y;
		x=A[j+k],y=A[j+k+i]*w[len/i/2*k];
		A[j+k]=x+y;
		A[j+k+i]=x-y;
	}
	if(typ==-1)for(int re i=0;i<len;++i)A[i]=A[i]/len;
}
inline void init_rev(int len){
	static int last;if(last==len)return ;last=len;
	for(int re i=0;i<len;++i)
	r[i]=r[i>>1]>>1|((i&1)*(len>>1)),
	w[i]=Complex(cos(2*PI*i/len),sin(2*PI*i/len));
}

inline void BS(Complex *a,int l,int typ){
	static Complex A[N<<2|1],B[N<<2|1];
	int l2=l<<1;
	for(int re i=0;i<l;++i)A[i]=a[i]*Complex(cos(PI*i*i/l),typ*sin(PI*i*i/l));
	for(int re i=0;i<l2;++i)B[i]=Complex(cos(PI*(i-l)*(i-l)/l),-typ*sin(PI*(i-l)*(i-l)/l));
	int len=1;while(len<l*3)len<<=1;
	for(int re i=l;i<len;++i)A[i]=Complex(0,0);
	for(int re i=l2;i<len;++i)B[i]=Complex(0,0);
	init_rev(len);
	FFT(A,len,1),FFT(B,len,1);
	for(int re i=0;i<len;++i)A[i]=A[i]*B[i];
	FFT(A,len,-1);
	for(int re i=0;i<l;++i)a[i]=A[i+l]*Complex(cos(PI*i*i/l),typ*sin(PI*i*i/l));
	if(typ==-1)for(int re i=0;i<l;++i)a[i]=a[i]/l;
}

int n;
Complex a[N],b[N],c[N];

signed main(){
	scanf("%d",&n);
	for(int re i=0;i<n;++i)scanf("%lf",&b[i].x);
	for(int re i=0;i<n;++i)scanf("%lf",&c[i].x);
	BS(b,n,1),BS(c,n,1);
	for(int re i=0;i<n;++i)a[i]=c[i]/b[i];
	BS(a,n,-1);
	for(int re i=0;i<n;++i)printf("%.4f\n",a[i].x);
	return 0;
}