1 Forewords
卷積,但不止卷積 - FFT 漫談
- 先有 FT,再有 DFT,才有 FFT
- 時頻轉換是最初的用途
- 發現單位根優秀性質,James Cooley, John Tukey 發明現代 FFT 加速 DFT,但此前相似的發現早已有之
- 後來將 DFT 與卷積定理聯絡,FFT 才被用於計算多項式乘法
- 複數運算精度誤差推動了 NTT 的發展
- 應用:任何需要頻率和卷積的地方.頻譜、濾波器、音樂、雷達、影像處理……
- OI/XCPC 中主要關心卷積
推薦食用方法
初步要求
- 知道 DFT、FFT 可用於快速計算多項式卷積
- 掌握 FFT 加速 DFT 計算的原理和實現
- 會應用結論改動 FFT 加速 NTT 計算
- 見識一些卷積解決的基本問題,初步瞭解生成函式在組合計數中的應用
- 題目可選擇性完成,請多花時間消化原理和思想
學有餘力 / 集訓後繼續消化
- 系統學習生成函式
- 實現多項式全家桶
- 對原理感興趣的同學可對數學部分做進一步研究.本講內容是線性代數、抽象代數、數論等多領域的綜合應用.歡迎討論.
- 學習集合冪級數相關知識點(FMT,FWT,……),體會其思想與 FFT 的同與異
- 學習 FFT 在訊號、頻譜等非演算法競賽向實際問題中的應用
Learn for fun :)
記號說明
記
,此時可用 代替下標取值範圍 的記法.- 集合論中已定義
,這裡的中括號是為了強調其集合含義.
- 集合論中已定義
使用 Iverson 括號約定:設
是一個命題,記多項式的規模定義為多項式的次數加一.特別的,零多項式的規模為
.- 以後會混用
次多項式和規模為 的多項式的說法.
- 以後會混用
2 FFT/NTT in a nutshell
2.1 FFT
多項式卷積
給定兩個至多
係數 - 點值 - 係數
代入任意
可得到多項式在 處的點值點值意義下的多項式乘法是
的 點確定一個至多 次的多項式1
若計算至多
係數 - 點值 - 係數 - 快速轉換?
樸素計算任意指定
個位置點值需要 .Lagrange[1] 插值給出了
將任意位置 個點值還原為多項式係數的演算法.能否選取
個特殊的點值使係數 - 點值、點值 - 係數的變換支援快速計算?
Discrete Fourier Transform
離散傅立葉變換(Discrete Fourier Transform, DFT)接受一個至多
得益於單位根的特殊運算性質,二者均有被稱為快速傅立葉變換(Fast Fourier Transform, FFT)的快速演算法.
複數域單位根
複數域上的
所有單位根模長均為
複數域單位根 - 三個重要性質3
Theorem 1 (消去引理)
Theorem 2 (折半引理)
消去 / 折半引理將在 FFT 的推導中使用.
複數域單位根 - 三個重要性質
Theorem 3 (求和引理)
求和引理的證明使用了等比數列求和公式.將在 IDFT 的推導中用到.
Fast Fourier Transform
考慮將至多
DFT 的矩陣表示
記
IDFT
由單位根的消去引理可證,DFT 矩陣
故快速計算 IDFT 的方法與 FFT 幾乎一致,只需將計算 DFT 時使用的本原單位根
梳理
怎麼計算卷積?
- 把至多
次的多項式 和至多 次的多項式 寫成至多 次的多項式(高位補 ).為計算 FFT 方便,還要繼續補 至一個大於其次數的 的冪. - 對
和 多點求值. - 把兩個多項式的點值逐點相乘.
- 多點插值還原
的係數.
- 把至多
怎麼快速求值?
- 選點選單位根就是 DFT.
- 消去引理和折半引理使我們可以用 FFT 演算法遞迴地計算 DFT.
- 推導已經給出了遞迴的寫法,之後還會介紹常數更優的迭代實現.
怎麼快速插值?
- 求和引理給出了 DFT 矩陣的逆矩陣.
- 計算方法很相似,最後逐項除掉一個規模.
FFT 遞迴實現 - DFT 部分
方便起見,我們只處理
#include<bits/stdc++.h>
#include<complex>
using namespace std;
typedef long long ll; typedef complex<double> CP;
const ll MXN=4E6+5; const double PI=3.14159265358979323846l;
[MXN];
CP tmpvoid _DFT(CP A[],ll n,ll typ){
/=2; if(n==0) return;
nfor(ll k=0;k<n;k++) tmp[k]=A[2*k],tmp[n+k]=A[2*k+1];
for(ll k=0;k<2*n;k++) A[k]=tmp[k];
(A,n,typ); _DFT(A+n,n,typ);
_DFT(cos(2*PI/(2*n)),typ*sin(2*PI/(2*n))), wk=1;
CP wfor(ll k=0;k<n;k++){
[ k]=A[k]+wk*A[n+k];
tmp[n+k]=A[k]-wk*A[n+k];
tmp*=w;
wk} for(ll k=0;k<2*n;k++) A[k]=tmp[k];
} void DFT(CP A[],ll n,ll typ){
(A,n,typ); if(typ==-1) for(ll i=0;i<n;i++) A[i]*=1.0l/n;
_DFT}
FFT 遞迴實現 - 卷積部分
// alternatively, use std::__lg() in GCC
(ll n){ll cnt=0; for(ll t=1;t<n;t<<=1) cnt++; return cnt;}
ll log2ceil[MXN],B[MXN],C[MXN]; ll outC[MXN];
CP A* conv(ll inA[],ll aN,ll inB[],ll bN){
ll=1LL<<log2ceil(aN+bN-1);
ll nfor(ll i=0;i<aN;i++) A[i]=inA[i];
for(ll i=0;i<bN;i++) B[i]=inB[i];
(A,n,1); DFT(B,n,1);
DFTfor(ll i=0;i<n;i++) C[i]=A[i]*B[i];
(C,n,-1); for(ll i=0;i<n;i++) outC[i]=round(C[i].real());
DFTreturn outC;
}
Drawbacks?
- 遞迴實現慢
- 臨時陣列醜
- 封裝性為零
FFT 迭代
迭代地實現 FFT 不僅在常數上更加優秀,亦更便於使用 C++ 的容器進行封裝.這並不困難,只需自底向上模擬 FFT 遞迴過程即可.
唯一的問題——最底層的順序?
來觀察一輪
你發現了什麼?
蝶形運算
在計算點值前,
我們有
void spawnrev(ll n){
[0]=0;
revfor(ll i=1;i<(1<<n);i++)
[i]=(rev[i>>1]>>1)+((i&1)<<(n-1));
rev}
FFT 迭代實現
void DFT(CP A[],ll n,ll typ){ // rev[] should be spawned in advance
for(ll i=0;i<n;i++) if(i<rev[i]) swap(A[i],A[rev[i]]); // a one-to-one permutation
for(ll hf=1;hf<n;hf*=2){
(cos(2*PI/(2*n)),typ*sin(2*PI/(2*n))), wk=1;
CP wfor(ll i=0;i<n;i+=hf*2){
=1;
CP wkfor(ll k=0;k<hf;k++){
=A[i+k],y=wk*A[i+hf+k];
CP x[i+k]=x+y; A[i+hf+k]=x-y;
A=wk*w;
wk}
}
}
if(typ==-1) for(ll i=0;i<n;i++) A[i]*=1.0l/n;
}
請自行實現更易用的容器封裝版本.
2.2 NTT
NTT 速成5
FFT 的缺點?浮點數帶來的大常數與精度問題.
我們指出,在係數和點值模
我們指出,對於滿足
只需修改單位根定義,把複數運算改為整數取模,就得到了能算 NTT 的 FFT 的實現.
const PR=3,MOD=998244353;
=qpow(PR,(MOD-1)/(hf*2)); if(typ==-1) w=inv(w); ll w
FFT/NTT in a nutshell - 小結:概念區分
關於 DFT
- Discrete Fourier Transform, DFT, 離散傅立葉變換
- Fast Fourier Transform, FFT, 快速傅立葉變換
- FFT 是計算 DFT 的快速演算法
關於 NTT
- Number Theoretic Transform, NTT, 數論變換
- FFT 在複數域上的多項式環
中進行,而 NTT 在模 剩餘類域 上進行 - 快速計算 DFT / IDFT, NTT 都用 FFT,故一般不使用”FNTT”的說法
3 Applications
3.1 基本應用
基本應用
- 大整數乘法
-
十進位制數可拆解為多項式表示,計算卷積後處理進位即可.由於數字最大隻是
,合理資料範圍下捲起來不會爆模數,直接用 NTT 實現即可. - 揹包計數
-
兩個揹包的合併就是多項式卷積.
- 濾波器
-
反轉多項式的係數陣列再做卷積,可以快速得到兩個多項式滑動視窗式的內積.
- 位運算
-
有些位運算可以寫成卷積的形式.模
意義下”異或”是加法,“與”是乘法,“或”可以透過取反轉化為”與”. - 字串
-
透過巧妙設定字串距離函式,FFT 可解決更廣泛的字串匹配問題.
- 值域次數化
-
當值域較小時,將待計算的值放在多項式次數上統計貢獻次數,可以繞開某些極難求解的數值問題.例如 Vandermonde 行列式的快速計算.
基本應用 - 分治 FFT
對多個長度相同的多項式的卷積,分治地捲起來可降低時間複雜度.長度不一時,挑小的先卷也可減小常數(用堆維護).
另有一種 CDQ 風格的分治 FFT.CDQ 長於處理帶偏序的二元點對貢獻,在處理形如
當卷積的前後項存在依賴關係時,也可使用此法保證處理順序恰當.然而此類依賴問題往往也可透過解生成函式方程的方法求得封閉形式.
時間複雜度均為
3.2 生成函式初步
生成函式初步 - 導言
生成函式是一種對數列的操作技巧.透過將數列表示為多項式或形式冪級數,數列間複雜的和式操作可用簡單函式的乘法、複合等運算進行表示,從而大大降低了數列變換技巧的使用門檻.
生成函式在組合數學中應用廣泛,且生成函式的部分操作在組合意義下也有較為直觀的理解.本節將帶大家初窺其中的奧妙.限於篇幅和主講人能力,我們僅以題帶點地講解,期冀為大家建立構造生成函式的直覺.請感興趣的同學下來做進一步研究.
熟悉 Taylor 展開的同學或能較快上手此部分內容.
Ordinary Generating Function
序列
Exercise 1 寫出下列數列的 OGF.下標從
OGF 組合意義
OGF 相乘,是揹包,是卷積,是兩塊無標號組合物件的有序拼接.
Exercise 2 寫出下列計數問題的 OGF,均以
- 在
個物品中選出 個的方案數. - 容量為
的揹包裝下體積分別為 的 個物品的方案數. - 同上,但每個物品有無限個.
個無標號球放入 個有標號盒的方案數,要求盒非空.- 同上,但盒可空.
- 將整數
分拆為若干正整數之和的方案數.
Exponential Generating Function
序列
Exercise 3 寫出下列序列的 EGF,下標從
EGF 組合意義
觀察兩個 EGF 的乘積
EGF 組合意義
Exercise 4 寫出下列計數問題的 EGF,均以
- 長度為
的排列的構型數. - 長度為
的圓排列的構型數. - 將
個有標號球放入 個有標號盒的方案數,要求盒非空. - 將
個有標號球放入 個無標號盒的方案數,要求盒非空. - 將
元集合劃分為 個等價類的方案數. - 劃分
元集合的方案數.
4 Mathematics behind
4.1 NTT 原理
NTT 原理 - 導言
本節介紹 NTT 的原理.
FFT 加速卷積演算法的核心,一是多項式的求值插值原理,二是單位根帶來的分治快速演算法.我們將在本節中證明,模
NTT 原理涉及原根等數論內容.本講的目標是建立 DFT 變換和 FFT 演算法的通用數學框架,而非具體研究其某一特例.故我們只講解 NTT 所需的基礎數論知識,無關的細節則略過處理.對數論感興趣的同學可前往 OI Wiki 學習.
模 剩餘類域
模
模
模
上的多點插值
之後記係數均在
Lemma 1 (
上的多點插值
Theorem 4 (Lagrange 定理) 設
Corollary 1 設
推論告訴我們,欲確定
本原單位根
回顧 DFT 中複數域
事實上所有的
我們把這一類重要的單位根稱為本原單位根.抽象的來說,
上的階
如何定義
稱在模
- 階最高有多高?
Theorem 5 (Fermat 小定理) 若
Theorem 6 (Euler 定理) 若
Euler 定理的證明 - 簡化剩餘系
簡化剩餘系的大小即
設
的取值 - , 互質時
補充討論
當
的取值 - 一般情況9
當
原根10
需要注意的是,Euler 定理只給出了
稱
Theorem 7 (原根存在定理)
求階和原根
- 求階,用定理
,求因子 + 快速冪即可 . - 找原根,從小到大用上述求階方法暴力即可.最小原根一般不會太大.
- 更快的方法請參考洛谷求原根模板題題解.
下面再介紹原根的兩個定理.
Theorem 8 (原根判定定理) 若
必要性顯然.充分性,反證出所有
Theorem 9 (原根個數定理) 若
上的本原單位根構造
設
使用數論中 Bézout 定理,
遺憾的是,一般的
上的 NTT - 是質數
求和引理是保障 NTT 逆變換對應矩陣確為
為保證
綜上,逆變換矩陣確為
上的 FFT 與卷積加速 - 998244353
為滿足 FFT 對
int
型中的單次加減操作不會溢位,是 OI/XCPC 計數題中不可多得的優秀模數12.
結合前述關於
上的 NTT - 對一般的
當
我們找到一篇有關該問題的參考文獻[7, Secs. 3 and appendix B],但尚不確定其證明的正確性.友情提示讀者:
從分析到代數
後續數學內容導讀 NTT 原理雖已非常”數學”,但也只是 DFT 在有限域上的一個例項.本節往後,我們要嘗試為多項式係數位於複數域
後續數學內容不再要求掌握.望同學們在紛繁的定理定義中抓住要旨,窺見抽象數學背後蘊藏的規律.熟悉高等代數和抽象代數的同學或會對某些內容感到熟悉.抽象地討論 FFT 的資料並不多見,後續內容多為主講人的新進探索,或有謬誤,望不吝指正.
4.2 求值與插值
求值與插值 - 導言
本節將重新審視已經熟知的多項式,把抽象的、代數的多項式和具體的、分析的多項式函式區分開來.我們指出,多項式和多項式函式不同但關聯緊密,形式冪級數與冪級數亦有此類聯絡.這些抽象的討論將幫助我們剖析多項式求值插值的基本原理.
除常見代數書目(如[8]),也推薦參考 OI Wiki 的多項式基礎[9]和 Wikipedia 的形式冪級數[10].
多項式
設(無窮)序列
多項式環上的加法、乘法的定義已經為大家所熟知.係數所處的環保證了多項式加法和乘法的良定義,而在這兩種運算下,
Remark (群、環、域). 群、環、域都是常見的代數結構,其中的元素在給定運算下封閉,並滿足特定的運算性質.簡單來說,環13上定義了加法和可能不可逆的、不一定交換的乘法,域上定義了加減乘除所有四則運算.交換環中的乘法滿足交換律.除環中的所有元素有乘法逆元.域是交換除環.
多項式
習慣上也會將多項式
多項式
下面額外為多項式定義一種新的運算.多項式
形式冪級數
形式冪級數定義與多項式的唯一區別是其不要求
由於涉及無限次運算,形式冪級數的複合運算需考慮環
Remark. DFT, NTT 與多項式環 DFT/FFT 加速的多項式卷積在複數域
帶餘除法
整環是無零因子的交換么環.所謂無零因子,即環中任意元素
定義有帶餘除法的環被稱為 Euclid 整環.域上的多項式環都是 Euclid 整環[8, 第 7 章第 2 節定理 3, p. 11].值得注意的是,
如在帶餘除法中保證除式是首一多項式17,則帶餘除法的良定義和進行過程也均可在整環上實現.
多項式函式 / 冪級數
剛剛強調,多項式 / 形式冪級數只是定義了加法和乘法的序列.現在介紹多項式函式和冪級數.它們不是序列,而是對映18.
多項式
與多項式 / 形式冪級數不同,這裡的
多項式函式的加法和乘法定義為函式的加法和乘法,即
求值
多項式和多項式函式似乎在許多情況下有著平行的關係.下面介紹一個較直觀的結論.
若環
多項式和多項式函式的這一關係為多項式在任意點的求值操作提供了理論基礎.在多項式函式的
插值
需要注意的是,前述結論的逆命題不一定成立,即環
該逆命題的本質是透過多項式函式的所有函式值反過來確定多項式(係數)的過程.而如果在確定時只使用一部分函式值,就是所謂的多點插值過程.很多時候待求多項式的次數是已知的,這在相當程度上縮小了待定多項式的範圍.我們指出,只要
證明的關鍵是用帶餘除法24討論多項式根與其一次因式的關係,即多項式餘式定理或小 Bézout 定理[8, 第 7 章第 6 節定理 6, p. 35].
需要強調,上述結論只能說明,只有那些確實可透過多項式
求值與插值的線性表示
到這裡,多項式的求值與插值的線性表示已經呼之欲出了.若將交換環
若將對
求值與插值 - 小結
本節的核心是多項式和多項式函式的區別與聯絡,兩個方向的”確定”分別給出了多項式多點求值和多點插值的理論基礎.
由於多項式和多項式函式這種若即若離的關係,往往在記號上也有意無意地混淆了它們,某些情況下加大了區分的難度.本篇使用的記號體系將盡量用單個字母
4.3 環上的 DFT
環上的 DFT - 導言
前面討論了在多項式任意點處求值插值的基本原理,但 DFT/FFT 的執行只需在單位根處求值和插值.本節將進一步放寬對多項式環的限制,介紹定義在有主要單位根的環上的一般的 DFT 及其快速演算法 FFT.
本節內容主要參考了 [12, Secs. 2, pp. 983–984] 和 [13].
環上的單位根 - 兩個定義
定義環
定義環
本原單位根在許多情況下與主要單位根等價,但亦非完全相同.
環上的單位根 - 區別與聯絡
Proposition 1 若環
Proof. 對任意正整數
環上的單位根 - 區別與聯絡
Proposition 2 若環
Proof. 反證.若存在一正整數
環上的單位根 - 其它性質
分別根據定義和數論中的 Bézout 定理,容易證明主要單位根和本原單位根的消去引理:設
注意到
環上的
環上的 DFT
設
環上的 FFT
除無法在任意環上使用
設
考慮對模
環上的 FFT
注意到
因此,令
4.4 迴圈卷積與卷積定理
迴圈卷積與卷積定理 - 導言
上一節中,我們建立了在有主要單位根的環上的 DFT 及其快速演算法 FFT 的相關理論,但由於放寬了對多項式環
在求值與插值部分已經介紹求值變換
迴圈卷積
設
設
令
迴圈卷積
卷積定理
前述討論已經證明29
我們指出,當要求所作變換可逆時,卷積定理反過來也要求所作變換是一類似 DFT 對映的變換30.
Exercise 5 嘗試將任意序列 DFT 兩次,觀察結果.證明你的結論.
Exercise 6 求 DFT 矩陣的行列式.儘可能縮小可行解範圍.
5 Advanced Operations
多項式全家桶 - 序言
有哪些?
- 求逆、開根、對數、指數、快速冪、複合等
咋推的?
- 求解思路幾乎都是倍增,時間複雜度幾乎都是大常數
- Newton 迭代法是推導全家桶迭代公式的通法
- 嚴格化需要進一步的形式冪級數理論,主講人不會
- 求解思路幾乎都是倍增,時間複雜度幾乎都是大常數
有啥用?
- 當你一波操作化出生成函式發現不會求係數
怎麼講?
- 受限於篇幅和主講人能力,我們講不完
- 只講求逆和對數,其餘請左轉 OI Wiki
多項式求逆
給定一多項式
多項式逆元存在的充分必要條件是其常數項非零(這是因為邊界條件
多項式求逆
以下簡記
兩邊平方得
多項式求逆 - 實現
(const Poly &A){
Poly inv=A.len(); Poly B(1); B[0]=inv(A[0]);
ll nfor(ll hf=1;hf<n;hf<<=1){
=B*2-B*B*A.subpoly(0,hf*2); B.resize(hf*2);
B} B.resize(n);
return B;
}
- 常數巨大的寫法,僅作演示,請勿學習.
- 正確的寫法是在 DFT 後的點值上操作多項式,請小心實現封裝.
多項式
給定一多項式
次數為
推導是容易的.方程兩側同時求導得
多項式 - 實現
(Poly A){ // derivative
Poly drvfor(ll i=0;i<A.len();i++) A[i]=(i+1)*A[i+1]_;
.pop_back(); return A;
A}
(Poly A,ll c){ // integral
Poly itg.push_back(0); for(ll i=A.len();i>=1;i--) A[i]=A[i-1]*inv(i)_;
A[0]=c; return A;
A}
(const Poly &A){
Poly lnreturn itg((drv(A)*inv(A)).subpoly(0,A.len()-1),0/*log(A[0])*/);
}
Acknowledgements
感謝
keke_046
學長教授 FFT、集合冪級數與生成函式.微言大義,博大精深,至今仍在消化.感謝隊友
ItzDesert
提供位運算典題一道並提供內容編排建議.
題單
主講人練題少,僅供參考.
- 模板
- 大整數乘法
-
- 洛谷 P1919 【模板】A*B Problem 升級版(FFT 快速傅立葉變換)
- 基礎計數
-
SPOJ-TSUM Triple Sums
BZOJ3513-MUTC2013 Idiots
上面兩道題都是 OGF 消序,較 EGF 消序困難.一般的方法是使用 Polya 計數原理.
- 位運算
-
ABC291G OR Sum
也是濾波器的應用.
- 分治 FFT
-
百度之星 2023 初賽第二場 T8
容斥後需要計算若干一次多項式乘積,分治 NTT 即可.
-
求
.值域次數化後 CDQ 處理偏序. 洛谷 P4721 【模板】分治 FFT
CDQ 偏序化處理前後項依賴.也可解生成函式方程再多項式求逆.
- 字串
- 全家桶
- 其它
References
Footnotes
一種常見的證法是使用 Vandermonde 行列式證明矩陣可逆.後面會介紹多項式環風格的證明.↩︎
有的文獻定義
,或是因為訊號處理領域常用 IDFT 將訊號時域取樣資料變為頻域資訊.事實上,DFT/IDFT 的說法也常有反轉,但這只是形式問題.↩︎這也表明適當歸一化後的 DFT 矩陣是一個酉矩陣.↩︎
NTT 原理需較多筆墨,稍後介紹.↩︎
否則只有
個不同元素的 中根本取不到 個不同位置的點值.後面會深入討論.↩︎多項式求逆等多項式進階操作,我們後續講解.↩︎
對質數
, .故 Euler 定理是 Fermat 小定理的一個推廣.↩︎這是後文所述定理“整環上的本原單位根也是主要單位根”在一般環上的一個反例.↩︎
網傳此模數由 UOJ 站長 vfleaking 提出並推廣.在所有需要取模的題目中使用該模數,可使選手無法透過模數判斷題目的做法.↩︎
本篇中環的定義包含乘法單位元,即么環.↩︎
無零因子的交換么環,稍後介紹.↩︎
可用前述多項式乘積次數公式證明.↩︎
試試在
上用 去除 .↩︎首項為
的多項式.↩︎函式和對映幾乎是等價名詞.有時函式特指值域包含於複數域
的對映.↩︎為良定義
,環 必須有單位元.↩︎這裡再次涉及環
上的收斂問題.由於實踐中只關心形式冪級數的前有限項,後續討論係數-點值-係數法轉換卷積時不需要用到冪級數理論,可以迴避.↩︎這種保持結構不變的對映被稱為同態(homomorphism).↩︎
其證明了域上一元多項式環的通用性質.仿照該證明應可證明環上的版本,從而證明這一同態關係.↩︎
該定理是下方高亮定理的一個自然的推論.↩︎
由於一次因式均為首一多項式,可以在整環上對其使用帶餘除法.↩︎
模是定義在環上的“線性空間”.↩︎
Lagrange 插值的構造用到了除法,且行列式非零推出矩陣可逆僅在域上的線性空間中適用,因此必須要求
是域.↩︎使得
的最小正整數 .不存在則記為 . , .可以證明域的特徵一定是 或一質數[8, 第 7 章第 11 節定理 3,p. 70].↩︎證明使用主要單位根的定義(求和引理)即可.↩︎
雖然前面用到了
可逆的要求,但該定理在 不可逆時也成立.只需類似地驗證兩邊相等即可.↩︎具體來說,該變換隻能是 DFT 矩陣的某個行置換.在
上的證明可參見[14],主講人目前正在研究整環上的版本,歡迎討論.↩︎