例题：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("%.4fn",a[i].x);
return 0;
}

最后

以上就是追寻吐司为你收集整理的【模板】Bluestein's algorithm 求循环卷积例题：POJ2821解析：的全部内容，希望文章能够帮你解决【模板】Bluestein's algorithm 求循环卷积例题：POJ2821解析：所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。