bluestein算法

咱們熟知的FFT算法其實是將一個多項式在2ⁿ個單位根處展開，將其點值對應相乘，並進行逆變換。然而，因爲單位根具備「旋轉」的特徵（即$w_{m}^{j}=w_{m}^{j+m}$），若多項式次數大於二分之長度，FFT將進行一次長度爲2ⁿ的循環卷積。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 com
11 {
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         {
34             com 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