ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

Bluestein's Algorithm学习笔记

2021-07-30 21:03:42  阅读:179  来源: 互联网

标签:Algorithm int sum 笔记 leq len 2n omega Bluestein


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

直接上式子:

\[\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;
}

标签:Algorithm,int,sum,笔记,leq,len,2n,omega,Bluestein
来源: https://www.cnblogs.com/SYDevil/p/15081285.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有