bluestein演算法
我們熟知的FFT演算法實際上是將一個多項式在2n個單位根處展開，將其點值對應相乘，並進行逆變換。然而，由於單位根具有「旋轉」的特徵（即$w_{m}^{j}=w_{m}^{j+m}$），若多項式次數大於二分之長度，FFT將進行一次長度為2n的循環卷積。bluestein的演算法是解決了在任意長度上的循環卷積問題。
我們知道，任何一個n次多項式都可以被n+1個點值進行表示，因此如果我們選取所有形如$w_{n+1}^{i}$的單位根並帶入多項式，進行類似於FFT的變化（這裡沒有證明），理應得到正確的結果。
設多項式A為$\sum_{i=0}^{n}{a_i*x^i}$，$F_k$為$A(w_{n+1}^{k})$，則$F_k=\sum_{i=0}^{n}{a_i*w_{n+1}^{ik}}$
考慮ik的另外一種組合含義，即有兩個盒子，每個盒子分別有i個球和k個球，求有多少種拿出兩個球且分別屬於兩個盒子的方法，因此$ik=\tbinom{i+k}{2}-\tbinom{i}{2}-\tbinom{k}{2}$。它的意義在下面推導中可見。
因此$F_k=\sum_{i=0}^{n}{a_i*w_{n+1}^{\tbinom{i+k}{2}-\tbinom{i}{2}-\tbinom{k}{2}}}=w_{n+1}^{-\tbinom{k}{2}}\sum_{i=0}^{n}{a_i*w_{n+1}^{-\tbinom{i}{2}}*w_{n+1}^{\tbinom{i+k}{2}}}$
注意到(i+k)-(i)=k，令$A_{-i}=a_i*w_{n+1}^{-\tbinom{i}{2}}$,$B_i=w_{n+1}^{\tbinom{i}{2}}$。因此，A和B的卷積的第k項即為$F_k$。由於A的下標為負數，我們將A的下標集體加上n。於是，一次bluestein操作花了三次長度為4n的FFT操作。
將多項式轉化為點值表達後，我們依葫蘆畫瓢地將對應位置相乘、進行相應的逆變換（即取單位根的共軛）。而此部分正確性的證明過程是與FFT類似的。
例題：poj2821
1 // 2821 2 #include<cstdio> 3 #include<math.h> 4 #include<cstring> 5 #include<iomanip> 6 #define mod 998244353 7 using namespace std; 8 typedef double ld; 9 const int maxn=(1<<19)+5; 10 const int LIMIT=1<<19; 11 const ld pi=acos(-1); 12 struct com 13 { 14 ld x,y; 15 com(ld a=0,ld b=0):x(a),y(b){} 16 com operator+(const com&A){return com(x+A.x,y+A.y);} 17 com operator-(const com&A){return com(x-A.x,y-A.y);} 18 com operator*(const com&A){return com(x*A.x-y*A.y,x*A.y+y*A.x);} 19 com operator/(const ld&d){return com(x/d,y/d);} 20 com operator/(const com&A){return com(x,y)*com(A.x,-A.y)/(A.x*A.x+A.y*A.y);} 21 void operator/=(const ld&d){x/=d,y/=d;} 22 }; 23 int r[maxn]; 24 inline void DFT(com*A,int limit,int type) 25 { 26 for(int i=1;i<limit;++i) 27 { 28 r[i]=(r[i>>1]>>1)|((i&1)?(limit>>1):0); 29 if(i<r[i]) 30 swap(A[i],A[r[i]]); 31 } 32 for(int len=2;len<=limit;len<<=1) 33 { 34 com w; 35 if(type==1) 36 w=com(cos(pi*2/len),sin(pi*2/len)); 37 else 38 w=com(cos(pi*2/len),-sin(pi*2/len)); 39 for(int i=0;i<limit;i+=len) 40 { 41 com d(1,0); 42 for(int j=0,p1=i,p2=i+len/2;j<len/2;++j,++p1,++p2) 43 { 44 com a=A[p1],b=A[p2]*d; 45 A[p1]=a+b; 46 A[p2]=a-b; 47 d=d*w; 48 } 49 } 50 } 51 } 52 com tmp1[maxn],tmp2[maxn]; 53 54 inline void bluestein(com*A,int n,int type) // n already stands for the number of terms 55 { 56 int limit=1; 57 while(limit<4*n) // 4 times !!!!!!! 58 limit<<=1; 59 for(int i=0;i<limit;++i) 60 tmp1[i]=tmp2[i]=0; 61 for(int i=0;i<n;++i) 62 tmp1[i]=A[i]*com(cos(pi*i*i/n),type*sin(pi*i*i/n)); 63 for(int i=0;i<n*2;++i) 64 tmp2[i]=com(cos(pi*(i-n)*(i-n)/n),-type*sin(pi*(i-n)*(i-n)/n)); 65 DFT(tmp1,limit,1); 66 DFT(tmp2,limit,1); 67 for(int i=0;i<limit;++i) 68 tmp1[i]=tmp1[i]*tmp2[i]; 69 DFT(tmp1,limit,-1); 70 for(int i=0;i<n;++i) 71 A[i]=tmp1[i+n]*com(cos(pi*i*i/n),type*sin(pi*i*i/n))/limit; // dont forget this !!! 72 } 73 com A[maxn],B[maxn],C[maxn]; 74 int n; 75 int main() 76 { 77 scanf("%d",&n); 78 --n; 79 for(int i=0;i<=n;++i) 80 scanf("%lf",&A[i].x); 81 for(int i=0;i<=n;++i) 82 scanf("%lf",&B[i].x); 83 bluestein(A,n+1,1); 84 bluestein(B,n+1,1); 85 for(int i=0;i<n+1;++i) 86 A[i]=B[i]/A[i]; 87 bluestein(A,n+1,-1); 88 for(int i=0;i<=n;++i) 89 A[i].x/=(n+1); 90 for(int i=0;i<=n;++i) 91 printf("%.4f\n",A[i].x); 92 return 0; 93 }
View Code
