Bluestein's Algorithm学习笔记

说是学习笔记,其实也没什么可写的

直接上式子:

\[\begin{aligned} \hat{a_i}&=\sum_{j=0}^n\omega_n^{ij}a_j\\ &=\sum_{j=0}^n\omega_{2n}^{2ij}a_j\\ &=\sum_{j=0}^n\omega_{2n}^{-(i-j)^2+i^2+j^2}a_j\\ &=\omega_{2n}^{i^2}\sum_{j=0}^n(\omega_{2n}^{j^2}a_j)(\omega_{2n}^{-(i-j)^2}) \end{aligned} \]

这就已经是卷积的形式了,\(\tt FFT/NTT\)即可

由于\(-n\leq i-j\leq n\),所以我们需要向右平移\(n\)位,最后再还原

void FFT(Complex *a,int len,int Ty){
	int *r=new int[len],L=log2(len);r[0]=0;
	for(int i=0;i<len;++i){
		r[i]=(r[i>>1]>>1)|((i&1)<<L-1);
		if(i<r[i])swap(a[i],a[r[i]]);
	}//puts("OK");
	for(int i=1;i<len;i<<=1){
		Complex T(cos(Pi/i),Ty*sin(Pi/i));
		for(int W=i<<1,j=0;j<len;j+=W){
			Complex	t(1,0);
			for(int k=0;k<i;++k,t*=T){
				Complex x(a[j+k]),y(t*a[i+j+k]);
				a[j+k]=x+y;
				a[i+j+k]=x-y;
			}
		}
	}delete[] r;
	if(Ty==-1)
		for(int i=0;i<len;++i)
			a[i]/=len;
}
void BS(Complex *a,int len,int Ty){
	static Complex b[N],c[N];
	for(int i=0;i<len;++i)b[i]=a[i]*Complex(cos(Pi*i*i/len),Ty*sin(Pi*i*i/len));
	for(int i=0;i<len*2;++i)c[i]=Complex(cos(Pi*(len-i)*(len-i)/len),-Ty*sin(Pi*(len-i)*(len-i)/len));
	int s=1;while(s<(len*3))s<<=1;
	for(int i=len;i<s;++i)b[i]=Complex(0,0);
	for(int i=len*2;i<s;++i)c[i]=Complex(0,0);
	FFT(b,s,1);FFT(c,s,1);
	for(int i=0;i<s;++i)b[i]*=c[i];
	FFT(b,s,-1);
	for(int i=0;i<len;++i)a[i]=b[i+len]*Complex(cos(Pi*i*i/len),Ty*sin(Pi*i*i/len));
	if(Ty==-1)for(int i=0;i<len;++i)a[i]/=len;
}
上一篇:Google资源管理器简析


下一篇:模拟赛2021.4.5