FFT&NTT(以及擴展)
預備知識：用於NTT
NTT/FFT其實本質相同，用途是快速求解 多項式乘積
前言
FT: 傅里葉變換:
這是一個工程上的概念，可以簡述為：一個周期性的信號波段可以用 若干個正弦曲線 的帶權和表示
DFT: 離散傅里葉變換，這是傅里葉變換在離散情況下的變種
FFT: 快速傅里葉變換
NTT: 快速數論變換
談及核心思想
1.單位根:
構造\(\omega_n\)為\(n\)階單位根(不知道\(\omega_n\)的值域)，滿足性質\(\omega_n^n=\omega_n^0=1\)
對於\(2|n\),\(\omega _n^{\frac{n}{2}}=-1\)
顯然\(\omega_n\)滿足一個非常簡單的性質：折半引理 \(\begin{aligned} \forall 2|i\and 2|n , \omega_n^i=\omega_{\frac{n}{2}}^{\frac{i}{2}}\end{aligned}\)
\(\omega_n\)實際上是一個在冪次上呈現\(n\)元循環的數值
2.多項式與點值式的轉化
一個\(n\)階多項式最普通的表示就是\(F(x)=\sum_{i=0}^{n-1} a_ix^i\)
然而，多項式也可以用\(n\)個互不相關的點表示，即\((x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1})\)
兩者可以互相轉化
對於同\(x_i\)的點值，兩個多項式卷積時，其\(y_i\)可以直接對應相乘
FFT/NTT的核心過程是
多項式\(\longrightarrow\) 點值式\(\longrightarrow\)點值式對應相乘\(\longrightarrow\)多項式
而用單位根來構造快速的多項式與點值式的轉化
3.分治思想
用於降低多項式與點值式轉換的復雜度
FFT的單位根
\((x,y)\)指復數\(i=\sqrt{-1},(x,y)=x+yi\)
基本運算\((x,y)+(a,b)=(x+a,y+b),(x,y)\cdot (a,b)=(ax-by,ay+bx)\)
FFT的單位根是:\(\omega_n\)=\((cos(\frac{2\pi}{n}),sin(\frac{2\pi}{n}))\)
而\(\omega _n^i=(cos(\frac{2\pi}{n}\cdot i),sin(\frac{2\pi}{n}\cdot i))\) (展開發現就是三角函數求和公式)
顯然滿足單位根的性質
(實際上可以發現，這個說是點值其實就是信號序列的三角函數表示)
NTT
相信您已經了解了原根的一些性質，\(\text{NTT}\)的單位根常用原根構造
\(\text{NTT}\)的單位根實際有較大的局限性，對於質數\(P\)只能構造出\(n|P-1,\omega_n=g^{\frac{P-1}{n}}\)
計算在模意義下就能滿足單位根的性質
通常我們\(P\)取\(998244353\)，\(2^{23}|(P-1)\)，它的一個原根是3
實際上，為了滿足下面分治需要，構造的模數通常滿足\(P-1=s\cdot 2^t\)的\(t\)較大，這類模數我們常稱作\(\text{NTT}\)模數
以上部分均為基礎知識，相對來說應該不會太難，下面是主要難點
多項式轉點值式
接下來我們考慮如何將多項式轉化為點值式
對於點值式，我們構造的點橫坐標為\(x_i=\omega_n^i\)
具體目標是對於函數\(F(x)\)，求出在\(x_0,x_1,\cdots ,x_{n-1}\)上的函數值
即求出\(F(x_i)=a_0\omega_n^0+a_1\omega_n^{i}+a_2\omega_n^{2i}+\cdots\)
接下來就是核心的分治思想，注意，這里的分治是子問題嚴格等大的
對於當前問題，分成兩部分子問題求解(實際是可以分成多部分的，但是這個是特殊情況暫時不予討論)，即求解
令\(m=\frac{n}{2}\)
\(\begin{aligned} 2|i,G(x_i)=a_0\omega_{m}^0+a_2\omega_{m}^{\frac{i}{2}}+a_4\omega_{m}^{\frac{i}{2}\cdot 2}+\cdots\end{aligned}\)
\(\begin{aligned} 2|i,H(x_i)=a_1\omega_{m}^0+a_3\omega_{m}^{\frac{i}{2}}+a_5\omega_{m}^{\frac{i}{2}\cdot 2}+\cdots\end{aligned}\)
更簡潔的描述為
\(i<m,G(x_i')=a_0\omega_{m}^0+a_2\omega_{m}^{i}+a_4\omega_{m}^{2i}+\cdots\)
\(i<m,H(x_i')=a_1\omega_{m}^0+a_3\omega_{m}^{i}+a_5\omega_{m}^{2i}+\cdots\)
由於\(G(x'_i),H(x'_i)\)計算的是\([0,m-1]\)項，而求\(F(x_i)\)時用到的是\(0,2,4,\cdots\)項，實際需要訪問\(G(x^2_i),H(x^2_i)\)
和\(F(x_i)\)的式子比較，我們得到合並的式子為
\(F(x_i)=G(x^2_i)+x_i H(x^2_i)\)
帶入折半引理，實際等價於
\(F(x_i)=G(x'_i)+x_i H(x'_i)\)
注意\(x_i=x'_{i\mod m}\)
為了保證復雜度，盡量使得每次分治的子問題都分為兩部分，這樣的復雜度為\(O(n\log n)\)
附:實際上，分為\(d\)個子問題時，每次合並的復雜度為\(O(n\cdot d)\)，因此復雜度為
保證每次分治為兩個嚴格等大的子問題，可以從一開始就把\(n\)擴充為\(2\)的冪次
int N=1;
while(N<=n+m) N<<=1;
附:\(d\)個子問題時，設子問題答案為\(G_j(x_i)\)，則合並的式子為
\(\begin{aligned} F(x_i)=\sum_{j=0}^{d-1}x_i^jG_j(x_i^d)=\sum_{i=0}^{d-1}x_i^jG_j(x'_{i\mod \frac{n}{d}})\end{aligned}\)
點值式轉多項式
一個簡單的性質：單位根反演 \(\sum_{j=0}^{n-1}\omega_n^{ij}= \left\{\begin{aligned} \frac{\omega_n^{in}-1}{\omega_n^i-1}=0 && i\ne 0\\ n && i=0\end{aligned} \right.\)
設點值式對應\(y_i\)的序列為\(b_i\)
則\(n\cdot a_i=\sum_{j=0}^{n-1}\omega_n^{-ij} b_j\)
證明
\(\begin{aligned} \sum_{j=0}^{n-1}\omega_n^{-ij}b_j=\sum_{j=0}^{n-1} \omega_n^{-ij}(\sum_{k=0}^{n-1}a_k\omega_n^{jk})\end{aligned}\)
\(\begin{aligned} \sum_{j=0}^{n-1}\omega_n^{-ij}b_j= \sum_{k=0}^{n-1}a_k\sum_{j=0}^{n-1}\omega_n^{j(k-i)} \end{aligned}\)
由上面的式子，發現只有\(k-i=0\)時右邊的求和式有值，所以上式成立
因此點值式轉多項式直接把系數改為\(\omega_n^{-i}\)即可
Tips:
1.由於單位根的循環特性，溢出會直接溢出到本來的式子里
因此，如果乘法過后的多項式產生了超過\(>n\)的項\(x^i\)，會溢出到\(x^{i\mod n}\)
2.點值式並不是不滿足除法，只是除法得到的多項式並不一定是一個\(n\)元以內的多項式，除了恰好整除的情況，得到的通常是一個無窮級數的式子，如\(\begin{aligned} \frac{1}{1-x}=\frac{1-x^{\infty}}{1-x}=\sum_{i=0}^{\infty}x^i\end{aligned}\)
真正要求除法，通常是求前\(n\)項的結果，即需要用到多項式乘法逆
代碼實現與優化
然后我們得到一份優美的代碼(FFT)
(Complex是C++庫自帶的復數,M_PI是C++自帶\(\pi\)常量)
void FFT(int n,Complex *a,int f) {
if(n==1) return;
Complex tmp[N];
int m=n/2;
rep(i,0,m-1) tmp[i]=a[i<<1],tmp[i+m]=a[i<<1|1]; // 按照奇偶分類
memcpy(a,tmp,sizeof(Complex) * n);
FFT(m,a,f),FFT(m,a+m,f); // 分兩半，算g(x),h(x)
Complex w(cos(2*M_PI/n),f*sin(2*M_PI/n)),e(1,0); // w=x^1,e=x^i
rep(i,0,m-1) {
tmp[i]=a[i]+e*a[i+m]; // f(x_i)=g(x_i)+e*h(x_i)
e=e*w;
}
rep(i,m,n-1) {
tmp[i]=a[i-m]+e*a[i];
e=e*w;
}
memcpy(a,tmp,sizeof(T)*n);
}
由於\((\omega_n)^{\frac{n}{2}}=-1\)，所以還可以簡化為
Complex w(cos(2*M_PI/n),f*sin(2*M_PI/n)),e(1,0);
rep(i,0,m-1) {
tmp[i]=a[i]+e*a[i+m];
tmp[i+m]=a[i]-e*a[i+m];
e=e*w;
}
由於用了double，最后輸出要取整
蝴蝶優化
我們加一點優化，取代遞歸的分治過程
可以看到，分治時我們按照\(i \mod 2\)分成兩組，然后繼續分
這個過程中，實際上我們就是將\(i\)的二進制位前后翻轉
所以我們可以暴力處理出\(i\)分治底層的位置
rep(i,0,n-1) {
int x=i,s=0;
for(int j=1;(j<<c)<=n;++j) {
s=(s<<1)|(x&1);
x>>=1;
} // s就是最終位置
}
當然也是有\(O(n)\)處理方法的
int N=1,c=-1;
while(N<=n+m) N<<=1,c++;
rep(i,1,N-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<c);
(建議自己模擬一下)
有了這個翻轉數組，我們可以直接從分治底層開始解決整個問題，每次合並操作完全相同
每次分治問題的大小，依次合並每一個子問題區間即可
為了在一個數組上完成操作，還需要注意合並順序
代碼解釋\(i\):分治子問題大小為\(2i\)，\(l\):合並區間的左端點為\(l\)，右端點為\(l+2i\)，\(j\)枚舉合並位置
void FFT(int n,Complex *a){
for(int i=1;i<n;i<<=1) {
Complex w(cos(2*M_PI/n),f*sin(2*M_PI/n));
for(int l=0;l<n;l+=i*2) {
Complex e(1,0);
for(int j=l;j<l+i;++j,e=e*w) {
Complex t=a[j+m]*e; // a'[j]=a[j]+e*a[j+m]
// a'[j+i]=a[j]-e*a[j+m]
a[j+m]=a[j]-t;
a[j]=a[j]+t;
}
}
}
}
事實上我們還有更快的寫法，就是將\(\omega_n^i\)預處理出來
(注意這個\(\text{FFT}\)的預處理很考驗double精度，不能每次都直接累乘上去，隔幾個就要重新調用依次三角函數)
當然如果自己寫復數會更快
關於點值式轉多項式的優化
由於每次求得點值是\(\omega_n^{-i}=\omega_n^{n-i}\)
所以可以直接用 多項式轉點值式的函數， 最后把\([1,n-1]\)這一段翻轉，每個數除掉\(n\)即可
對於加減運算取模的優化
三目運算
a+=b,a=a>=P?a-P:a;
a-=b,a=a<0?a+P:a;
邏輯運算優化(原理是邏輯預算會在第一個確定表達式值的位置停下)
a+=b,((a>=P)&&(a-=P));
a-=b,((a<0)&&(a+=P));
關於系數預處理優化(以NTT為例)
帶入上面已經提到的優化，無預處理系數的\(\text{NTT}\)大概是這樣的
ll qpow(ll x,ll k=P-2){
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
int a[N];
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1) {
int w=qpow(3,(P-1)/i/2);
for(int l=0;l<n;l+=i*2) {
int e=1;
for(int j=l;j<l+i;++j,e=1ll*e*w%P) {
int t=1ll*a[j+i]*e%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}
一種簡單的預處理是，每次對於每個分治大小，預處理依次系數
ll qpow(ll x,ll k=P-2){
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
int a[N],e[N];
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
e[0]=1;
for(int i=1;i<n;i<<=1) {
int w=qpow(3,(P-1)/i/2);
for(int j=1;j<i;++j) e[j]=1ll*e[j-1]*w%P;
//for(int j=i-2;j>=0;j-=2) e[j+1]=1ll*w*(e[j]=e[j>>1])%P;
//這個版本是沿用上一次預處理的結果，實際(只有)用這種預處理方法可以極大程度上加強FFT的精度
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}
另一種是在一開始就把所有的系數用一個數組存下來，具體過程可以描述為
對於每個分治長度\(n\)，我們只需要訪問\(\omega_n^{0},\omega_n^{1},\cdots,\omega_n^{\frac{n}{2}-1}\)
那么對於分治長度\(n\)，我們在\(w\)數組的第\(\frac{n}{2}\) ~ \(n-1\)項依次存儲這些值
優化:我們只需要對於最大的分治長度處理，剩下的部分發現可以直接用折半引理訪問得到
ll qpow(ll x,ll k=P-2){
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
const int N=1<<21;
int a[N],w[N];
void Init(){
w[N>>1]=1;
int t=qpow(3,(P-1)/N);
rep(i,(N>>1)+1,N-1) w[i]=1ll*w[i-1]*t%P;
drep(i,(N>>1)-1,1) w[i]=w[i<<1];
}
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1) {
int *e=w+i;
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}
三份代碼在duck.ac上的評測結果表明，不預處理系數將近慢一倍
單組數據來看，預處理系數會慢一點
多組來看，預處理系數會快
實際差距不大，都可以使用
但是在某些層面來說，下面這份板子才是最好的(適用NTT,FFT且精度較高)，不需要預處理
ll qpow(ll x,ll k=P-2){
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
int a[N],e[N];
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
e[0]=1;
for(int i=1;i<n;i<<=1) {
int w=qpow(3,(P-1)/i/2);
for(int j=i-2;j>=0;j-=2) e[j+1]=1ll*w*(e[j]=e[j>>1])%P;
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}
拓展
1.分治+NTT
常用於處理多個計數背包的快速合並 (實際無權值01背包也是可以的)
我們可以用NTT\(n\log n\)合並兩個大小為\(n\)的背包
分治時，每次合並兩個分治子問題，總共的時間就是\(\sum size\log n\)
每個背包的\(size\)會被計算\(\log n\)次，所以總共復雜度是\(n \log ^2 n\)
2.CDQ+NTT
對於形如\(dp_i=\sum_{j=0}^{i-1}dp_jg_{i-j}\)的\(dp\)轉移(就是dp轉移與差值有關)
由於求\(dp_i\)時，需要保證\(dp_0,dp_1,\cdots,dp_{i-1}\)才能卷積，這個限制，我們可以用CDQ分治解決
對於當前分治區間\([L,R]\)
依次考慮\([L,mid]\)內部轉移，\([L,mid]\)向\([mid+1,R]\)的轉移(用FFT/NTT解決)，\([mid+1,R]\)內部轉移
算法流程
void Solve(l,r){
if(l==r) return;
mid=(l+r)>>1;
Solve(l,mid);
(l,mid)->(mid+1,r);
Solve(mid+1,r);
}
3.MTT(任意模數NTT)
4.\(n\)元點值式
練習建議：
1.高精度乘法
6.\(CDQ\)分治+\(FFT\)：HDU-5730 題解
7.\(CDQ\)+NTT/降次前綴和優化\(dp\)：HDU-5332 題解
9.圖上\(dp\)：
聯通圖個數：BZOJ-3456 題解
森林數量和帶限制森林數量：HDU - 5279 題解
10.點分治+FFT：CodeChef-PRIMEDST 題解
更多應用和優化參見毛嘯2016論文
(如：兩次FFT做卷積，4次FFT做MTT。。。)