FFT是幹什麼的
FFT在演算法競賽中就有一個用途:加速多項式乘法(暴言)
簡單來說,形如 a0X0+a1X1+a2X2+⋯+anXna0X0+a1X1+a2X2+⋯+anXn 的代數表示式叫做多項式,可以記作f(X)=a0X0+a1X1+a2X2+⋯+anXnf(X)=a0X0+a1X1+a2X2+⋯+anXn,其中a0,a1,⋯,ana0,a1,⋯,an叫做多項式的係數,XX是一個不定元(就是不可以合併),不表示任何值,不定元在多項式中最大項的次數稱作多項式的次數
如果我們當前有兩個多項式f(X),g(X)f(X),g(X),現在要把他們乘起來(求卷積),最樸素的做法就是
∑i=02n−1(∑j=0iai∗bi−j)∗xi
這樣的複雜度是Θ(n2)Θ(n2)的,十分不美觀,FFT就是要將這個過程最佳化為Θ(nlogn)Θ(nlogn)
前置技能
多項式複數
複數形如a+bia+bi,其中i=−1−−−√i=−1
aa叫作複數的實部,bibi叫做複數的虛部
複數(a1+b1i)∗(a2+b2i)(a1+b1i)∗(a2+b2i)相乘的值,即a1a2−b1b2+(a1b2+a2b1)ia1a2−b1b2+(a1b2+a2b1)i也是一個複數,同時我們也得到了複數的乘法法則
複數c+dic+di可以用這種方式表示出來
(c+di)
複數乘法的在複平面中表現為輻角相加,模長相乘
單位根
複數ww滿足wn=1wn=1稱作ww是nn次單位根,下圖包含了所有的88次單位根(圖中圓的半徑是1)
8次單位根
同樣的,下圖是所有的4次單位根
4次單位根
聰明的你也許已經發現了單位根的些許性質,即
w2m2n=wmnw2n2m=wnm
wmn=−wm+n2nwnm=−wnm+n2
這兩個要記住,一會很有用
多項式的係數表達法
我們有多項式f(X)=a0X0+a1X1+a2X2+⋯+anXnf(X)=a0X0+a1X1+a2X2+⋯+anXn,令a⃗ =(a0,a1,a2,...,an)a→=(a0,a1,a2,...,an),則稱A(X)A(X)為多項式f(X)f(X)的係數表示法
在係數表示法下,計算多項式乘法是Θ(n2)Θ(n2)的
多項式的點值表達法
任取n+1n+1個互不相同的S={p1,p2,...,pn+1}S={p1,p2,...,pn+1},對f(X)f(X)分別求值得到f(p1),f(p2),...,f(pn+1)f(p1),f(p2),...,f(pn+1),那麼稱A(X)={(p1,f(p1)),(p2,f(p2)),...,(pn+1,f(pn+1))}A(X)={(p1,f(p1)),(p2,f(p2)),...,(pn+1,f(pn+1))}為多項式f(X)f(X)在SS下的點值表示法
可以把多項式想象成一個nn次函式,點值表示法就是取SS下每一個橫座標時對應的點,因為nn次函式可以由n+1n+1個點確定下來(可以將每一個點列一個nn次方程),所以nn維點值與nn維繫數一一對應
更重要的一點,點值表示法下的乘法運算獲得了簡化
兩個多項式PP,QQ分別取點(x,y1)(x,y1)和(x,y2)(x,y2),P∗QP∗Q就會取到點(x,y1∗y2)(x,y1∗y2);令C=P∗QC=P∗Q,因為C(X)=P(X)∗Q(X)C(X)=P(X)∗Q(X),所以C(x)=P(x)∗Q(x)C(x)=P(x)∗Q(x),即C(x)=y1∗y2C(x)=y1∗y2
所以在點值表示法下,計算多項式乘法是Θ(n)Θ(n)的
FFT的具體過程
FFT就是將係數表示法轉化成點值表示法相乘,再由點值表示法轉化為係數表示法的過程,第一個過程叫做求值(DFT),第二個過程叫做插值(IDFT)
求值
還記得我們之前提到的單位根嗎?再回顧一下:
想要求出一個多項式的點值表示法,需要選出n+1n+1個數分別帶入到多項式裡面,帶入一個數的複雜度是Θ(n)Θ(n)的,那麼總複雜度就是Θ(n2)Θ(n2)的,因為單位根有上面兩個優美的性質,所以我們嘗試可以取nn次單位根組成SS,看看能不能加速我們的運算
設A0(X)A0(X)為A(X)A(X)偶次項的和,設A1(X)A1(X)為A(X)A(X)奇次項的和,即
A0(X)=a0x0+a2x1+...+anxn/2A0(X)=a0x0+a2x1+...+anxn/2
A1(X)=a1x0+a3x1+...+an−1xn/2A1(X)=a1x0+a3x1+...+an−1xn/2
因為A(wmn)=a0w0n+a1wmn+a2w2mn+a3w3mn+...+an−1w(n−1)∗mn+anwnmnA(wnm)=a0wn0+a1wnm+a2wn2m+a3wn3m+...+an−1wn(n−1)∗m+anwnnm
所以有
A(wmn)==A0((wmn)2)+wmnA1((wmn)2)A0(wmn2)+wmnA1(wmn2)A(wnm)=A0((wnm)2)+wnmA1((wnm)2)=A0(wn2m)+wnmA1(wn2m)
A(wm+n2n)==A0((wmn)2)+wm+n2nA1((wmn)2)A0(wmn2)−wmnA1(wmn2)A(wnm+n2)=A0((wnm)2)+wnm+n2A1((wnm)2)=A0(wn2m)−wnmA1(wn2m)
也就是說,只要有了A0(X)A0(X)和A1(X)A1(X)的點值表示,就能在Θ(n)Θ(n)時間算出A(X)A(X)的點值表示,對於當前層確定的位置ii,就可以用下一層的兩個值更新當前的值,我們稱這個操作為“蝴蝶變換”
因為這個過程一定要求每層都可以分成兩大小相等的部分,所以多項式最高次項一定是2p2p(p∈Np∈N)次方,如果不是的話,直接在最高次項補零就可以啦
於是我們有了遞迴的寫法:
void FFT(Complex* a,int len){
if(len==1) return;
Complex* a0=new Complex[len/2];
Complex* a1=new Complex[len/2];
for(int i=0;i<len;i+=2){
a0[i/2]=a[i];
a1[i/2]=a[i+1];
}
FFT(a0,len/2);FFT(a1,len/2);
Complex wn(cos(2*Pi/len),sin(2*Pi/len));
Complex w(1,0);
for(int i=0;i<(len/2);i++){
a[i]=a0[i]+w*a1[i];
a[i+len/2]=a0[i]-w*a1[i];
w=w*wn;
return;
}
但遞迴版的FFT常數巨大,實現起來比較複雜,於是又有了迭代的寫法
重新考慮下遞迴FFT的過程,在第ii次求解中,我們將所有元素二進位制ii位為00的放在了左面,ii位為11的放在了右面,事實上,每個元素最終到的是他二進位制顛倒過來的位置
這裡寫圖片描述
再拿一張別人的圖
迭代寫法:
inline void DFT(Complex a[]){
for(int i=0;i<len;i++)///pos[i]代表反轉後的位置
if(i<pos[i])
swap(a[i],a[pos[i]]);
for(int i=2,mid=1;i<=len;i<<=1,mid<<=1){///len代表多項式最高次項
///第一層i是列舉合併到了哪一層。
Complex wm(cos(2.0*pi/i),sin(2.0*pi/i));
for(int j=0;j<len;j+=i){///第二層j是列舉合併區間。
for(int k=j;k<j+mid;k++,w=w*wm){///第三層k是列舉區間內的下標。
Complex l=a[k],r=w*a[k+mid];
a[k]=l+r;a[k+mid]=l-r;
插值
有人證出來插值只要將所有wmnwnm換成wm+n2nwnm+n2,也就是所有的虛部取相反數,再將最終結果除以lenlen,就是插值的過程了
究竟為什麼?我覺得人生一定要有點遺憾可以參考這裡,我就不多說了
支援插值的迭代寫法:
const double DFT=2.0,IDFT=-2.0;///進行求值,第二個引數傳DFT,插值傳IDFT
inline void FFT(Complex a[],double mode){
for(int i=0;i<len;i++)
for(int i=2,mid=1;i<=len;i<<=1,mid<<=1){
Complex wm(cos(2.0*pi/i),sin(mode*pi/i));
for(int j=0;j<len;j+=i){
for(int k=j;k<j+mid;k++,w=w*wm){
if(mode==IDFT)
a[i].x/=len;
現在,你已經會寫FFT了!
FFT是幹什麼的
FFT在演算法競賽中就有一個用途:加速多項式乘法(暴言)
簡單來說,形如 a0X0+a1X1+a2X2+⋯+anXna0X0+a1X1+a2X2+⋯+anXn 的代數表示式叫做多項式,可以記作f(X)=a0X0+a1X1+a2X2+⋯+anXnf(X)=a0X0+a1X1+a2X2+⋯+anXn,其中a0,a1,⋯,ana0,a1,⋯,an叫做多項式的係數,XX是一個不定元(就是不可以合併),不表示任何值,不定元在多項式中最大項的次數稱作多項式的次數
如果我們當前有兩個多項式f(X),g(X)f(X),g(X),現在要把他們乘起來(求卷積),最樸素的做法就是
∑i=02n−1(∑j=0iai∗bi−j)∗xi
∑i=02n−1(∑j=0iai∗bi−j)∗xi
這樣的複雜度是Θ(n2)Θ(n2)的,十分不美觀,FFT就是要將這個過程最佳化為Θ(nlogn)Θ(nlogn)
前置技能
多項式複數
複數形如a+bia+bi,其中i=−1−−−√i=−1
aa叫作複數的實部,bibi叫做複數的虛部
複數(a1+b1i)∗(a2+b2i)(a1+b1i)∗(a2+b2i)相乘的值,即a1a2−b1b2+(a1b2+a2b1)ia1a2−b1b2+(a1b2+a2b1)i也是一個複數,同時我們也得到了複數的乘法法則
複數c+dic+di可以用這種方式表示出來
(c+di)
複數乘法的在複平面中表現為輻角相加,模長相乘
單位根
複數ww滿足wn=1wn=1稱作ww是nn次單位根,下圖包含了所有的88次單位根(圖中圓的半徑是1)
8次單位根
同樣的,下圖是所有的4次單位根
4次單位根
聰明的你也許已經發現了單位根的些許性質,即
w2m2n=wmnw2n2m=wnm
wmn=−wm+n2nwnm=−wnm+n2
這兩個要記住,一會很有用
多項式的係數表達法
我們有多項式f(X)=a0X0+a1X1+a2X2+⋯+anXnf(X)=a0X0+a1X1+a2X2+⋯+anXn,令a⃗ =(a0,a1,a2,...,an)a→=(a0,a1,a2,...,an),則稱A(X)A(X)為多項式f(X)f(X)的係數表示法
在係數表示法下,計算多項式乘法是Θ(n2)Θ(n2)的
多項式的點值表達法
任取n+1n+1個互不相同的S={p1,p2,...,pn+1}S={p1,p2,...,pn+1},對f(X)f(X)分別求值得到f(p1),f(p2),...,f(pn+1)f(p1),f(p2),...,f(pn+1),那麼稱A(X)={(p1,f(p1)),(p2,f(p2)),...,(pn+1,f(pn+1))}A(X)={(p1,f(p1)),(p2,f(p2)),...,(pn+1,f(pn+1))}為多項式f(X)f(X)在SS下的點值表示法
可以把多項式想象成一個nn次函式,點值表示法就是取SS下每一個橫座標時對應的點,因為nn次函式可以由n+1n+1個點確定下來(可以將每一個點列一個nn次方程),所以nn維點值與nn維繫數一一對應
更重要的一點,點值表示法下的乘法運算獲得了簡化
兩個多項式PP,QQ分別取點(x,y1)(x,y1)和(x,y2)(x,y2),P∗QP∗Q就會取到點(x,y1∗y2)(x,y1∗y2);令C=P∗QC=P∗Q,因為C(X)=P(X)∗Q(X)C(X)=P(X)∗Q(X),所以C(x)=P(x)∗Q(x)C(x)=P(x)∗Q(x),即C(x)=y1∗y2C(x)=y1∗y2
所以在點值表示法下,計算多項式乘法是Θ(n)Θ(n)的
FFT的具體過程
FFT就是將係數表示法轉化成點值表示法相乘,再由點值表示法轉化為係數表示法的過程,第一個過程叫做求值(DFT),第二個過程叫做插值(IDFT)
求值
還記得我們之前提到的單位根嗎?再回顧一下:
w2m2n=wmnw2n2m=wnm
wmn=−wm+n2nwnm=−wnm+n2
想要求出一個多項式的點值表示法,需要選出n+1n+1個數分別帶入到多項式裡面,帶入一個數的複雜度是Θ(n)Θ(n)的,那麼總複雜度就是Θ(n2)Θ(n2)的,因為單位根有上面兩個優美的性質,所以我們嘗試可以取nn次單位根組成SS,看看能不能加速我們的運算
設A0(X)A0(X)為A(X)A(X)偶次項的和,設A1(X)A1(X)為A(X)A(X)奇次項的和,即
A0(X)=a0x0+a2x1+...+anxn/2A0(X)=a0x0+a2x1+...+anxn/2
A1(X)=a1x0+a3x1+...+an−1xn/2A1(X)=a1x0+a3x1+...+an−1xn/2
因為A(wmn)=a0w0n+a1wmn+a2w2mn+a3w3mn+...+an−1w(n−1)∗mn+anwnmnA(wnm)=a0wn0+a1wnm+a2wn2m+a3wn3m+...+an−1wn(n−1)∗m+anwnnm
所以有
A(wmn)==A0((wmn)2)+wmnA1((wmn)2)A0(wmn2)+wmnA1(wmn2)A(wnm)=A0((wnm)2)+wnmA1((wnm)2)=A0(wn2m)+wnmA1(wn2m)
A(wm+n2n)==A0((wmn)2)+wm+n2nA1((wmn)2)A0(wmn2)−wmnA1(wmn2)A(wnm+n2)=A0((wnm)2)+wnm+n2A1((wnm)2)=A0(wn2m)−wnmA1(wn2m)
也就是說,只要有了A0(X)A0(X)和A1(X)A1(X)的點值表示,就能在Θ(n)Θ(n)時間算出A(X)A(X)的點值表示,對於當前層確定的位置ii,就可以用下一層的兩個值更新當前的值,我們稱這個操作為“蝴蝶變換”
因為這個過程一定要求每層都可以分成兩大小相等的部分,所以多項式最高次項一定是2p2p(p∈Np∈N)次方,如果不是的話,直接在最高次項補零就可以啦
於是我們有了遞迴的寫法:
void FFT(Complex* a,int len){
if(len==1) return;
Complex* a0=new Complex[len/2];
Complex* a1=new Complex[len/2];
for(int i=0;i<len;i+=2){
a0[i/2]=a[i];
a1[i/2]=a[i+1];
}
FFT(a0,len/2);FFT(a1,len/2);
Complex wn(cos(2*Pi/len),sin(2*Pi/len));
Complex w(1,0);
for(int i=0;i<(len/2);i++){
a[i]=a0[i]+w*a1[i];
a[i+len/2]=a0[i]-w*a1[i];
w=w*wn;
}
return;
}
但遞迴版的FFT常數巨大,實現起來比較複雜,於是又有了迭代的寫法
重新考慮下遞迴FFT的過程,在第ii次求解中,我們將所有元素二進位制ii位為00的放在了左面,ii位為11的放在了右面,事實上,每個元素最終到的是他二進位制顛倒過來的位置
這裡寫圖片描述
再拿一張別人的圖
這裡寫圖片描述
迭代寫法:
inline void DFT(Complex a[]){
for(int i=0;i<len;i++)///pos[i]代表反轉後的位置
if(i<pos[i])
swap(a[i],a[pos[i]]);
for(int i=2,mid=1;i<=len;i<<=1,mid<<=1){///len代表多項式最高次項
///第一層i是列舉合併到了哪一層。
Complex wm(cos(2.0*pi/i),sin(2.0*pi/i));
for(int j=0;j<len;j+=i){///第二層j是列舉合併區間。
Complex w(1,0);
for(int k=j;k<j+mid;k++,w=w*wm){///第三層k是列舉區間內的下標。
Complex l=a[k],r=w*a[k+mid];
a[k]=l+r;a[k+mid]=l-r;
}
}
}
return;
}
插值
有人證出來插值只要將所有wmnwnm換成wm+n2nwnm+n2,也就是所有的虛部取相反數,再將最終結果除以lenlen,就是插值的過程了
究竟為什麼?我覺得人生一定要有點遺憾可以參考這裡,我就不多說了
支援插值的迭代寫法:
const double DFT=2.0,IDFT=-2.0;///進行求值,第二個引數傳DFT,插值傳IDFT
inline void FFT(Complex a[],double mode){
for(int i=0;i<len;i++)
if(i<pos[i])
swap(a[i],a[pos[i]]);
for(int i=2,mid=1;i<=len;i<<=1,mid<<=1){
Complex wm(cos(2.0*pi/i),sin(mode*pi/i));
for(int j=0;j<len;j+=i){
Complex w(1,0);
for(int k=j;k<j+mid;k++,w=w*wm){
Complex l=a[k],r=w*a[k+mid];
a[k]=l+r;a[k+mid]=l-r;
}
}
}
if(mode==IDFT)
for(int i=0;i<len;i++)
a[i].x/=len;
return;
}
現在,你已經會寫FFT了!