淺談FFT和NTT
前言
概述
NOIP2018退役了以後滾回去學文化課了，順便看點新的東西。
FFT和NTT都是用來快速解決多項式乘法的。使用FFT或NTT可以做一些DP優化，以及優化高精度乘法。
FFT適用於所有情況，但是容易丟精度，而且比較慢。
NTT是模意義下的多項式乘法，而且對模數有一定要求，速度比較快。
多項式乘法簡介
一個n次多項式可以表示為A(x)可以一般地表示為：
A(x)=\displaystyle\sum_{i=0}^{n} a_ix^i
同樣，對於另一個n次多項式
B(x)=\displaystyle\sum_{i=0}^{n} b_ix^i
將A(x)和B(x)相加可以得到一個n次多項式，設其為C(x)
則C(x)=\displaystyle\sum_{i=0}^{n} (a_i+b_i)x^i
這就是多項式加法（如果AB的次數不相等，可以在次數少的多項式的高次項上補0）我們可以直接求解，時間複雜度\Theta(n) 。顯然不可能更優。
但是對於多項式乘法，設D(x)=A(x)*B(x) 我們得到的式子是這樣的：
D(x)=\displaystyle\sum_{i=0}^{n+n}(\displaystyle\sum_{j=0}^{n}a_jb_{i-j})x^i
對於這個式子，我們需要
\Theta(n^2)
的時間來求解。
快速傅立葉變換(FFT)
概述
快速傅立葉變換（Fast Fourier Transform）可在$\Theta(nlgn)$的時間內完成多項式乘法。這裡，我們假定n-1為多項式次數，n是2的冪次（實際應用中可通過補0實現）。
點值表達法
之前的多項式表達方式我們稱之為係數表達法，顧名思義就是通過n項的係數來表示多項式。而另一種表達方式稱點值表達，點值表達通過座標軸上該多項式函式影象中的n個點來表示一個多項式。
（證明點值表達的唯一性較為困難，這裡不作證明，但對證明FFT的正確性無影響）
FFT的基本思路是：
A(x)與B(x)的係數表達->A(x)與B(x)的點值表達->
D(x)的點值表達->D(x)的係數表達
如果我們已知A(x)與B(x)的點值表達，那麼我們可以非常輕鬆地用$\Theta(n)$的時間求出D(x)的點值表達。但是對於一般的求值點，我們無法在短時間內實現點值表達與係數表達的快速轉換。
單位複數根
FFT選取的求指點分別是n個n次單位複數根。
【尤拉公式】$e^{i\theta}=isin\theta+cos\theta$
證明：
設$x=isin\theta+cos\theta$
則$\frac{dx}{d\theta}=ix$
$\Rightarrow \frac{1}{x}dx=id\theta$
$\Rightarrow \int \frac{1}{x}dx=\int id\theta$
$\Rightarrow lnx=i\theta$
$\Rightarrow x=e^{i\theta}$
證畢。
設$\omega_n=e^{2\pi i/n}$
n次單位複數根定義為$x^n=1$的x，根據尤拉公式，我們可以找到n個n次單位複數根：$\omega_n^1,\omega_n^2,…,\omega_n^n$
根據尤拉公式，我們可以講這n個根想象成複平面上的n個點，這些點分別與原點連線依次形成n條線，那麼每一條線與下一條線的夾角全部相同。
【推論1】$\omega_{n/2}^{1/2}=\omega_n$
證明：$\omega_{n/2}^{1/2}=e^{i\pi/(n/2)}=e^{2i\pi}=\omega_n$
證畢。
求值
求值就是將係數表達轉化為點值表達。
FFT以這n個單位複數根作為求值點，我們可以進行快速求值。
我們考慮遞迴求解，使用$T(n)=2T(n/2)+O(n)$形式的遞迴式。
可以通過奇偶分類來完成操作：
偶數項分開來是這樣的：
$A_0(x)=a_0+a_2x^2+a_4x^4+…+a_{n-2}x^{n-2}$
同理奇數項是這樣的：
$A_1(x)=a_1x+a_3x^3+a_5x^5+…+a_{n-1}x^{n-1}$
則有$A(x)=A_0(x^2)+xA_1(x^2)$
設結果分別為$y_0,y_1,y_2…,y_{n-1}$
設偶數項結果為$y_{0,0},y_{0,1},y_{0,2},…,y_{0,n/2-1}$
設奇數項結果為$y_{1,0},y_{1,1},y_{1,2},…,y_{1,n/2-1}$
則根據推論1以及單位複數根迴圈的性質，對於任意i<n/2有：
$y_i=y_{0,i}+\omega_n^iy_{1,i}$
同樣，對稱地有：
$y_{i+n/2}=y_{0,i}-\omega_n^iy_{1,i}$
這樣我們就可以得到一個時間複雜度為$\Theta(nlgn)$的演算法。
插值
插值就是將點值表達轉化為係數表達。
對於插值問題，我們可以通過變換模型來把他轉換成一個求值的問題，然後再使用求值中的方法解決。
我們考慮求未知數矩陣的逆矩陣。設未知數矩陣為V，則$V_{i,j}=\omega_n^{ij\%n}=\omega_n^{ij}$
【推論2】$V_{i,j}^{-1}=\omega_n^{ij}/n$。
證明：
設矩陣X滿足：$X_{i,j}=\omega_n^{ij}/n$
設矩陣Y=VX。
則$Y_{i,j}=\displaystyle\sum_{k=0}^{n-1}(\omega_n^{ik}\omega_n^{kj}/n)$
$=\frac{1}{n}\displaystyle\sum_{k=0}^{n-1}\omega_n^{(i-j)k}$
顯然，當i=j時，該式值為1，否則，由於單位複數根的對稱性，該式值為0。
因此$Y=I_n$，則$V^{-1}=X$
證畢。
這樣我們就得到了一個類似於求值問題的式子，只不過所有未知數的係數都反了個頭。我們只需要翻轉一下點值表達的陣列中1~n-1位置的值，再講結果除以n即可。
程式碼
實際實現時只需要用自底向上迭代的方法模擬遞迴即可：
typedef complex<double> C;
const double pi=acos(-1);
void FFT(C* a,int sz){
int n=1<<sz,pos=0;
for(int i=0;i<n-1;++i){
if(i<pos) swap(a[i],a[pos]);
for(int p=sz-1;!(pos&(1<<p+1));--p)
pos^=1<<p;
}
for(int step=1;step<n;step<<=1){
C w(cos(pi/step),sin(pi/step));
for(int i=0;i<n;i+=step<<1){
C p(1,0);
for(int j=i;j<i+step;++j,p*=w){
C x=a[j],y=a[j+step]*p;
a[j]=x+y,a[j+step]=x-y;
}
}
}
}
void times(C* a,C* b,int sz){
FFT(a,sz),FFT(b,sz);
for(int i=0;i<(1<<sz);++i) a[i]*=b[i];
reverse(a+1,a+(1<<sz));
FFT(a,sz);
for(int i=0;i<(1<<sz);++i) a[i]/=1<<sz;
}
NTT
概述
快速數論變換（Number Theoretic Transform）由於解決特殊模數的模意義下的多項式乘法。
原根
對於一些特殊的模數，我們可以使用與單位複數根有相同性質的該模數的原根的若干次方來代替單位複數根。
特殊模數MOD需要滿足性質：(MOD-1)%n=0。設原根為g，那麼$g^((MOD-1)/n)$在模意義下同樣滿足單位複數根在FFT中需要用到的性質。
程式碼
typedef long long ll;
ll fast_pow(ll x,ll n,ll mod){
ll ret=1;
for(int p=1;p<=n;p<<=1,x=x*x%mod)
if(n&p) ret=ret*x%mod;
return ret;
}
void mod(int& x){ x>=MOD?x-=MOD:0; }
void NTT(int* a,int sz){
int n=1<<sz,pos=0;
for(int i=0;i<n-1;++i){
if(i<pos) swap(a[i],a[pos]);
for(int p=sz-1;!(pos&(1<<p+1));--p)
pos^=1<<p;
}
for(int step=1;step<n;step<<=1){
int w=fast_pow(3,(MOD-1)/(step+step),MOD);
for(int i=0;i<n;i+=step<<1){
int p=1;
for(int j=i;j<i+step;++j,p=(ll)p*w%MOD){
int x=a[j],y=(ll)a[j+step]*p%MOD;
mod(a[j]=x+y),mod(a[j+step]=x-y+MOD);
}
}
}
}
void times(int* a,int* b,int* c,int sz){
NTT(a,sz),NTT(b,sz);
for(int i=0;i<(1<<sz);++i) c[i]=(ll)a[i]*b[i]%MOD;
reverse(c+1,c+(1<<sz)),NTT(c,sz);
int inv=fast_pow(1<<sz,MOD-2,MOD);
for(int i=0;i<(1<<sz);++i) c[i]=(ll)c[i]*inv%MOD;
}