再探快速傅里葉變換(FFT)學習筆記(其三)(循環卷積的Bluestein算法+分治FFT+FFT的優化+任意模數NTT)
寫在前面
為了不使篇幅過長，預計將把學習筆記分為四部分:
- DFT,IDFT,FFT的定義,實現與證明:快速傅里葉變換(FFT)學習筆記(其一)
- NTT的實現與證明:快速傅里葉變換(FFT)學習筆記(其二)
- 任意模數NTT與FFT的優化技巧
- 多項式相關操作
一些約定
- \([p(x)]=\begin{cases}1,p(x)為真 \\ 0,p(x)為假 \end{cases}\)
- 本文中序列的下標從0開始
- 若\(s\)是一個序列,\(|s|\)表示\(s\)的長度
- 若大寫字母如\(F(x)\)表示一個多項式,那么對應的小寫字母如\(f\)表示多項式的每一項系數，即\(F(x)=\sum_{i=0}^{n-1} f_ix^i\)
循環卷積
DFT卷積的本質
考慮在(其一)中提到的卷積的定義式。
我們一般做FFT時忽略了式子中的\(\bmod\),其實它是在\(\bmod 2^q\)的意義下的循環卷積,只是因為\(|a|,|b|,|c|<2^q\),所以取不取模都沒什么影響。
如果序列長度\(n\)是2的整數次冪,那么直接做就可以了。
如果序列長度\(n\)不是2的整數次冪考慮暴力的做法:先做一次普通FFT,再把\(c_{k+n}\)加到\(c_k\)上。但是這樣顯然不是很優秀。下面給出了一種在\(O(n \log n)\)的時間內實現任意長度循環卷積的算法:Bluestein’s Algorithm
Bluestein’s Algorithm
注:原論文的推導可能有誤
考慮DFT的式子
不妨設
$y_j=\omega_n^{-\frac{j^2}{2}}= \cos \frac{\pi j^2}{n}-\text{i}\sin \frac{\pi j^2}{n} $
那么\(a_i'=\omega_n^{\frac{j^2}{2}}\sum_{j=0}^{n-1} x_j y_{i-j}\)
這已經很類似卷積的形式了,但是注意到\(j\)的上界是\(n-1\)而不是\(i\),\(j-i\)可能為負數。那么我們把\(y\)數組的長度擴大到\(2n\),定義:
$y_j=\omega_n^{-\frac{(j-n)^2}{2}}= \cos \frac{\pi (j-n)^2}{n}-\text{i}\sin \frac{\pi (j-n)^2}{n} $.
這樣\(j<n\)的時候就對應了\(j-i\)為負數的情形,\(j\geq n\)就對應了\(j-i\)為正的情形。然后對\(x\)和\(y\)用一般的FFT,最后的答案存儲在\(i+n\)的位置上，也就是說真正的\(a'_i\)實際上對應了乘積結果的\((x \cdot y)_{i+n}\)
這樣，我們就只做了一次FFT就求出了任意長度循環卷積。逆變換同理,只是換成共軛復數。注意到在上述的推導中我們沒有用到單位根\(\omega\)的任何性質,因此這里的\(\omega\)可以換成任意復數\(z\),這樣的變換稱為Chirp Z-Transform,CZT.可見,CZT實際上是DFT的廣義形式。
代碼實現見下方模板題
例題
這是Bluestein算法的模板題
[POJ 2821] 給出兩個長度為\(n\)的序列\(B,C\),已知\(A\)和\(B\)的循環卷積為\(C\),求\(A\).
\(n<2^{17}\)
代碼:
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#define maxn (1<<17)
const double pi=acos(-1.0);
using namespace std;
struct com{
double real;
double imag;
com(){
}
com(double _real,double _imag){
real=_real;
imag=_imag;
}
com(double x){
real=x;
imag=0;
}
void operator = (const com x){
this->real=x.real;
this->imag=x.imag;
}
void operator = (const double x){
this->real=x;
this->imag=0;
}
friend com operator + (com p,com q){
return com(p.real+q.real,p.imag+q.imag);
}
friend com operator + (com p,double q){
return com(p.real+q,p.imag);
}
void operator += (com q){
*this=*this+q;
}
void operator += (double q){
*this=*this+q;
}
friend com operator - (com p,com q){
return com(p.real-q.real,p.imag-q.imag);
}
friend com operator - (com p,double q){
return com(p.real-q,p.imag);
}
void operator -= (com q){
*this=*this-q;
}
void operator -= (double q){
*this=*this-q;
}
friend com operator * (com p,com q){
return com(p.real*q.real-p.imag*q.imag,p.real*q.imag+p.imag*q.real);
}
friend com operator * (com p,double q){
return com(p.real*q,p.imag*q);
}
void operator *= (com q){
*this=(*this)*q;
}
void operator *= (double q){
*this=(*this)*q;
}
friend com operator / (com p,double q){
return com(p.real/q,p.imag/q);
}
void operator /= (double q){
*this=(*this)/q;
}
friend com operator / (com p,com q){//復數的除法,類似解二元一次方程,代入復數乘法公式解出答案
return com((p.real*q.real+p.imag*q.imag)/(q.real*q.real+q.imag*q.imag),(p.imag*q.real-p.real*q.imag)/(q.real*q.real+q.imag*q.imag));
}
void print(){
printf("%lf + %lf i ",real,imag);
}
};
void fft(com *x,int *rev,int n,int type){
for(int i=0;i<n;i++) if(i<rev[i]) swap(x[i],x[rev[i]]);
for(int len=1;len<n;len*=2){
int sz=len*2;
com wn1=com(cos(2*pi/sz),type*sin(2*pi/sz));
for(int l=0;l<n;l+=sz){
int r=l+len-1;
com wnk=1;
for(int i=l;i<=r;i++){
com tmp=x[i+len];
x[i+len]=x[i]-wnk*tmp;
x[i]=x[i]+wnk*tmp;
wnk=wnk*wn1;
}
}
}
if(type==-1) for(int i=0;i<n;i++) x[i]/=n;
}
void bluestein(com *a,int n,int type){
static com x[maxn*4+5],y[maxn*4+5];
static int rev[maxn*4+5];
memset(x,0,sizeof(x));
memset(y,0,sizeof(y));
int N=1,L=0;
while(N<n*4){
L++;
N*=2;
}
for(int i=0;i<N;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
for(int i=0;i<n;i++) x[i]=com(cos(pi*i*i/n),type*sin(pi*i*i/n))*a[i];
for(int i=0;i<n*2;i++) y[i]=com(cos(pi*(i-n)*(i-n)/n),-type*sin(pi*(i-n)*(i-n)/n));
fft(x,rev,N,1);
fft(y,rev,N,1);
for(int i=0;i<N;i++) x[i]*=y[i];
fft(x,rev,N,-1);
for(int i=0;i<n;i++){
a[i]=x[i+n]*com(cos(pi*i*i/n),type*sin(pi*i*i/n));
if(type==-1) a[i]/=n;//一定記得除以n,因為做一次Bluestein相當於一次FFT,IFFT最后要除n,這里也要除n
}
}
void div(com *a,com *b,com *c,int n){//求解A*B=C
bluestein(b,n,1);
bluestein(c,n,1);
for(int i=0;i<n;i++) a[i]=c[i]/b[i];
bluestein(a,n,-1);
}
int n;
com a[maxn+5],b[maxn+5],c[maxn+5];
int main(){
scanf("%d",&n);
for(int i=0;i<n;i++) scanf("%lf",&b[i].real);
for(int i=0;i<n;i++) scanf("%lf",&c[i].real);
div(a,b,c,n);
for(int i=0;i<n;i++) printf("%.4f\n",a[i].real);
}
分治FFT
//填坑中
FFT的弱常數優化
下面介紹一些優化FFT的常數的技巧。雖然這些技巧都只是對FFT的一些小優化，但是在某些題目中優化效果極其明顯。
復雜算式中減少FFT次數
如果我們要計算一個復雜的多項式,如\(A(x)=B(x)C(x)+D(x)E(x)\)
最簡單的方法是分別計算\(B(x)C(x)\)和\(D(x)E(x)\),這樣需要做6次FFT. 但是如果先對\(B,C,D,E\)做DFT,然后直接用點值表達式計算\(a_i=b_ic_i+d_ie_i\),再把\(a\)IDFT回去。這樣只需要做5次FFT,且多項式越復雜，這樣的常數就越優秀。
例題
[BZOJ 3771] Triple(FFT+容斥原理+生成函數)
利用循環卷積
考慮對於兩個長度為\(n\)的序列\(a,b\),計算它們的卷積\(c\)的第\(0.5n\)項到第\(1.5n\)項。傳統的方法是補0擴充到\(2n\)的序列。但是因為FFT求得實際上是我們已經提到過的循環卷積,所以如果只補0到\(1.5n\)(上取整),對第\(0.5n\)項到第\(1.5n\)項無影響
在基於牛頓迭代的算法中，能起到較明顯的優化作用。會在(其四)中詳細介紹這些算法。
小范圍暴力
由於FFT的常數較大。在數據范圍較小的時候甚至不如\(O(n^2)\)的暴力卷積的優秀。因此在做多次FFT和分治FFT的時候,如果當前的序列長度較小，可以采用暴力算法。
例題
[BZOJ 3509] [CodeChef] COUNTARI (FFT+分塊)
快速冪乘法次數的優化
這個東西實際上比較雞肋。因為多項式快速冪可以通過多項式\(\ln\)和\(\exp\)優化到\(O(n \log n)\).但是為了應對考場上時間不夠的情況,我們來考慮如何通過簡單的實現來減少\(O(n \log^2n)\)的倍增快速冪的復雜度。
倍增法的思路是根據前面算過的乘積快速算出當前的乘積,如\(1 \to 2 \to 4 \to 8\).最壞情況下需要\(2 \log_2n+C\)次乘法。但這並不是下界。我們定義additional chain為一條鏈,最開始是1,后一個數減前一個數的差是鏈上這個是前面的某一個數。例如\(1 \to 2 \to 4 \to 6\).\(6-4=2\)在前面出現過,\(4-2=2\)在前面出現過。那么根據這條additional chain計算6次冪的時候，可以從1次冪出發,用1次冪乘1次冪得到2次冪,再乘2次冪得到4次冪，再乘2次冪得到6次冪。
很可惜，對於數\(k\)求出得到\(k\)的最短additional chain是NP-hard的。但是有很好的近似算法。近似算法基於BFS。每次我們對於隊頭的數\(x\),枚舉它對應的additional chain中的數\(y\),如果\(x+y\)還沒有訪問過那么將其入隊,並將\(x\)對應的鏈后面接上\(x+y\). 這個預處理是\(O(k)\)的,且對快速冪的常數優化很顯著。
如果\(k\)很大,比如\(10^{10000}\),可以采用十進制快速冪。但是用Method of Four Russians(俗稱四毛子算法),可以將乘法次數減少到\(\log_2n+O(\frac{\log n}{\log \log n})\).具體方法見2017年國家集訓隊論文《非常規大小分塊算法初探》