Bluestein's Algorithm學習筆記
阿新 • • 發佈:2021-07-31
說是學習筆記,其實也沒什麼可寫的
直接上式子:
\[\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; }
浮世蒼茫,不過瞬逝幻夢
善惡愛誑,皆有定數
於命運之輪中
吞噬於黃泉之冥暗
嗚呼,吾乃夢之戍人
幻戀之觀者
唯於萬華鏡中,永世長存