淺談FFT、NTT和MTT

前言

\(\text{FFT}\)（快速傅里葉變換）是 \(O(n\log n)\) 解決多項式乘法的一個算法，\(\text{NTT}\)（快速數論變換）則是在模域下的，而 \(\text{MTT}\)（毛神仙對\(\text{FFT}\)的精度優化算法）能夠針對任意模數。本文主要講解這三種算法，具體的應用還請參考我博客內的題解。

正文

FFT-快速傅里葉變換

學習這個算法能夠藉助《算法導論》，固然算導上的東西須要耐心才能啃下來。這裏只是歸納一下算導上的介紹，並加入一些我的的看法。下面逐步介紹這個算法。

複數

若是學過的話能夠跳過。實數能夠一一對應數軸上的點，那麼複數就能夠一一對應平面直角座標系上的點。對應 \(x\) 軸上的點的就是咱們熟悉的實數，而外面的就是虛數。其中 \((0,1)\) 這個點對應的數記做 \(i\) ，即 \(\sqrt{-1}\)，它表示虛數單位。複數能夠表示成 \(a+ib\) 的形式，其中 \(a,b\) 爲實數。

極座標表示法下的乘法

\((a,\alpha)\cdot(b,\beta)=(ab,\alpha+\beta)\)

證實以下：
\[ \begin{align}{} &(a,\alpha)\cdot(b,\beta) \notag\\ =&ab(\cos\alpha+i\sin\alpha)(\cos\beta+i\sin\beta)\notag\\ =&ab[\cos\alpha\cos\beta+i(\sin\alpha\cos\beta+\cos\alpha\sin\beta)-\sin\alpha\sin\beta]\notag\\ =&ab[\cos(\alpha+\beta))+i\sin(\alpha+\beta)]\notag\\ =&(ab,\alpha+\beta)\notag \end{align} \]

代數表示法下的乘法

\((a+ib)\cdot(c+id)=ac-bd+i(ad+bc)​\)

無需證實，肉眼化簡。

單位複數根

在單位圓上，咱們用 \(\omega_{n}^k\) 表示將單位圓 \(n\) 等分，取其第 \(k\) 條線對應的單位複數。其中 \(\omega_n^0=1\) ，逆時針方向編號，如圖所示：

單位復根有一些重要的性質。

消去引理

\(\omega_{dn}^{dk}=\omega_{n}^k\) 其中 \(n,k\geq 0,d>0\)

折半引理

\(\omega_n^{k+n/2}=(\omega_n^k)^2\) 其中\(n\geq0,k\geq 0\)

若是藉助向量去理解的話，理解起來很是方便。

多項式

一個形如 \(\displaystyle A(x)=\sum_{i=0}^{n-1}a_ix^i\) 的式子。

係數表示

直接列出 \(A(x)\) 的各項係數。這種表示方法能夠 \(O(n)\) 的實現多項式加法，但多項式乘法卻須要 \(O(n^2)\)

點值表示

經過帶入若干個特值肯定，顯然，一個最高次爲 \(n-1\) 的多項式須要 \(n\) 的特殊值便惟一肯定。這種表示方法能夠 \(O(n)\) 的加和乘，可是要轉化成係數表示才能體現出它做爲多項式的價值。

DFT

對於一個列向量 \(a=(a_0,a_1,\cdots,a_{n-1})\) ，以它爲係數的多項式 \(A(x)=\displaystyle\sum_{j=0}^{n-1}a_jx^j\)

如有一個列向量 \(y=(y_0,y_1,\cdots,y_{n-1})\) 知足 \(y_k=A(\omega_n^k)\) ，則\(y=\text{DFT}_n(a)\)

\(\text{DFT}\) 的全稱爲離散傅里葉變換，是將多項式的係數表達化做點值表達的一個變換。

同理 \(a=\text{DFT}_n^{-1}(y)\) ，\(\text{DFT}^{-1}\) 就是逆離散傅里葉變換，也稱 \(\text{IDFT}\)，咱們嘗試寫出 \(\text{DFT}^{-1}\) 的表達式。

寫出 \(y\) 與 \(a\) 的關係
\[ y_k=\sum_{j=0}^{n-1}a_j\omega^{kj} \]
而後咱們能夠矩陣乘積 \(y=V_na​\) 的形式表示向量 \(a​\) 到向量 \(y​\) 的變換。\(V_n​\) 爲由 \(\omega_n​\) 各項指數構成的範德蒙德矩陣。
\[ \begin{pmatrix} y_0\\ y_1\\ y_2\\ y_3\\ \vdots\\ y_{n-1}\\ \end{pmatrix} = \begin{pmatrix} \omega^0 & \omega^0 &\omega^0 & \omega^0 & \cdots & \omega^0 \\ \omega^0 & \omega^1 &\omega^2 & \omega^3 & \cdots & \omega^{(n-1)}\\ \omega^0 & \omega^2 &\omega^4 & \omega^6 & \cdots & \omega^{2(n-1)} \\ \omega^0 & \omega^3 &\omega^6 & \omega^9 & \cdots & \omega^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots &\vdots\\ \omega^{0} & \omega^{1(n-1)} &\omega^{2(n-1)} & \omega^{3(n-1)} & \cdots & \omega^{(n-1)(n-1)} \\ \end{pmatrix} \begin{pmatrix} a_0\\ a_1\\ a_2\\ a_3\\ \vdots\\ a_{n-1}\\ \end{pmatrix} \]
那咱們如今要求的就是 \(V_n^{-1}\) 的矩陣，即 \(V_n\) 的逆矩陣。

有以下定理：

對於 \(j,k\in[0,n)\) ，\(V_n^{-1}\) \((j,k)\) 處的元素爲 \(\omega_n^{-jk}/n\)

證實以下
\[ \begin{array}{} [V_nV_n^{-1}]_{jj'}&=\displaystyle\sum_{k=0}\omega_n^{jk}\omega^{-kj'}/n\\ &=\displaystyle{1\over n}\sum_{k=0}\omega_n^{k(j-j')} \end{array} \]
顯然，當 \(j=j'\) 時，\([V_nV_n^{-1}]_{jj'}\) 的值爲 \(1\) ，不然爲 \(0\) ，那麼 \([V_nV_n^{-1}]\) 是一個行列數爲 \(n\) 的單位矩陣，即得證 \(V^{-1}\) 爲 \(V\) 的逆矩陣。

那麼在做 \(\text{IDFT}​\) 的時候，只需將單位根換成 \({\omega_n^{-1}}​\) ，最後係數再除以 \(n​\) 便可。

固然，直接變換是 \(O(n^2)\) 的。咱們考慮用分治的思想進行變換。

FFT

首先觀察多項式 \(A(x)\) ，咱們將指數分奇偶兩類。偶數項以 \(\{a_0,a_2,...,a_{n-2}\}\) 構造一個新的多項式 \(\displaystyle A^{[0]}(x)=\sum_{j=0}^{n/2-1}a_{2j}x^j\)，奇數項同理爲 \(\displaystyle A^{[1]}(x)=\sum_{j=0}^{n/2-1}a_{2j+1}x^j\)。

那麼顯然有
\[ A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2) \]
咱們把 \(\omega_n^k\) 代入獲得
\[ A(\omega_n^k)=A^{[0]}(\omega_n^{2k})+\omega_n^kA^{[1]}(\omega_n^{2k}) \]
利用消去引理獲得
\[ A(\omega_n^k)=A^{[0]}(\omega_{n/2}^k)+\omega_n^kA^{[1]}(\omega_{n/2}^k) \]
那麼將 \(A^{[0]},A^{[1]}\) 的係數向量 \(a^{[0]},a^{[1]}\) 進行一次 \(\text{DFT}\) ，分別獲得 \(y^{[0]},y^{[1]}\) 。

有
\[ y^{[0]}_k=A^{[0]}(\omega_{n/2}^k)\\ y^{[1]}_k=A^{[1]}(\omega_{n/2}^k) \]
只要令 \(k<n/2​\) ，將 \(k\geq n/2​\) 的部分用折半引理便可。
\[ A(\omega_n^k)=A^{[0]}(\omega_{n/2}^k)+\omega_n^kA^{[1]}(\omega_{n/2}^k)\\ A(\omega_n^{k+n/2})=A^{[0]}(\omega_{n/2}^k)-\omega_n^kA^{[1]}(\omega_{n/2}^k) \]
推導不難，注意將在單位圓上的旋轉借用平面向量來理解。

用 \(y\) 代入，最終的表達式爲
\[ y_k=y^{[0]}_k+\omega_n^ky_k^{[1]}\\ y_{k+n/2}=y_k^{[0]}-\omega_n^ky_k^{[1]} \]
這樣就能夠分治求解了。

更高效的FFT

事實上 \(\text{FFT}\) 能夠迭代求解。先觀察一下遞歸求解的過程，如圖所示。

而後用人類智慧觀察，發現 \(a_i​\) 在底層是在的位置爲 \(i​\) 的二進制位翻轉。

發現只須要枚舉區間長度，掃整個序列，就能夠進行對區間進行合併。觀察遞歸求解的式子
\[ y_k=y^{[0]}_k+\omega_n^ky_k^{[1]}\\ y_{k+n/2}=y_k^{[0]}-\omega_n^ky_k^{[1]} \]

它的流程能夠用上圖來表示，上面操做叫做蝴蝶操做，其實和遞歸求解的流程類似。具體仍是看代碼，碼風仍是清晰的。

struct Complex
{
    double x,y;
    Complex operator +(const Complex &_){return (Complex){x+_.x,y+_.y};}
    Complex operator -(const Complex &_){return (Complex){x-_.x,y-_.y};}
    Complex operator *(const Complex &_){return (Complex){x*_.x-y*_.y,x*_.y+y*_.x};}
    Complex operator /(const int &_){return (Complex){x/_,y/_};}
};
namespace _Polynomial
{
    Complex A[N<<1],B[N<<1];
    Complex w[N<<1];int r[N<<1];
    void DFT(Complex *a,int op,int n)
    {
        FOR(i,0,n-1)if(i<r[i])swap(a[i],a[r[i]]);   //位翻轉
        for(int i=2;i<=n;i<<=1)     //合併出一個長i的區間
            for(int j=0;j<n;j+=i)       //區間開頭的位置
                for(int k=0;k<i/2;k++)      //蝴蝶操做
                {
                    Complex u=a[j+k],t=w[op==1?n/i*k:n-n/i*k]*a[j+k+i/2];
                    a[j+k]=u+t,a[j+k+i/2]=u-t;
                }
        if(op==-1)FOR(i,0,n-1)a[i]=a[i]/n;
    }
    void multiply(const int *a,const int *b,int *c,int n1,int n2)
    {
        int n=1;
        while(n<n1+n2-1)n<<=1;
        FOR(i,0,n1-1)A[i].x=a[i],A[i].y=0;
        FOR(i,0,n2-1)B[i].x=b[i],B[i].y=0;
        FOR(i,n1,n-1)A[i].x=A[i].y=0;
        FOR(i,n2,n-1)B[i].x=B[i].y=0;
        FOR(i,0,n-1)r[i]=(r[i>>1]>>1)|((i&1)*(n>>1));
        FOR(i,0,n)w[i]=(Complex){cos(2*PI*i/n),sin(2*PI*i/n)};
         
        DFT(A,1,n),DFT(B,1,n);
        FOR(i,0,n-1)A[i]=A[i]*B[i];
        DFT(A,-1,n);
        FOR(i,0,n1+n2-2)c[i]=A[i].x+0.5;
    }
};

顯而易見，因爲 \(\text{double}\) 的存在，精度多多少少會被卡一點。而具體的題目常常每每會給一個特殊的模數，這種時候就要用到接下來介紹的算法了。

NTT-快速數論變換

待補充。。。