爱站程序员基地我们熟知的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项即为$\\frac{F_k}{w_{n+1}^{-\\tbinom{k}{2}}}$。由于A的下标为负数，我们将A的下标集体加上n。于是，一次bluestein操作花了三次长度为4n的FFT操作。
将多项式转化为点值表达后，我们依葫芦画瓢地将对应位置相乘、进行相应的逆变换（即取单位根的共轭）。而此部分正确性的证明过程是与FFT类似的。
例题：poj2821。对于循环卷积B=A*C，求C。
1 #include<cstdio>2 #include<math.h>3 #include<cstring>4 #include<iomanip>5 using namespace std;6 typedef long long int ll;7 typedef double ld;8 const int maxn=(1<<19)+5;9 const ld pi=acos(-1);10 struct com11 {12 ld x,y;13 com(ld a=0,ld b=0):x(a),y(b){}14 com operator+(const com&A){return com(x+A.x,y+A.y);}15 com operator-(const com&A){return com(x-A.x,y-A.y);}16 com operator*(const com&A){return com(x*A.x-y*A.y,x*A.y+y*A.x);}17 com operator/(const ld&d){return com(x/d,y/d);}18 com operator/(const com&A){return com(x,y)*com(A.x,-A.y)/(A.x*A.x+A.y*A.y);}19 }A[maxn],B[maxn];20 int r[maxn];21 inline void DFT(com*A,int limit,int type)22 {23 for(int i=1;i<limit;++i)24 {25 r[i]=(r[i>>1]>>1)|((i&1)?(limit>>1):0);26 if(i<r[i])27 swap(A[i],A[r[i]]);28 }29 for(int len=2;len<=limit;len<<=1)30 {31 com d(cos(pi*2/len),sin(pi*2/len)*type);32 for(int i=0;i<limit;i+=len)33 {34ad8com w(1,0);35 for(int j=0,p1=i,p2=i+len/2;j<len/2;++j,++p1,++p2)36 {37 com x=A[p1],y=A[p2]*w;38 A[p1]=x+y;39 A[p2]=x-y;40 w=w*d;41 }42 }43 }44 }45 com tmp1[maxn],tmp2[maxn];46 inline com get(int k,int n,int type)47 {48 ll g=(ll)k*(k-1)/2;49 return com(cos(pi*2/(n+1)*g),sin(pi*2/(n+1)*g)*type);50 }51 inline void bluestein(com*A,int n,int type)52 {53 int limit=1;54 while(limit<4*n+1)55 limit<<=1;56 for(int i=0;i<limit;++i)57 tmp1[i]=tmp2[i]=com(0,0);58 for(int i=0;i<=n;++i)59 tmp1[n-i]=A[i]*get(i,n,-type);60 for(int i=0;i<=n*2;++i)61 tmp2[i]=get(i,n,type);62 DFT(tmp1,limit,1);63 DFT(tmp2,limit,1);64 for(int i=0;i<limit;++i)65 tmp1[i]=tmp1[i]*tmp2[i];66 DFT(tmp1,limit,-1);67 for(int i=0;i<=n;++i)68 A[i]=tmp1[i+n]/limit*get(i,n,-type);69 }70 int n;71 int main()72 {73 scanf(\"%d\",&n);74 --n;75 for(int i=0;i<=n;++i)76 scanf(\"%lf\",&A[i].x);77 for(int i=0;i<=n;++i)78 scanf(\"%lf\",&B[i].x);79 bluestein(A,n,1);80 bluestein(B,n,1);81 for(int i=0;i<=n;++i)82 A[i]=B[i]/A[i];83 bluestein(A,n,-1);84 for(int i=0;i<=n;++i)85 A[i].x/=(n+1);86 for(int i=0;i<=n;++i)87 printf(\"%.4f\\n\",A[i].x);88 return 0;89 }
View Code
