Algorithm: 多項式乘法 Polynomial Multiplication: 快速傅里葉變換 FFT / 快速數論變換 NTT
Intro:
本篇博客將會從樸素乘法講起，通過分治乘法，到達FFT和NTThtml
旨在可以讓讀者（也讓本身）充分理解其思想c++
模板題入口：洛谷 P3803 【模板】多項式乘法（FFT）git
樸素乘法
約定：兩個多項式爲\(A(x)=\sum_{i=0}^{n}a_ix^i,B(x)=\sum_{i=0}^{m}b_ix^i\)
Prerequisite knowledge:
初中數學知識（手動滑稽）算法
最簡單的多項式方法就是逐項相乘再合併同類項，寫成公式：數組
若\(C(x)=A(x)B(x)\)，那麼\(C(x)=\sum_{i=0}^{n+m}c_ix^i\)，其中\(c_i=\sum_{j=0}^ia_jb_{i-j}\)app
因而一個樸素乘法就產生了，見代碼（利用某種喪心病狂的方式省了\(b\)數組）性能
//This program is written by Brian Peng.
#pragma GCC optimize("Ofast","inline","no-stack-protector")
#include<bits/stdc++.h>
using namespace std;
#define Rd(a) (a=read())
#define Gc(a) (a=getchar())
#define Pc(a) putchar(a)
int read(){
register int x;register char c(getchar());register bool k;
while(!isdigit(c)&&c^'-')if(Gc(c)==EOF)exit(0);
if(c^'-')k=1,x=c&15;else k=x=0;
while(isdigit(Gc(c)))x=(x<<1)+(x<<3)+(c&15);
return k?x:-x;
}
void wr(register int a){
if(a<0)Pc('-'),a=-a;
if(a<=9)Pc(a|'0');
else wr(a/10),Pc((a%10)|'0');
}
signed const INF(0x3f3f3f3f),NINF(0xc3c3c3c3);
long long const LINF(0x3f3f3f3f3f3f3f3fLL),LNINF(0xc3c3c3c3c3c3c3c3LL);
#define Ps Pc(' ')
#define Pe Pc('\n')
#define Frn0(i,a,b) for(register int i(a);i<(b);++i)
#define Frn1(i,a,b) for(register int i(a);i<=(b);++i)
#define Frn_(i,a,b) for(register int i(a);i>=(b);--i)
#define Mst(a,b) memset(a,b,sizeof(a))
#define File(a) freopen(a".in","r",stdin),freopen(a".out","w",stdout)
#define N (2000010)
int n,m,a[N],b,c[N];
signed main(){
Rd(n),Rd(m);
Frn1(i,0,n)Rd(a[i]);
Frn1(i,0,m){Rd(b);Frn1(j,0,n)c[i+j]+=b*a[j];}
Frn1(i,0,n+m)wr(c[i]),Ps;
exit(0);
}
Time complexity: \(O(nm)\)，若是\(m=O(n)\)，則爲\(O(n^2)\)優化
Memory complexity: \(O(n)\)ui
看看效果spa
意料之中，因此必須優化
樸素分治乘法
P.s 這一部分講述了FFT的分治方法，與FFT仍是有區別的，若是已經理解的能夠跳過
約定：\(n\)爲同時屬於\(A(x),B(x)\)次數界的最小的\(2\)的正整數冪，並將兩個多項式設爲\(A(x)=\sum_{i=0}^{n-1}a_ix^i,B(x)=\sum_{i=0}^{n-1}b_ix^i\)，不存在的係數補零
次數界：嚴格\(>\)一個多項式次數的整數（E.g 多項式\(P(x)=x^2+x+1\)的次數界爲\(\geqslant3\)的全部整數）
Reference：
《算法導論》
Prerequisite knowledge:
分治思想
如今來考慮如何去優化乘法
嘗試將兩個多項式按照未知項次數的奇偶性分開：
\(A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2),B(x)=B^{[0]}(x^2)+xB^{[1]}(x^2)\)
其中\(A^{[0]}(x)=\sum_{i=0}^{n/2-1}a_{2i}x^i,A^{[1]}(x)=\sum_{i=0}^{n/2-1}a_{2i+1}x^i\)，\(B^{[0]}(x)\)與\(B^{[1]}(x)\)同理
因而兩個多項式就被拆成了兩個次數界爲\(n/2\)的四個多項式啦：
P.s 如下的公式中，用\(A\)表示\(A(x)\)，\(A^{[0]}\)和\(A^{[1]}\)分別表示\(A^{[0]}(x^2)\)和\(A^{[1]}(x^2)\)，\(B\)同理
\(AB=(A^{[0]}+xA^{[1]})(B^{[0]}+xB^{[1]})=A^{[0]}B^{[0]}+x(A^{[1]}B^{[0]}+A^{[0]}B^{[1]})+x^2A^{[1]}B^{[1]}\)
在此能夠發現一種分治算法：把兩個多項式折半，而後再遞歸算\(4\)次多項式乘法，最後合併加起來（反正多項式加法是\(O(n)\)的）
P.s 注意合併方式：\(A^{[0]}\)和\(A^{[1]}\)分別表示\(A^{[0]}(x^2)\)和\(A^{[1]}(x^2)\)，因此是交錯的，見代碼
（爲了省空間用了vector）
//This program is written by Brian Peng.
#pragma GCC optimize("Ofast","inline","no-stack-protector")
#include<bits/stdc++.h>
using namespace std;
#define Rd(a) (a=read())
#define Gc(a) (a=getchar())
#define Pc(a) putchar(a)
int read(){
register int x;register char c(getchar());register bool k;
while(!isdigit(c)&&c^'-')if(Gc(c)==EOF)exit(0);
if(c^'-')k=1,x=c&15;else k=x=0;
while(isdigit(Gc(c)))x=(x<<1)+(x<<3)+(c&15);
return k?x:-x;
}
void wr(register int a){
if(a<0)Pc('-'),a=-a;
if(a<=9)Pc(a|'0');
else wr(a/10),Pc((a%10)|'0');
}
signed const INF(0x3f3f3f3f),NINF(0xc3c3c3c3);
long long const LINF(0x3f3f3f3f3f3f3f3fLL),LNINF(0xc3c3c3c3c3c3c3c3LL);
#define Ps Pc(' ')
#define Pe Pc('\n')
#define Frn0(i,a,b) for(register int i(a);i<(b);++i)
#define Frn1(i,a,b) for(register int i(a);i<=(b);++i)
#define Frn_(i,a,b) for(register int i(a);i>=(b);--i)
#define Mst(a,b) memset(a,b,sizeof(a))
#define File(a) freopen(a".in","r",stdin),freopen(a".out","w",stdout)
typedef vector<int> Vct;
int n,m,s;
Vct a,b,c;
void add(Vct&a,Vct&b,Vct&c){Frn0(i,0,c.size())c[i]=a[i]+b[i];}
void mlt(Vct&a,Vct&b,Vct&c,int n);
signed main(){
Rd(n),Rd(m),a.resize(s=1<<int(log2(max(n,m))+1)),b.resize(s),c.resize(s<<1);
Frn1(i,0,n)Rd(a[i]);
Frn1(i,0,m)Rd(b[i]);
mlt(a,b,c,s);
Frn1(i,0,n+m)wr(c[i]),Ps;
exit(0);
}
void mlt(Vct&a,Vct&b,Vct&c,int n){
int n2(n>>1);
Vct a0(n2),a1(n2),b0(n2),b1(n2),ab0(n),ab1(n),abm(n);
if(n==1){c[0]=a[0]*b[0];return;}
Frn0(i,0,n2)a0[i]=a[i<<1],a1[i]=a[i<<1|1],b0[i]=b[i<<1],b1[i]=b[i<<1|1];
mlt(a0,b0,ab0,n2),mlt(a1,b1,ab1,n2);
Frn0(i,0,n)c[i<<1]=ab0[i]+(i?ab1[i-1]:0);
mlt(a0,b1,ab0,n2),mlt(a1,b0,ab1,n2),add(ab0,ab1,abm);
Frn0(i,0,n-1)c[i<<1|1]=abm[i];
}
看看效果
好像更慘……
爲何呢，由於這個算法的時間複雜度仍是\(O(n^2)\)的，具體證實以下
\(T(n)=4T(n/2)+f(n)\)，其中\(f(n)=O(n)\)（就是\(n\)位加法的時間）
運用主方法，\(a=4,b=2,log_ba=log_2 4=2>1\)，因此\(T(n)=O(n^{log_ba})=O(n^2)\)
並且不只複雜度高，常數因子也由於遞歸變高了
因此繼續優化吧……
分治乘法
接上上一部分的內容，考慮如何優化時間複雜度
先來一個小插曲：如何只作\(3\)次乘法，求出線性多項式\(ax+b\)與\(cx+d\)的乘積
先看看結果：\((ax+b)(cx+d)=acx^2+(ad+bc)x+bd\)，總共有\(4\)次乘法
因此若是隻用\(3\)次乘法，那麼\(ad+bc\)必須只能用一次乘法獲得
嘗試把\(3\)個係數加起來，就是\(ac+ad+bc+bd=(a+b)(c+d)\)
答案出來了，用\(3\)次乘法分別算出\(ac,bd\)與\((a+b)(c+d)\)，那麼中間項係數\(=(a+b)(c+d)-ac-bd\)
回到原題目
\(AB=(A^{[0]}+xA^{[1]})(B^{[0]}+xB^{[1]})=A^{[0]}B^{[0]}+x(A^{[1]}B^{[0]}+A^{[0]}B^{[1]})+x^2A^{[1]}B^{[1]}\)
因而中間項也可使用相似的方法：\(A^{[1]}B^{[0]}+A^{[0]}B^{[1]}=(A^{[0]}+A^{[1]})(B^{[0]}+B^{[1]})-A^{[0]}B^{[0]}-A^{[1]}B^{[1]}\)
成功減小一次乘法運算，見代碼
//This program is written by Brian Peng.
#pragma GCC optimize("Ofast","inline","no-stack-protector")
#include<bits/stdc++.h>
using namespace std;
#define Rd(a) (a=read())
#define Gc(a) (a=getchar())
#define Pc(a) putchar(a)
int read(){
register int x;register char c(getchar());register bool k;
while(!isdigit(c)&&c^'-')if(Gc(c)==EOF)exit(0);
if(c^'-')k=1,x=c&15;else k=x=0;
while(isdigit(Gc(c)))x=(x<<1)+(x<<3)+(c&15);
return k?x:-x;
}
void wr(register int a){
if(a<0)Pc('-'),a=-a;
if(a<=9)Pc(a|'0');
else wr(a/10),Pc((a%10)|'0');
}
signed const INF(0x3f3f3f3f),NINF(0xc3c3c3c3);
long long const LINF(0x3f3f3f3f3f3f3f3fLL),LNINF(0xc3c3c3c3c3c3c3c3LL);
#define Ps Pc(' ')
#define Pe Pc('\n')
#define Frn0(i,a,b) for(register int i(a);i<(b);++i)
#define Frn1(i,a,b) for(register int i(a);i<=(b);++i)
#define Frn_(i,a,b) for(register int i(a);i>=(b);--i)
#define Mst(a,b) memset(a,b,sizeof(a))
#define File(a) freopen(a".in","r",stdin),freopen(a".out","w",stdout)
typedef vector<int> Vct;
int n,m,s;
Vct a,b,c;
void add(Vct&a,Vct&b,Vct&c){Frn0(i,0,c.size())c[i]=a[i]+b[i];}
void mns(Vct&a,Vct&b,Vct&c){Frn0(i,0,c.size())c[i]=a[i]-b[i];}
void mlt(Vct&a,Vct&b,Vct&c);
signed main(){
Rd(n),Rd(m),a.resize(s=1<<int(log2(max(n,m))+1)),b.resize(s),c.resize(s<<1);
Frn1(i,0,n)Rd(a[i]);
Frn1(i,0,m)Rd(b[i]);
mlt(a,b,c);
Frn1(i,0,n+m)wr(c[i]),Ps;
exit(0);
}
void mlt(Vct&a,Vct&b,Vct&c){
int n(a.size()),n2(a.size()>>1);
Vct a0(n2),a1(n2),b0(n2),b1(n2),ab0(n),ab1(n),abm(n);
if(n==1){c[0]=a[0]*b[0];return;}
Frn0(i,0,n2)a0[i]=a[i<<1],a1[i]=a[i<<1|1],b0[i]=b[i<<1],b1[i]=b[i<<1|1];
mlt(a0,b0,ab0),mlt(a1,b1,ab1);
Frn0(i,0,n)c[i<<1]=ab0[i]+(i?ab1[i-1]:0);
add(a0,a1,a0),add(b0,b1,b0),mlt(a0,b0,abm),mns(abm,ab0,abm),mns(abm,ab1,abm);
Frn0(i,0,n-1)c[i<<1|1]=abm[i];
}
看看效果
比樸素分治乘法好一點，可是仍是沒樸素乘法強，仍是很慘
看看這個算法的時間複雜度：
\(T(n)=3T(n/2)+f(n)\)，其中\(f(n)=O(n)\)
運用主方法，\(a=3,b=2,\log_ba=\log_2 3\approx1.58>1\)，因此\(T(n)=O(n^{\log_ba})=O(n^{\log_2 3})\)
額那不是應該比樸素算法要好嗎，這是什麼狀況
Reason 1. 分治乘法的常數因子太大
Reason 2. 打開\(\#5\)數據一看，\(n=1,m=3e6\)，那麼\(O(n^{\log_2 3})\)的分治乘法也頂不過\(O(nm)\)的樸素乘法啊……
因此就要請上本文的主角了
快速傅里葉變換 FFT (Fast Fourier Transform)
Fairly Frightening Transform
約定：\(n\)爲屬於\(A(x),B(x)\)的乘積\(C(x)\)次數界的最小的\(2\)的正整數冪（即知足\(>\)輸入\(n+m\)的最小的\(2\)的正整數冪），並一樣將兩個多項式設爲\(A(x)=\sum_{i=0}^{n-1}a_ix^i,B(x)=\sum_{i=0}^{n-1}b_ix^i\)
Reference：
《算法導論》
Prerequisite knowledge:
分治思想
複數的基本知識
線性代數的基本知識
Part 1: 多項式的兩種表示方式
1. 係數表達
對一個次數界爲\(n\)的多項式\(A(x)=\sum_{i=0}^{n-1}a_ix^i\)，其係數表達是向量\(\pmb{a}=\left[\begin{matrix}a_0\\a_1\\\vdots\\a_{n-1}\end{matrix} \right]\)
使用係數表達時，下列操做的時間複雜度：
求值\(O(n)\)
加法\(O(n)\)
乘法樸素\(O(n^2)\)，優化\((n^{\log_2 3})\)（即分治乘法）
P.s 當多項式\(C(x)=A(x)B(x)\)時，\(\pmb{c}\)被稱爲\(\pmb{a}\)與\(\pmb{b}\)的卷積(convolution)，記爲\(\pmb{c}=\pmb{a}\bigotimes\pmb{b}\)
2. 點值表達
一個次數界爲\(n\)的多項式\(A(x)\)的點值表達是一個有\(n\)個點值對所組成的集合\(\{(x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1})\}\)
進行\(n\)次求值就能夠把係數表達轉化爲點值表達，總時間\(O(n^2)\)，用公式表示就是：
\(\left[\begin{matrix}1&x_0&x_0^2&\cdots&x_0^{n-1}\\1&x_1&x_1^2&\cdots&x_1^{n-1}\\\vdots&\vdots&\vdots&\ddots&\vdots\\1&x_{n-1}&x_{n-1}^2&\cdots&x_{n-1}^{n-1}\end{matrix} \right]\left[\begin{matrix}a_0\\a_1\\\vdots\\a_{n-1}\end{matrix} \right]=\left[\begin{matrix}y_0\\y_1\\\vdots\\y_{n-1}\end{matrix} \right]\)
其中左邊的矩陣表示爲\(V(x_0,x_1,\cdots,x_{n-1})\)稱爲範德蒙德矩陣，因而能夠將公式簡化爲\(V(x_0,x_1,\cdots,x_{n-1})\pmb{a}=\pmb{y}\)
使用拉格朗日公式，能夠在\(O(n^2)\)時間將點值表達轉化爲係數表達，該過程稱爲插值
對於兩個在相同位置求值的點值表達多項式，下列操做的時間複雜度：
加法\(O(n)\)（只要將各個位置的\(y\)值相加便可）
乘法\(O(n)\)（同理）
因此這就是使用FFT的緣由：經過精心選取\(x\)值，能夠在\(O(n\log n)\)時間完成求值，再\(O(n)\)乘法，最後\(O(n\log n)\)插值
傅里葉大神究竟選了什麼神奇的\(x\)值呢，請看
Part 2: 單位複數根及其性質
\(n\)次單位複數根是知足\(\omega^n=1\)的複數\(\omega\)，正好有\(n\)個，記爲：
\(\omega_n^k=e^{2\pi ik/n}=\cos(2\pi k/n)+i\sin(2\pi k/n)\)
其中\(\omega_n=e^{2\pi i/n}=\cos(2\pi /n)+i\sin(2\pi /n)\)被稱爲主\(n\)次單位根，全部其餘\(n\)次單位根都是\(\omega_n\)的冪次
能夠把\(n\)個單位根看做是複平面上以單位圓的\(n\)個等分點爲終點的向量，具體緣由就是複數乘法「模長相乘，輻角相加」的規律
如圖表示的是\(8\)次單位複數根在複平面上的位置
因而就能夠獲得規律：\(\omega_n^j\omega_n^k=\omega_n^{j+k}=\omega_n^{(j+k)\mod n}\)，相似地\(\omega_n^{-1}=\omega_n^{n-1}\)
接下來的三個引理就是FFT的重頭戲啦
1. 消去引理：對任何整數\(n\geqslant 0,k\geqslant 0,d>0\)，有\(\omega_{dn}^{dk}=\omega_n^k\)
Proof: \(\omega_{dn}^{dk}=(e^{2\pi i/dn})^{dk}=(e^{2\pi i/n})^k=\omega_n^k\)
2. 折半引理：對任何偶數\(n\)和整數\(k\)，有\((\omega_n^k)^2=(\omega_n^{k+n/2})^2=\omega_{n/2}^k\)
Proof: \((\omega_n^k)^2=\omega_n^{2k},(\omega_n^{k+n/2})^2=\omega_n^{2k+n}=\omega_n^{2k}\)，最後用消去引理，\(\omega_n^{2k}=\omega_{n/2}^k\)
3. 求和引理：對任何整數\(n\geqslant 0\)與非負整數\(k:n\nmid k\)，有\(\sum_{j=0}^{n-1}(\omega_n^k)^j=0\)
Proof: 利用等比數列求和公式，\(\sum_{j=0}^{n-1}(\omega_n^k)^j=\frac{1-(\omega_n^k)^n}{1-\omega_n^k}=\frac{1-\omega_n^{nk}}{1-\omega_n^k}=\frac{1-1}{1-\omega_n^k}=0\)，爲了使分母\(1-\omega_n^k\neq 0\)，必須知足\(\omega_n^k\neq 1\implies n\nmid k\)
Part 3: 離散傅里葉變換 DFT (Discrete Fourier Transform)
DFT就是將次數界爲\(n\)的多項式\(A(x)\)在\(n\)次單位複數根上求值的過程
簡化一下表示方法：
\(V_n=V(\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1})=\left[\begin{matrix}1&1&1&1&\cdots&1\\1&\omega_n&\omega_n^2&\omega_n^3&\cdots&\omega_n^{n-1}\\1&\omega_n^2&\omega_n^4&\omega_n^6&\cdots&\omega_n^{2(n-1)}\\1&\omega_n^3&\omega_n^6&\omega_n^9&\cdots&\omega_n^{3(n-1)}\\\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\1&\omega_n^{n-1}&\omega_n^{2(n-1)}&\omega_n^{3(n-1)}&\cdots&\omega_n^{(n-1)(n-1)}\end{matrix} \right]\)
用公式表示就是\(V_n\pmb{a}=\pmb{y}\)，也記爲\(\pmb{y}=DFT_n(\pmb a)\)
另外，能夠發現\([V_n]_{ij}=\omega_n^{ij}\implies y_i=\sum_{j=0}^{n-1}[V_n]_{ij}a_j=\sum_{j=0}^{n-1}\omega_n^{ij}a_j\)
終於能夠看看具體操做了
Part 4: FFT
FFT利用單位根的特殊性質把DFT優化到了\(O(n\log n)\)
和分治乘法同樣，按未知項次數的奇偶性分開：\(A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)\)
其中\(A^{[0]}(x)=\sum_{i=0}^{n/2-1}a_{2i}x^i,A^{[1]}(x)=\sum_{i=0}^{n/2-1}a_{2i+1}x^i\)
這時，求\(A(x)\)在\(\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}\)的值變成了：
1. 求\(A^{[0]}(x)\)與\(A^{[1]}(x)\)在\((\omega_n^0)^2,(\omega_n^1)^2,\cdots,(\omega_n^{n-1})^2\)的值
根據折半引理，\((\omega_n^0)^2,(\omega_n^1)^2,\cdots,(\omega_n^{n-1})^2\)中兩兩重複，其實就是\(n/2\)個\(n/2\)次單位根
因此只要對拆開的兩個多項式分別作\(DFT_{n/2}\)便可，獲得\(\pmb y^{[0]}\)與\(\pmb y^{[1]}\)
2. 合併答案
\(\omega_n^{n/2}=e^{2\pi i (n/2)/n}=e^{\pi i}=-1\)（根據傳說中的最美公式\(e^{i\pi}+1=0\)）
因此\(\omega_n^{k+n/2}=\omega_n^k\omega_n^{n/2}=-\omega_n^k\)
因此\(y_i=y^{[0]}_i+\omega_n^i y^{[1]}_i,y_{i+n/2}=y^{[0]}_i-\omega_n^i y^{[1]}_i,i=0,1,\cdots,n/2-1\)
具體運行時，就每次循環結束時讓一個初始爲\(1\)的變量\(*\omega_n\)便可
遞歸邊界：\(n=1\)，那麼\(w_1^0 a_0=a_0\)，因此直接返回自身
計算一下時間複雜度
\(T(n)=2T(n/2)+f(n)\)，其中\(f(n)=O(n)\)（合併答案）
運用主方法，\(a=2,b=2,\log_ba=\log_2 2=1\)，因此\(T(n)=O(n^{\log_ba}\log n)=O(n\log n)\)（皆大歡喜）
Part 5: 離散傅里葉逆變換
可別高興太早，還有插值哦
由於\(\pmb{y}=DFT_n(\pmb{a})=V_n\pmb{a}\)，因此\(\pmb{a}=V_n^{-1}\pmb{y}\)，記爲\(\pmb{a}=DFT_n^{-1}(\pmb{y})\)
定理：對\(i,j=0,1,\cdots,n-1\)，有\([V_n^{-1}]_{ij}=\omega_n^{-ij}/n\)
Proof: 證實\(V_n^{-1}V_n=I_n\)便可
\([V_n^{-1}V_n]_{ij}=\sum_{k=0}^{n-1}(\omega_n^{-ik}/n)\omega_n^{kj}=\frac{\sum_{k=0}^{n-1}\omega_n^{-ik}\omega_n^{kj}}{n}=\frac{\sum_{k=0}^{n-1}\omega_n^{(j-i)k}}{n}\)
若是\(i=j\)，則該值\(=\frac{\sum_{k=0}^{n-1}\omega_n^0}{n}=n/n=1\)
不然，由於\(n\nmid k\)，根據求和引理，該值\(=0/n=0\)，因此構成了\(I_n\)
接下來\(\pmb{a}=DFT_n^{-1}(\pmb{y})=V_n^{-1}\pmb{y}\implies a_i=\sum_{j=0}^{n-1}[V_n^{-1}]_{ij}y_j=\sum_{j=0}^{n-1}(\omega_n^{-ij}/n)y_j=\frac{\sum_{j=0}^{n-1}\omega_n^{-ij}y_j}{n}\)
比較一下DFT中\(y_i=\sum_{j=0}^{n-1}\omega_n^{ij}a_j\)
只要運算時把\(\omega_n\)換成\(-\omega_n\)，而後將最終答案\(/n\)，就把DFT變成逆DFT了
終於能夠來到激動人心的實現環節了
Part 6: 遞歸實現
根據前文，只要將分治乘法的代碼修改一下便可
能夠作到直接在原址進行FFT，就是將分開的兩個多項式分置在左右兩邊
STL提供了現成的complex類可供使用
代碼中用iv表示是否爲逆DFT，用o存儲主單位根，用w累積
P.s 最後別忘了\(/n\)，並且\(+0.5\)爲了四捨五入提升精度
//This program is written by Brian Peng.
#pragma GCC optimize("Ofast","inline","no-stack-protector")
#include<bits/stdc++.h>
using namespace std;
#define Rd(a) (a=read())
#define Gc(a) (a=getchar())
#define Pc(a) putchar(a)
int read(){
register int u;register char c(getchar());register bool k;
while(!isdigit(c)&&c^'-')if(Gc(c)==EOF)exit(0);
if(c^'-')k=1,u=c&15;else k=u=0;
while(isdigit(Gc(c)))u=(u<<1)+(u<<3)+(c&15);
return k?u:-u;
}
void wr(register int a){
if(a<0)Pc('-'),a=-a;
if(a<=9)Pc(a|'0');
else wr(a/10),Pc((a%10)|'0');
}
signed const INF(0x3f3f3f3f),NINF(0xc3c3c3c3);
long long const LINF(0x3f3f3f3f3f3f3f3fLL),LNINF(0xc3c3c3c3c3c3c3c3LL);
#define Ps Pc(' ')
#define Pe Pc('\n')
#define Frn0(i,a,b) for(register int i(a);i<(b);++i)
#define Frn1(i,a,b) for(register int i(a);i<=(b);++i)
#define Frn_(i,a,b) for(register int i(a);i>=(b);--i)
#define Mst(a,b) memset(a,b,sizeof(a))
#define File(a) freopen(a".in","r",stdin),freopen(a".out","w",stdout)
double const Pi(acos(-1));
typedef complex<double> Cpx;
#define N (2100000)
Cpx o,w,a[N],b[N],tmp[N],x,y;
int n,m,s;
bool iv;
void fft(Cpx*a,int n);
signed main(){
Rd(n),Rd(m),s=1<<int(log2(n+m)+1);
Frn1(i,0,n)Rd(a[i]);
Frn1(i,0,m)Rd(b[i]);
fft(a,s),fft(b,s);
Frn0(i,0,s)a[i]*=b[i];
iv=1,fft(a,s);
Frn1(i,0,n+m)wr(a[i].real()/s+0.5),Ps;
exit(0);
}
void fft(Cpx*a,int n){
if(n==1)return;
int n2(n>>1);
Frn0(i,0,n2)tmp[i]=a[i<<1],tmp[i+n2]=a[i<<1|1];
copy(tmp,tmp+n,a),fft(a,n2),fft(a+n2,n2);
o={cos(Pi/n2),(iv?-1:1)*sin(Pi/n2)},w=1;
Frn0(i,0,n2)x=a[i],y=w*a[i+n2],a[i]=x+y,a[i+n2]=x-y,w*=o;
}
Time complexity: \(O(n\log n)\)
Memory complexity: \(O(n)\)
看看效果
性能已經超過了樸素乘法（必然的），可是仍是沒有AC
注意到\(n,m\leqslant 1e6\)，因此不只要讓時間複雜度至少\(O(n\log n)\)，還要保持小的常數因子，總之遞歸還不夠快
Part 6: 迭代實現
設\(l=\lceil\log_2(n+m+1)\rceil,s=2^l\)，那麼\(A(x),B(x),A(x)B(x)\)都是次數界爲\(s\)的多項式
如今須要尋找到一種迭代的方式，使答案自底向上合併以減小常數因子
仍是像遞歸版同樣，把\(A^{[0]}(x)\)放在左邊，\(A^{[1]}(x)\)放在右邊
觀察每一層遞歸時各個係數所在位置的規律，以\(s=8\)爲例
0-> 0 1 2 3 4 5 6 7 1-> 0 2 4 6|1 3 5 7 2-> 0 4|2 6|1 5|3 7 end 0|4|2|6|1|5|3|7
沒看出來？那就拆成二進制看看
0-> 000 001 010 011 100 101 110 111 1-> 000 010 100 110|001 011 101 111 2-> 000 100|010 110|001 101|011 111 end 000|100|010|110|001|101|011|111
顯然地在最後一層遞歸時，係數編號正好是位置編號的反轉（更準確的說是前\(l\)位的反轉）
一個較爲感性的Proof: 由於是按照奇偶性分類，也就是說在第\(i\)層遞歸時判斷的是該編號二進制第\(i\)位（從零開始），爲\(0\)放左邊，\(1\)放右邊，而放右邊的結果就是它的位置編號的二進制第\(l-i-1\)位是\(1\)
因此到了遞歸最底層，位置編號的二進制就正好是係數編號二進制前\(l\)位的反轉啦
構造數組\(r_{0..s-1}\)，其中\(r_i\)表示\(i\)二進制前\(l\)位的反轉，能夠利用遞推完成，具體請本身思考（或看代碼）
蝴蝶操做 (Butterfly Operation)
其實在遞歸版代碼中已經出現，可是這裏再詳細說明一下
還記得\(y_i=y^{[0]}_i+\omega_n^i y^{[1]}_i,y_{i+n/2}=y^{[0]}_i-\omega_n^i y^{[1]}_i,i=0,1,\cdots,n/2-1\)嗎？
可是如今不使用\(\pmb{y}\)，而是\(\pmb{a}\)直接合並
由於按照奇偶性分置在兩邊，因此\(a^{[0]}_i=a_i,a^{[1]}_i=a_{i+n/2}\)
設\(x=a_i,y=\omega_n^i a_{i+n/2}\)
那麼新的\(a_i=x+y,a_{i+n/2}=x-y\)
這就是蝴蝶操做啦
有了蝴蝶操做，只要將全部係數按照\(r\)數組的位置排列，再迭代合併，就完成了FFT
在代碼中，用\(i,i_2\)表示當前合併產生的和開始的序列長度，\(j\)表示合併序列的開頭位置，\(k\)控制每一位的合併，上代碼
//This program is written by Brian Peng.
#pragma GCC optimize("Ofast","inline","no-stack-protector")
#include<bits/stdc++.h>
using namespace std;
#define Rd(a) (a=read())
#define Gc(a) (a=getchar())
#define Pc(a) putchar(a)
int read(){
register int u;register char c(getchar());register bool k;
while(!isdigit(c)&&c^'-')if(Gc(c)==EOF)exit(0);
if(c^'-')k=1,u=c&15;else k=u=0;
while(isdigit(Gc(c)))u=(u<<1)+(u<<3)+(c&15);
return k?u:-u;
}
void wr(register int a){
if(a<0)Pc('-'),a=-a;
if(a<=9)Pc(a|'0');
else wr(a/10),Pc((a%10)|'0');
}
signed const INF(0x3f3f3f3f),NINF(0xc3c3c3c3);
long long const LINF(0x3f3f3f3f3f3f3f3fLL),LNINF(0xc3c3c3c3c3c3c3c3LL);
#define Ps Pc(' ')
#define Pe Pc('\n')
#define Frn0(i,a,b) for(register int i(a);i<(b);++i)
#define Frn1(i,a,b) for(register int i(a);i<=(b);++i)
#define Frn_(i,a,b) for(register int i(a);i>=(b);--i)
#define Mst(a,b) memset(a,b,sizeof(a))
#define File(a) freopen(a".in","r",stdin),freopen(a".out","w",stdout)
double const Pi(acos(-1));
typedef complex<double> Cpx;
#define N (2100000)
Cpx a[N],b[N],o,w,x,y;
int n,m,l,s,r[N];
void fft(Cpx*a,bool iv);
signed main(){
Rd(n),Rd(m),s=1<<(l=log2(n+m)+1);
Frn1(i,0,n)Rd(a[i]);
Frn1(i,0,m)Rd(b[i]);
Frn0(i,0,s)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
fft(a,0),fft(b,0);
Frn0(i,0,s)a[i]*=b[i];
fft(a,1);
Frn1(i,0,n+m)wr(a[i].real()+0.5),Ps;
exit(0);
}
void fft(Cpx*a,bool iv){
Frn0(i,0,s)if(i<r[i])swap(a[i],a[r[i]]);
for(int i(2),i2(1);i<=s;i2=i,i<<=1){
o={cos(Pi/i2),(iv?-1:1)*sin(Pi/i2)};
for(int j(0);j<s;j+=i){
w=1;
Frn0(k,0,i2){
x=a[j+k],y=w*a[j+k+i2];
a[j+k]=x+y,a[j+k+i2]=x-y,w*=o;
}
}
}
if(iv)Frn0(i,0,s)a[i]/=s;
}
Time complexity: \(O(n\log n)\)
Memory complexity: \(O(n)\)
看看效果
終於……
到如今爲止FFT的內容已經所有結束啦，下面是拓展部分
Extension: 快速數論變換 NTT (Number Theoretic Transform)
雖然FFT具備優秀的時間複雜度，但由於用到了複數，不可避免會出現精度問題
若是多項式係數和結果都是必定範圍非負整數，能夠考慮使用NTT來優化精度和時空常數
Reference：
《算法導論》
Prerequisite knowledge:
FFT（必須知道的）
模運算基本知識
原根的性質
如今考慮全部運算都在\(mod P\)意義下
設有正整數\(g\)，若是在\(g\)的次冪可以獲得\(<P\)的任何正整數，那麼稱\(g\)是\(Z_P^*\)的原根，其中\(Z_P^*\)是模\(P\)乘法羣，在這裏很少做解釋
E.g 對於\(P=7\)，計算全部\(<P\)的正整數的次冪構成的集合
1-> {1}
2-> {1,2,4}
3-> {1,2,3,4,5,6}
4-> {1,2,4}
5-> {1,2,3,4,5,6}
6-> {1,6}
因此\(3,5\)就是\(Z_7^*\)的原根
在代碼中，通常使用大質數\(P=998244353,g=3\)
原根的特色就是它的次冪以長度爲\(\phi(P)\)循環
E.g \(P=7,g=3\)
那麼\(g\)的次冪（從\(g^0\)）開始分別是：\(1,3,2,6,4,5,1,3,2,6,4,5,\cdots\)
這個特性和單位根很是類似
可是要徹底替換單位根，還差一步
單位根的代替品
在FFT中使用的是循環長度爲\(n\)，且知足消去引理和求和引理的\(n\)次單位複數根（折半引理由消去引理推出，故不考慮）
因此爲了讓循環長度爲\(n\)，咱們不直接使用原根，而是原根的次冪
離散對數定理：若是\(g\)是\(Z_P^*\)的一個原根，則\(x\equiv y(\mod\phi(P))\iff g^x\equiv g^y(\mod P)\)
Proof: 設\(x\equiv y(\mod\phi(P))\)，則對某個整數\(k\)有\(x=y+k\phi(P)\)
所以\(g^x\equiv g^{y+k\phi(P)} \equiv g^y (g^{\phi(P)})^k \equiv g^y 1^k \equiv g^y (\mod P)\)（根據歐拉定理）
反過來，由於循環長度是\(\phi(P)\)，一定有\(x\equiv y(\mod\phi(P))\)
如今考慮有一個\(g\)的次冪\(g^q\)知足\(g^q\)的次冪以長度\(n\)循環
也就是說對任意整數\(x\geqslant0,0<y<n\)，有\(g^{qx}\equiv g^{q(x+n)}\not\equiv g^{q(x+y)}(\mod P)\)
即\(qx\equiv q(x+n)\not\equiv q(x+y)(\mod \phi(P))\)
即\(0\equiv qn\not\equiv qy(\mod \phi(P))\)
可得\(\phi(P)|qn,\phi(P)\nmid qy\)
那麼爲了使\(q\)的因數數量最小化，\(q=\phi(P)/n\)
此時\(qy-\phi(P)y/n\)
由於\(y<n\)，因此\(qy<\phi(P)\)，一定有\(\phi(P)\nmid qy\)
因此\(q=\phi(P)/n\)是可取的
那麼問題來了，萬一\(\phi(P)/n\)不是整數怎麼辦？
這就引出了大質數\(P=998244353\)的另外一個性質：\(\phi(P)=P-1=998244352=2^{23}\cdot 7\cdot 17\)
而根據數據範圍\(n,m\leqslant1e6\)，可知\(l\leqslant 21\)，因此\(q\)老是整數（真是一個神奇的數字）
總結一下，\(g^{\phi(P)/n}=g^{\frac{P-1}{n}}\)就是\(\omega_n\)的代替者
消去引理（只需考慮\(d=2\)的狀況）和求和引理就請你們本身證實了（其實道理都很是類似）
最後只要把\(\omega_n\)替換掉，運算都改成模\(P\)意義下的運算便可，在算\(-\omega_n\)時要用到\(g^{-1}=332748118\)，還有最後答案別忘了\(*s^{-1}\)，上代碼
//This program is written by Brian Peng.
#pragma GCC optimize("Ofast","inline","no-stack-protector")
#include<bits/stdc++.h>
using namespace std;
#define int long long
#define Rd(a) (a=read())
#define Gc(a) (a=getchar())
#define Pc(a) putchar(a)
int read(){
register int u;register char c(getchar());register bool k;
while(!isdigit(c)&&c^'-')if(Gc(c)==EOF)exit(0);
if(c^'-')k=1,u=c&15;else k=u=0;
while(isdigit(Gc(c)))u=(u<<1)+(u<<3)+(c&15);
return k?u:-u;
}
void wr(register int a){
if(a<0)Pc('-'),a=-a;
if(a<=9)Pc(a|'0');
else wr(a/10),Pc((a%10)|'0');
}
signed const INF(0x3f3f3f3f),NINF(0xc3c3c3c3);
long long const LINF(0x3f3f3f3f3f3f3f3fLL),LNINF(0xc3c3c3c3c3c3c3c3LL);
#define Ps Pc(' ')
#define Pe Pc('\n')
#define Frn0(i,a,b) for(register int i(a);i<(b);++i)
#define Frn1(i,a,b) for(register int i(a);i<=(b);++i)
#define Frn_(i,a,b) for(register int i(a);i>=(b);--i)
#define Mst(a,b) memset(a,b,sizeof(a))
#define File(a) freopen(a".in","r",stdin),freopen(a".out","w",stdout)
#define P (998244353)
#define G (3)
#define Gi (332748118)
#define N (2100000)
int n,m,l,s,r[N],a[N],b[N],o,w,x,y,siv;
int fpw(int a,int p){return p?a>>1?(p&1?a:1)*fpw(a*a%P,p>>1)%P:a:1;}
void ntt(int*a,bool iv);
signed main(){
Rd(n),Rd(m),siv=fpw(s=1<<(l=log2(n+m)+1),P-2);
Frn1(i,0,n)Rd(a[i]);
Frn1(i,0,m)Rd(b[i]);
Frn0(i,0,s)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
ntt(a,0),ntt(b,0);
Frn0(i,0,s)a[i]=a[i]*b[i]%P;
ntt(a,1);
Frn1(i,0,n+m)wr(a[i]),Ps;
exit(0);
}
void ntt(int*a,bool iv){
Frn0(i,0,s)if(i<r[i])swap(a[i],a[r[i]]);
for(int i(2),i2(1);i<=s;i2=i,i<<=1){
o=fpw(iv?Gi:G,(P-1)/i);
for(int j(0);j<s;j+=i){
w=1;
Frn0(k,0,i2){
x=a[j+k],y=w*a[j+k+i2]%P;
a[j+k]=(x+y)%P,a[j+k+i2]=(x-y+P)%P,w=w*o%P;
}
}
}
if(iv)Frn0(i,0,s)a[i]=a[i]*siv%P;
}
Time complexity: \(O(n\log n)\)
Memory complexity: \(O(n)\)
看看效果
時間上的提高效果不大，可是空間少了一半（由於用了int而不是complex）
Conclusion:
打了一天的博客終於寫完了（好累）
可是對FFT和NTT的理解也加深了很多
這個算法對數學知識和分治思想的要求都很高
本蒟蒻花了近一年的時間才真正理解
若是有錯誤和意見請大佬多多指教
那麼本篇博客就到這裏啦，謝謝各位大佬的支持！ありがとう！
- 1. 快速傅里葉變換FFT& 數論變換NTT
- 2. 快速傅里葉變換（FFT）與快速數論變換（NTT）+例題
- 3. 快速傅里葉變換與快速數論變換（一）
- 4. 快速傅里葉變換（FFT）和逆快速傅里葉變換（IFFT）
- 5. 快速傅里葉變換與快速數論變換
- 6. 初探FFT（快速傅里葉變換）
- 7. 快速傅里葉變換(FFT)詳解
- 8. 快速傅里葉變換
- 9. 快速數論變換NTT
- 10. 快速傅里葉變換（FFT）詳解
- 更多相關文章...
- • Eclipse 快速修復 - Eclipse 教程
- • ASP 快速參考 - ASP 教程
- • 使用阿里雲OSS+CDN部署前端頁面與加速靜態資源
- • IntelliJ IDEA 代碼格式化配置和快捷鍵
-
每一个你不满意的现在，都有一个你没有努力的曾经。