快速傅里葉變換|快速數論變換

FFT——快速傅里葉變換

什么是FFT

\(FFT\)全稱(\(Fast \ Fourier \ Transformation\))，中文名：快速傅里葉離散變換。

這個名字聽起來很高級，實際上也很高級，它是\(DFT\)和\(IDFT\)的加速版本，用於加速多項式乘法。

• 接下來我們把某個多項式記為“大寫字母+元素”，如多項式：

\(A(x)=a_nx^n+a_{n-1}x^{n-1}+…+a_2x^2+a_1x^1+a_0\)

\(B(x)=b_mx^m+b_{m-1}x^{m-1}+…+b_2x^2+b_1x^1+b_0\)

你也可以理解為\(A/B\)是關於\(x\)的函數

現在，我們要求\(A(x)*B(x)\)，即兩個多項式的乘積(卷積)，那么我們很容易想到的就是\(O(N^2)\)的做法：把每一項依次分別相乘再相加。

\[\sum^{n+m-1}_{i=0}(\sum^{i}_{j=0}a_i*b_{i-j})*x^i \]

這樣的時間復雜度明顯不盡人意，那么\(FFT\)就是把這種算式計算優化到\(NlogN\)的算法。

---

前置知識

• 復數

我們定義虛數 \(i\)， 且 \(i^2 = -1\)，即 \(i = \sqrt {-1}\)，那么復數就形如 \(a+bi\)，\(i\) 是一個虛數單位，\(a\) 叫做復數的實部， \(b\) 叫做復數的虛部，那么 \(a+bi\) 在復平面（高中會講）可以類似於平面直角坐標系，表示為點 \((a, b)\)。

復數的運算法則

1. 加法：實部和虛部分別相加。

2. 減法：實部和虛部分別相減。

3. 乘法：像一次多項式一樣相乘。

4. 除法：分子分母同乘分母的共軛。

共軛：兩個復數實部相同，虛部相反。

上面除法分子分母同乘分母的共軛，把分母有理化成實數。

    struct complex { //復數 
        long double a, b;
        complex() { a = b = 0;}
        //重載復數四則運算 
        complex operator + (const complex x) const {
            complex temp;
            temp.a = a + x.a;
            temp.b = b + x.b;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.a = a - x.a;
            temp.b = b - x.b;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.a = a * x.a - b * x.b;
            temp.b = a * x.b + b * x.a;
            return temp; 
        }

        complex operator / (const complex x) const {
            double div = x.a * x.a + x.b * x.b;
            complex temp;
            temp.a = (a * x.a + b * x.b) / div;
            temp.b = (b * x.a - a * x.b) / div;
            return temp; 
        }
    };

這邊大家直接記一個重要結論：

復數相乘，輻角相加，模長相乘。

---

單位根

對於復數\(w\),若整數 \(n\) 滿足 \(w^n=1\)，則稱 \(w\) 是 \(n\) 次單位根。

扔個單位圓：

為什么扔個單位圓呢？因為單位圓上任意一點到原點的距離，即模長都為\(1\)。

那么，對於一個復數，只有模長為\(1\)的復數才有可能成為\(n\)次單位根!

現在，我們在復平面上把復數鎖定在了這個單位圓上。

設復數\(w\)模長為\(q\)，其幅角為\(mπ\)，\(m\)為\([0,2]\)間實數，那么，其\(n\)次方幅角為\(nmπ\)。

則有\(nmπ=kπ,k∈N\)，即\(nm=k\)。

1. 若\(k=0\),因為\(n>0\)，所以方程有且只有\(m=0\)一解。

2. 若\(k>0\)，則\(k=1,2,3,…\)，易知\(m\)有\(n-1\)個不重解

綜上，復數\(w\)的\(n\)次單位根有\(n\)個。

並且，復數的\(n\)次單位根\(n\)等分單位圓。

單位根的書寫：

單位根用符號\(ω\)表示，從1開始，沿單位圓逆時針把單位根從\(0\)開始標號。

然后單位根還有個神奇的東西\(ω^{k}_{n}=ω_{n}^{k\mod n}\)

單位根的性質

1. 對於任意\(n\)，\(ω_n^0=1\)。

2. \(ω_n^k=(ω_n^1)^k\)。

3. \((ω_n^m)^k=(ω_n^k)^m\)

4. \(ω_n^i*ω_n^j=ω_{n}^{i+j}\)。

5. \(ω_{2n}^{2m}=ω_n^m\)。

6. 若 \(2 \mid n\)，\(ω_n^{m+\frac{n}{2}}=-ω_n^m\)

---

多項式的表達方法

1. 系數表達法：這個如同我們上面寫的多項式即為系數表達法。

2. 點值表達法

定義多項式\(F(x)\)：

我們知道，\(1\)個n元一次方程最少需要\(n\)個式子求解，而我們把系數表達\(F(x)\)抽象為一個函數\(F(x)\)，那么這個函數至少需要\(n+1\)個點確定下來，那么我們任取\(n+1\)個不同的數構成集合\(U={k_1,k_2,…,k_{n+1}}\)，對\(F(x)\)分別求值得到一個新的集合，那么將兩個集合合並成點集，有：

\[T={(k_1,F(k_1)),k_2,F(k_2),…,(k_{n+1},F(n+1))} \]

這便是\(F(x)\)的點值表達式，並且，點值表達法下的乘法運算更加簡單：

多項式\(P,Q\)分別取點\((x,y_1)\)和\((x,y_2)\)，則\(P*Q\)可得到點\((x,y_1*y_2)\)。令\(S=P*Q\)，則\(C(x_i)=y_{1_i}*y_{2_i}\)。

可見，點值表達法下，多項式乘法是\(O(N)\)的。

• （補充）乘法本質：卷積

根據我們上面暴力計算\(A(x)*B(x)\)，不妨設\(C=A*B\)，那么：

\(C[k]=\sum^{k}_{i=0}A[i]B[k-i]\)

\(A[i],B[k-1]\)分別為\(A\)中\(x^i\)的系數和\(B\)中\(x^{k-i}\)的系數，那么我們知道：\(x^i*x^{k-i}=x^k\)，故所有的\(x^i*x^{k-i}\)的系數對\(x\)的次項系數均有貢獻，需要累加。

那么就可以寫成：\(C[k]=\sum\limits_{i+j=k}A[i]B[j]\)

那么形如\(C[k]=\sum\limits_{i⊕j=k}A[i]B[j]\)(\(⊕\)為某種運算)的式子稱為卷積。

那么，多項式的乘法就是對於加法的卷積。

我們仍然可以對多項式相乘的模板打個暴力：

    #include<bits/stdc++.h>
    using namespace std;

    const int SIZE = 1e6 + 10;

    long long A[SIZE], B[SIZE], C[SIZE];
    int n, m;

    int main() {
        scanf("%d %d", &n, &m);
        n++, m++;
        for (int i = 0; i < n; i++) scanf("%lld", A + i);
        for (int i = 0; i < m; i++) scanf("%lld", B + i);
        for (int k = 0; k < n+m-1; k++)
            for (int i = 0; i <= k; i++)
                C[k] += A[i] * B[k - i];
        for (int i = 0; i< n + m - 1; i++) 
            printf("%lld ", C[i]);
        return 0;
    }

模板題親測 44 分。

---

算法理論

上面的主要部分還是寫了一個\(O(N^2)\)的暴力，那我們要尋找一個更快的辦法：

就是它！\(FFT\)算法

前面我們已經說過：\(FFT\)算法是一個\(O(NlogN)\)的算法，但常數較大(相信你們會卡常)，並且，\(FFT\)是\(DFT\)和\(IDFT\)的快速算法：

---

思想：

首先你要會點值表達式，當然我們已經講過。

然后我們來看\(A(x)=x^2+3x-2\)（紅色）與\(B(x)=-x^2+6x+4\)(綠色)的取點得到點的點值表達式:

在\(A(x)\)上取三個點（紫色）\((-4,2)(-2,-4)(1,2)\)；在\(B(x)\)上取三個點\((0,4)(2,12)(6,4)\)，然后轉為系數表達\(O(N^2)\)相乘然后再帶入\(x\)得到\(C\)的點值表達式？這時絕對不行的，我們上面說過，點值表達式相乘可以在\(O(n)\)內解決，那怎么辦呢？

我們可以取\(3\)個\(x\)相同的點：

這樣就行了？那是不可能的。我們知道，兩個二次項式相乘得到的是一個四次項式，而\(C(x)=y_1*y_2\)，只能得到三個點，這肯定是不行的。

這很明顯我們需要\(2n+1\)個點，不夠怎么辦呢？再添兩個不就好了？
下面添加了\(I,J\)(橫坐標相同)，\(G,H\)(橫坐標相同)，這樣就行了。

那么，我們只需要進行\(2n+1\)次乘法即可，省略常數（都說了常數巨大了），復雜度\(O(N)\)

神仙問題：\(tql\) \(but\) 這有用嗎

當然有用啊：

算法流程：

\[\text{系數表達式}\implies \text{點值表達式} \implies \text{系數表達式} \]

這不就有個用了嗎？

但這個思路僅僅是\(DFT+IDFT\)，而\(FFT\)是用來加速\(DFT\)和\(IDFT\)的。

如果讓我們求點值表達式，那我們肯定會帶整數進去，因為這樣好計算。但是傅里葉，（不知道是不是腦抽了），把復數帶進了多項式中，然后神奇的事情就發生了……

我們把多項式如下拆開：

\(F(x)=(a_0+a_2x^2+…+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3+…+a_{n-1}x^{n-1})\)

（保證\(n=2^k,k∈N^*\))

然后設兩個有\(n/2\)次項的多項式\(FL(x)\)和\(FR(x)\)

\[FL(x)=a_0+a_2x+…+a_{n-2}x^{\frac{n}{2}-1} \]

\[FR(x)=a_1+a_3x+…+a_{n-1}x^{\frac{n}{2}-1} \]

\[F(x)=FL(x^2)+xFR(x^2) \]

設\(k<\frac{n}{2}\)，把\(ω_n^k\)代入\(F(x)\)：\(F(ω_n^k)=FL(ω_n^{2k})+ω_n^kFR(ω_n^{2k})\)

接着，就有\(F(ω_n^k)=FL(ω_{n/2}^{k})+ω_n^kFR(ω_{n/2}^{k})\)

那么，如果我們知道\(FL(x),FR(x)\)分別在\(ω^0_{n/2},ω^1_{n/2},…,ω^{n/2-1}_{n/2}\)處的點值表示，那么我們就可以\(O(N)\)內求出\(F(x)\)在\(ω^0_{n},ω^1_{n},…,ω^{n/2-1}_{n}\)處的點值表示了。

接下來我們的目標就是求\(F(x)\)在\(ω^{n/2}_{n},ω^{n/2+1}_{n},…,ω^{n-1}_{n}\)處的點值表示。

然后……

繼續設\(k<n/2\)，然后把\(ω^{k+n/2}_{n}\)代入\(F(x)\)中。

運用一系列單位根的性質，讀者可以自行證明出：

\[F(ω^{k+n/2}_{n})=FL(ω^{k}_{n/2})-ω^{k}_{n}FR(ω^{k}_{n/2}) \]

然后我們就可以求出\(F(x)\)在\(ω^0_{n},ω^1_{n},…,ω^{n-1}_{n}\)處的點值表示，接着，我們就實現了\(DFT\)（聽着好累的樣子……）

但是我們還不知道\(FL(x),FR(x)\)在\(ω^0_{n},ω^1_{n},…,ω^{n/2-1}_{n}\)處的點值表示。難道要暴力代入嗎？這樣反而會使時間復雜度上升。

但是我們想想：\(FL(x)\)和\(FR(x)\)是由原來問題\(F(x)\)轉化來的，怎么轉化的？看到這里，你應該想到了分治。我們知道，分治的時間復雜度是\(logN\)，好像完成了一半了？

然后，因為\(FL(x),FR(x)\)是\(F(x)\)的子問題，那么我們可以對它們繼續進行分解成一個個小的子問題，然后再合並，類似於歸並排序。

\(FFT\)分治的合並過程根據：

1. \(F(ω^{k}_{n})=FL(ω^{k}_{n/2})+ω^{k}_{n}FR(ω^{k}_{n/2})\)

2. \(F(ω^{k+n/2}_{n})=FL(ω^{k}_{n/2})-ω^{k}_{n}FR(ω^{k}_{n/2})\)

從理論到實現

上文的所有證明都基於\(2=2^k,k∈N^*\)，但是我們做題的時候並沒有這一點，所以我們可以補一些系數為\(0\)的項，這對計算結果沒有影響。

---

1. 預處理單位根

我們知道一個性質\(ω_n^k=(ω_n^1)^k\)

那么，我們只要知道\(ω_n^1\)就可以求出\(ω_n^0,ω_n^1,ω_n^2,…,ω_n^{n-1}\)。

怎么求\(ω_n^1\)？

解三角形不就好了？已知\(|OA|=1\),幅角是\(\frac{1}{n}\)圓周，那么：

(由於\(c++\)下三角函數用弧度制表示，以下全部使用弧度制。)

設幅角\(α=\frac{2π}{n}\),，則\(ω_n^1(cos\big(\frac{2π}{n} \big),sin\big(\frac{2π}{n} \big))\)。

這邊要牢記，\(x\)是實軸，\(y\)是虛軸：

    struct complex { //復數 
        long double a, b;
        complex() { a = b = 0;}
            //重載復數四則運算 
        complex operator + (const complex x) const {
            complex temp;
            temp.a = a + x.a;
            temp.b = b + x.b;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.a = a - x.a;
            temp.b = b - x.b;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.a = a * x.a - b * x.b;
            temp.b = a * x.b + b * x.a;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.a * x.a + x.b * x.b;
            complex temp;
            temp.a = (a * x.a + b * x.b) / div;
            temp.b = (b * x.a - a * x.b) / div;
            return temp; 
        }
    }w[SIZE],temp, buf;

    void unit_roots(int n) {
        const long double pi = 3.14159265358979;
        //complex ;
        temp.a = cos(2*pi/n), temp.b = sin(2*pi/n);
        w[0].a = 1, w[0].b = 0;
        buf.a = 1, buf.b = 0;
        for(int i = 1; i < n ; i++) {
            buf *= temp;
            w[i] = buf; 
        }
        for(int i = 0; i < n; i++)
            printf("w_%d %.3Lf %.3Lf\n", i, w[i].a, w[i].b);	
    } 

然后我們就可以實現\(DFT\)了

因為還沒講\(IDFT\)，先寫\(DFT\):

遞歸\(DFT\)實現（復數板子上面有了下面就不寫了）

    int n, m;
    Complex f[SIZE];
    void dft(Complex *f, int len) { //當前子問題長度
        if(len == 0) return;
        Complex fl[len + 10], fr[len + 10];
        for(int i = 0 ; i < len; i++)
            fl[i] = f[i << 1], fr[i] = f[i << 1 | 1];
        dft(fl, len >> 1);
        dft(fr, len >> 1);
        len *= 2;
        Complex temp, buf;
        temp.a = cos(2*pi/n), temp.b = sin(2*pi/n);
        buf.a = 1;
        for(int i = 0 ; i < len / 2; i++) {
            f[i] = fl[i] + buf * fr[i];
            f[i + len / 2] = fl[i] - buf * fr[i];
            buf *= temp;
        }
    }
    int main() {

        scanf("%d", &n);
        for(int i = 0; i < n; i++) scanf("%Lf", &f[i].a);
        for(m = 1; m < n; m <<= 1);//補到2的冪次方
        dft(f, m >> 1);
        for(int i = 0; i < m; i++) printf("%.3Lf ", f[i].a);
        return 0;
    }

事實上虛數可以用(小心TLE)：

            complex <double> f;

然后我們得到了一堆點值表達……

這並沒有什么用

於是，我們現在來學習\(IDFT\)——\(DFT\)的逆運算。

剛剛的多項式為\(F(x)\)，系數是\(A\)數列，點值表示成\(M\)，\(M[k]=\sum \limits_{i=0}^{n-1}(ω_n^{k})^i*F[i]\)

然后就有\(F[k]=\frac{1}{n}\sum \limits_{i=0}^{n-1}(ω_n^{-k})^i*M[i]\)

結論記住就好了，可以自己嘗試證明。

我們可以先求出和式，然后再除以\(n\)。

接下來我們求出\(ω_n^0,ω_n^{-1},…ω_n^{-n+1}\)

並且仍有：\(ω_n^{-j}=(ω_n^{-1})^j\)

所以我們還是求出\(ω_0^{-1}\)即可，並且，\(ω_0^{-1}\)與\(ω_0^{1}\)關於\(x\)軸對稱。

    #include<iostream>
    #include<cstdio>
    #include<cmath>
    using namespace std;
    const long double pi = 3.141592653589793;
    const int SIZE = 4e6 + 10;
    struct complex { //復數 
        long double real, imag;
        complex() { real = imag = (double)0;}
            //重載復數四則運算 
        complex operator + (const complex x) const {
            complex temp;
            temp.real = real + x.real;
            temp.imag = imag + x.imag;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.real = real - x.real;
            temp.imag = imag - x.imag;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.real = real * x.real - imag * x.imag;
            temp.imag = real * x.imag + imag * x.real;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.real * x.real + x.imag * x.imag;
            complex temp;
            temp.real = (real * x.real + imag * x.imag) / div;
            temp.imag = (imag * x.real - real * x.imag) / div;
            return temp; 
        } 
        void conj() { imag = -imag; } 
    }a[SIZE], b[SIZE];

    complex omega(int n, int k) {
        complex temp;
        temp.real = cos(2 * pi * k / n);
        temp.imag = sin(2 * pi * k / n);
        return temp;
    }	
    int n, m, p;		
    void fft(complex *f, int len, int inv) {
        if(len == 1) return;
        int mid = len >>1;
        complex fl[mid + 10], fr[mid + 10];
        for(int i = 0; i < mid; i++)
            fl[i] = f[i << 1], fr[i] = f[i << 1 | 1];
        fft(fl, mid, inv);
        fft(fr, mid, inv);
        complex temp, buf; 
        temp.real = cos(2 * pi / len);
         temp.imag = sin(2 * pi / len);
        buf.real = 1, buf.imag = 0;
        if(inv) temp.conj();
        for(int i = 0; i < mid; i++, buf *= temp) {
            complex s = buf * fr[i];
            f[i] = fl[i] + s;
            f[i + mid] = fl[i] - s;
        }
    }

    int main() {

        scanf("%d %d", &n, &m);
        for(p = 1; p < n + m + 1; p <<= 1);
        for(int i = 0; i <= n; i++) 
            scanf("%Lf", &a[i].real);
        for(int i = 0; i <= m; i++) 
            scanf("%Lf", &b[i].real);
        fft(a, p, 0);
        fft(b, p, 0);
        for(int i = 0; i < p; i++) 
            a[i] *= b[i];
        fft(a, p, 1);
        for(int i = 0; i <= n + m; i++)
            printf("%d ", (int)(a[i].real / p + 0.5));
	return 0;
}

接着TLE了……(絲毫沒有卡常的欲望)

因為遞歸常數巨大，並且很容易爆空間。(但實際上是因為寫了\(long \ double\)炸了)

下面是AC的遞歸代碼，以后不要寫long double 了

    #include<iostream>
    #include<cstdio>
    #include<cmath>
    using namespace std;
    const double pi = 3.141592653589793;
    const int SIZE = 4e6 + 10;
    struct complex { //復數 
        double real, imag;
        complex() { real = imag = (double)0;}
            //重載復數四則運算 
        complex operator + (const complex x) const {
            complex temp;
            temp.real = real + x.real;
            temp.imag = imag + x.imag;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.real = real - x.real;
            temp.imag = imag - x.imag;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.real = real * x.real - imag * x.imag;
            temp.imag = real * x.imag + imag * x.real;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.real * x.real + x.imag * x.imag;
            complex temp;
            temp.real = (real * x.real + imag * x.imag) / div;
            temp.imag = (imag * x.real - real * x.imag) / div;
            return temp; 
        } 

        void conj() { imag = -imag;}
    }a[SIZE], b[SIZE];

    complex omega(int n, int k) {
        complex temp;
        temp.real = cos(2 * pi * k / n);
        temp.imag = sin(2 * pi * k / n);
        return temp;
    }	
    int n, m, p;		
    void fft(complex *f, int len, int inv) {
        if(len == 1) return;
        int mid = len >>1;
        complex fl[mid + 10], fr[mid + 10];
        for(int i = 0; i < mid; i++)
            fl[i] = f[i << 1], fr[i] = f[i << 1 | 1];
        fft(fl, mid, inv);
        fft(fr, mid, inv);
        complex temp, buf; 
        temp.real = cos(2 * pi / len);
         temp.imag = sin(2 * pi / len);
        buf.real = 1, buf.imag = 0;
        if(inv) temp.conj();
        for(int i = 0; i < mid; i++, buf *= temp) {
            complex s = buf * fr[i];
            f[i] = fl[i] + s;
            f[i + mid] = fl[i] - s;
        }
    }

    int main() {

        scanf("%d %d", &n, &m);
        for(p = 1; p < n + m + 1; p <<= 1);
        for(int i = 0; i <= n; i++) 
            scanf("%lf", &a[i].real);
        for(int i = 0; i <= m; i++) 
            scanf("%lf", &b[i].real);
        fft(a, p, 0);
        fft(b, p, 0);
        for(int i = 0; i < p; i++) 
            a[i] *= b[i];
        fft(a, p, 1);
        for(int i = 0; i <= n + m; i++)
            printf("%d ", (int)(a[i].real / p + 0.5));
        return 0;
    }

實際上跑的很快的遞歸評測記錄。

那有沒有更優化的辦法呢？這是肯定的。

---

蝴蝶變換

先看一張圖

這張圖是什么？這不是我們分治的順序嗎？

我在一開始和結束都打上了二進制，你有沒有什么發現？

是的，每一個數的二進制都反過來了。

然后我們的任務就轉換成求二進制反序：

    rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << l)//l是二進制位數
    //rev[i]是i翻轉得到的

結論牢記即可。

還記得\(DP\)嗎？（貌似沒什么關系）

我們用循環實現\(FFT\)，也就像\(DP\)一樣定義狀態、階段決策：

• 第一層循環枚舉變換區間大小（階段），第二層循環枚舉開頭（決策），第三層處理區間。

    #include<iostream>
    #include<cstdio>
    #include<cmath>
    using namespace std;
    const double pi = 3.141592653589793;
    const int SIZE = 4e6 + 10;
    struct complex { //復數 
        double real, imag;
        complex() { real = imag = (double)0;}
            //重載復數四則運算 
        complex operator + (const complex x) const {
            complex temp;
            temp.real = real + x.real;
            temp.imag = imag + x.imag;
            return temp; 
        }

        complex operator += (const complex &x) {
            *this = *this + x;
            return *this;
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.real = real - x.real;
            temp.imag = imag - x.imag;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.real = real * x.real - imag * x.imag;
            temp.imag = real * x.imag + imag * x.real;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.real * x.real + x.imag * x.imag;
            complex temp;
            temp.real = (real * x.real + imag * x.imag) / div;
            temp.imag = (imag * x.real - real * x.imag) / div;
            return temp; 
        } 
    }a[SIZE], b[SIZE];
    int n, m, p, L = -1;
    int rev[SIZE];

    void fft(complex *f, int inv) {
        for(int i = 0 ; i < p ; i++) 
            if(i < rev[i]) 
                swap(f[i], f[rev[i]]);
        for(int len = 2; len <= p; len <<= 1) {
            int length = len >> 1;
            complex temp;
            temp.real = cos(pi / length);
            temp.imag = sin(pi / length);
            if(inv) temp.imag = -temp.imag;
            for(int l = 0; l < p; l += len) {
                complex buf; buf.real = 1;
                for(int k = l; k < l + length; k++) {
                    complex t = buf * f[k + length];
                    f[k + length] = f[k] - t;
                    f[k] += t;
                    buf *= temp;
                }
            }
        }
    }

    int main() {
        scanf("%d %d", &n, &m);
        for(int i = 0; i <= n; i++)
            scanf("%lf", &a[i].real);
        for(int i = 0; i <= m; i++)
            scanf("%lf", &b[i].real);
        for(p = 1; p < n + m + 1; p <<= 1, L++);
        //l是位數，因為會多算一位，一開始初值直接給-1
        for(int i = 1; i < p; i++)
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L);
        fft(a, 0), fft(b, 0);
        for(int i = 0; i < p; i++)
            a[i] *= b[i];
        fft(a, 1);
        for(int i = 0; i <= n + m; i++)
            printf("%d ", (int)(a[i].real / p + 0.5));
        return 0;
    } 

評測記錄

NTT——快速數論變換

• 首先，學習\(NTT\)算法你要學會\(FFT\)（理解也可）（剛剛我們才講完你不會忘了吧……）（我認為我寫的很詳細了呢\(QAQ\)）。

• 接下來，你要確保你會原根這個東西（起碼要知道），不過我在后面還會再提，具體可以先看我的這一篇博客：同余相關，在比較下面的部分（由於還待完善，所以只是寫了一點）。

• 然后就是希望大家看了這一篇博客對\(NTT\)有所了解，可能我講的不好，但是起碼你要知道\(NTT\)是什么。

---

前置知識

（怎么又是前置知識）

\(FFT\)

上面已經說了先理解\(FFT\)，並且我也給了寫的比較詳細的博客，這里不再展開。

但其實\(FFT\)並不是重點，而是里面用於優化的很重要的一個東西：\(\color{red}{\text{單位根}}\)。

---

• 原根

但是如果 \(FFT\) 要求取模呢？

\[\color{red}{What?\text{你不會？}} \]

好像因為 \(FFT\) 用的單位根不太好取模（復數），所以我們可能需要找一個東西來替代下單位根，那么就要求這個東西和單位根具有類似的性質。然后科學家們就開始尋找了……（也不知道是誰找到的）接下來他們找到了原根。

---

首先我們要了解一下啥是原根。（如果你看了上面給出的博客你應該知道了）

原根的定義

其實在那一篇博客里的原根定義很容易讓人摸不着頭腦……為啥呢？

\[\color{green}{δ_m(a)=φ(m) \text{？} what \ is \ this?} \]

因為它很不好求……那我們不妨換一種說法：

設 \(p\) 為質數，則根據 \(\color{orange}{\text{費馬小定理}}\)，有：

\[\text{對於任意 }p\nmid a,a^{p-1}\equiv1 \pmod{p} \]

\[\color{red}{\text{那么稱 $g$ 為模 $p$ 的原根，當且僅當}\{g^0,g^1,…,g^{p-2}\}\text{對} \ p \ \text{的模數互不相同。}} \]

這樣可能會更好理解些？至於為什么是上面那個綠綠的，我也不知道啊……

• 補充：本原單位根

其實定義很簡單，對於一些比較特殊的\(n\)次單位根，它的\(0,1,…,n-1\)次冪恰好生成了所有的\(n\)次單位根，我們把這樣的單位根稱為本原單位根。

在復數域 \(C\) 上，\(\omega_n=\exp\frac{2\pi i}{n}=\cos\frac{2\pi}{n}+i\sin\frac{2\pi}{n}\) 是一個本原單位根。

（以上不是重點，下面這句話才是重點）

在有限域（即模質數）上，本原單位根和數論中的原根有關！

這也是我補充本原單位根的原因了……因為后面我們還會提到它。

並且我們可以證明的是，模質數的原根總是存在的（不知道為啥你就問度娘吧）。

並且原根的性質和單位根非常相似哦~

\(NTT\) 里原根的性質

換句話說，在模 \(p\) 意義下，\(g\) 可以被看做一個 \(p-1\) 次本原單位根。

反正感性理解就好，結論不就是用來記的嘛……

設 \(n\mid p-1\) ,我們不妨令 \(\omega_n=g^{\frac{p-1}{n}}\)（這個其實很像單位根的求法），然后你就可以得到一大串東西了：（好像也不多）

\[\omega_n^n=(g^{\frac{p-1}{n}})^n=g^{p-1}\equiv1\pmod p \]

並且因為單位根的定義，你會知道 \(\omega_n^0,\omega_n^0,…,\omega_n^{n-1}\) 互不相同……這不就完了嘛 $QAQ $

因此此時的 \(\omega_n\) 在模 \(p\) 意義下具有 \(n\) 次本原單位根的性質，所以我們利用這個東西來定義 \(DFT\) 即可，這被稱為 \(\color{red}{NTT}\)。

原根的求法

求模質數原根的辦法： 將 \(p-1\) 分解質因數，即 \(p-1=p_1^{a_1}p_2^{a2}…p_{n}^{a_n}\) 為標准的\(p-1\)分解式，那么若 \(g^{\frac{p-1}{p_i}}\ne 1\pmod p\) ，則\(g\)是模 \(p\) 意義下的一個原根。

求模合數原根的辦法： 其實你只需要把上面的 \(p-1\) 換成 \(φ(p)\) 即可。

證明：（假設法）
  先對於第一種模質數的情況：
  假設存在一個 \(t<φ(p)=p-1\) 使得 \(g^t\equiv1\pmod p\) 且 \(\forall i∈[1,n],g^{\frac{p-1}{p_i}}\ne 1 \pmod p\)。由裴蜀定理可知：一定存在一組 \((a,b)\) 使得：\(a*t=(p-1)*b+gcd(t,p-1)\) 成立。由歐拉定理/費馬小定理：\(g^{p-1}\equiv 1 \pmod p\)。所以 \(g^{a*t}\equiv g^{(p-1)*b+gcd(t,p-1)}\equiv1\pmod p\)。因為 \(t<p-1\)，所以又有 \(gcd(t,p-1)<p-1\)。又因為 \(gcd(t,p-1)\mid p-1\)，於是 \(\exists i ∈[1,n],gcd(t,p-1)\mid{g^{\frac{p-1}{p_i}}}\)。那么就可以得到 \(g^{\frac{p-1}{p_i}}\equiv g^{gcd(t, p-1)}\equiv 1 \pmod p\)。
 因此假設不成立。所以原命題成立。
\(Q.E.D.\)

那么上述結論自然可以推廣到模合數而把 \(p-1\) 改為 \(φ(p)\) 了。

以上的證明是針對於這一句話的：\(δ_m(a)=φ(m)\)，也就是上述求法成立的條件。我們從另一個性質出發，即 \(g^0,g^1,…,g^{p-1}\) 模\(p\)的余數各不相同，那么可以得知這些余數的集合為：\(\{1,2,…p-1\}\)（不分順序）。我們知道取模是滿足乘法性質的。也就是說，如果\(\exists i∈[1,n]\)使得\(g^{\frac{p-1}{p^i}}=1\)，我們知道唯有 \(g^0 \bmod p=1 \bmod p = 1\) ，那么很明顯 \(\frac{p-1}{p^i}\) 是不會為 \(0\) 的，那么一旦出現第一個數字的重復，也就是 \(1\)，那么很明顯不滿足原根的上述性質。

那么我就很好奇，為什么重復的\(1\)一定會出現在某個\(g^{\frac{p-1}{p^i}}\)當中呢？上面其實也已經證明，如果存在某個 \(t<φ(p)\) 使得 \(g^t\equiv1\pmod p\) 成立，那么必定會存在一個 \(g^\frac{p-1}{p^i}\bmod p=1\) 。並且我認為這樣可能更好理解。

反正證明這種東西感性理解一下就好了……

然后原根你也會求了，所以用原根 \(g\) 去代替 \(FFT\) 中的 \(\omega_n\) (即單位根),就可以完成 \(FFT\) 了。

NTT 模數

首先根據 \(FFT\) 我們知道 \(g^{\frac{p-1}{n}}\)中的 \(n\) 是 \(2\) 的冪次，所以在做 \(FFT\) 的時候我們最好構造形如 \(p=a*2^k+1\)的質數，這樣就可以滿足上面的需求了。這樣的質數最好取大一點（精度），所以這樣的質數可以取：

\(p=a*2^k+1\)	$\ $ \(a\)	\(k\)	\(g\)
\(3\)	\(1\)	\(1\)	\(2\)
\(5\)	\(1\)	\(2\)	\(2\)
\(17\)	\(1\)	\(4\)	\(3\)
\(97\)	\(3\)	\(5\)	\(5\)
\(193\)	\(3\)	\(6\)	\(5\)
\(257\)	\(1\)	\(8\)	\(3\)
\(7681\)	\(15\)	\(9\)	\(17\)
\(12289\)	\(3\)	\(12\)	\(11\)
\(40961\)	\(5\)	\(13\)	\(3\)
\(65537\)	\(1\)	\(16\)	\(3\)
\(786433\)	\(3\)	\(18\)	\(10\)
\(5767169\)	\(11\)	\(19\)	\(3\)
\(7340033\)	\(7\)	\(20\)	\(3\)
\(23068673\)	\(11\)	\(21\)	\(3\)
\(104857601\)	\(25\)	\(22\)	\(3\)
\(167772161\)	\(5\)	\(25\)	\(3\)
\(469762049\)	\(7\)	\(26\)	\(3\)
\(\color{red}{998244353}\)	\(119\)	\(23\)	\(3\)
\(\color{red}{1004535809}\)	\(479\)	\(21\)	\(3\)
\(1998585857\)	\(953\)	\(21\)	\(3\)
\(2013265921\)	\(15\)	\(27\)	\(31\)
\(2281701377\)	\(17\)	\(27\)	\(3\)
\(3221225473\)	\(3\)	\(30\)	\(5\)
\(75161927681\)	\(35\)	\(31\)	\(3\)
\(77309411329\)	\(9\)	\(33\)	\(7\)
\(206158430209\)	\(3\)	\(36\)	\(22\)
\(2061584302081\)	\(15\)	\(37\)	\(7\)
\(2748779069441\)	\(5\)	\(39\)	\(3\)
\(6597069766657\)	\(3\)	\(41\)	\(5\)
\(39582418599937\)	\(9\)	\(42\)	\(5\)
\(79164837199873\)	\(9\)	\(43\)	\(5\)
\(263882790666241\)	\(15\)	\(44\)	\(7\)
\(1231453023109121\)	\(35\)	\(45\)	\(3\)
\(1337006139375617\)	\(19\)	\(46\)	\(3\)
\(3799912185593857\)	\(27\)	\(47\)	\(5\)
\(4222124650659841\)	\(15\)	\(48\)	\(10\)
\(7881299347898369\)	\(7\)	\(50\)	\(6\)
\(31525197391593473\)	\(7\)	\(52\)	\(3\)
\(180143985094819841\)	\(5\)	\(55\)	\(6\)
\(1945555039024054273\)	\(27\)	\(56\)	\(5\)
\(4179340454199820289\)	\(29\)	\(57\)	\(3\)

理論到實現

第一件事，我們先來回憶一下 \(FFT\)，而 \(FFT\) 分為兩個步驟：\(DFT\) 與 \(IDFT\)。

然后我們知道單位根的定義是對於復數域 \(C\)，有 \(z^n=1\) 的復數就是一個 \(n\)次單位根，\(n\) 次單位根有 \(n\) 個，為：\(e^{\frac{2\pi k}{n}i}, \ k = \{0,1,2…,n-1\}\)（\(i\)是虛數單位）。（你用三角表示也是可以的，但是不好理解）

然后這里就用 \(g^{\frac{p-1}{n}}\) 代替 上面那個就可以，因為你會發現他們具有相同的性質。

唯一具有卷積性質的變換就是 \(DFT\)，按照上面的式子，\(DFT\) 的變換式就是：

\[X_k=\sum_{n=0}^{N-1}x_ne^{-\frac{2\pi}{N}nki} \ \ \ \ \ \ k=0,1,2…,N-1 \]

\(DFT\) 反演（也就是逆變換）公式有：

\[x_n=\frac{1}{N}\sum_{k=0}^{N-1}X_ke^{\frac{2\pi}{N}nki} \ \ \ \ \ \ n=0,1,2…N-1 \]

那么由於 \(g^{\frac{p-1}{n}}\) 和 \(e^{\frac{2\pi k}{n}i}, \ k = \{0,1,2…,n-1\}\) 等價，所以：

\[e^{-\frac{2\pi}{N}nki}\equiv g^{\frac{p-1}{N}}\pmod p \]

接着我們就得到了快速數論變換的公式：

\[X_k=\sum_{n=0}^{N-1}x_ng^{nk}\pmod p \ \ \ \ \ \ k=0,1,2…,N-1 \]

反演：

\[x_n=\frac{1}{N}\sum_{k=0}^{N-1}X_kg^{-nk}\pmod p\ \ \ \ \ \ n=0,1,2…N-1 \]

這樣我們就成功把復數域轉到了整數域，然后剩下的東西在\(\mod p\)意義下考慮即可。上面也已經講過素數的構造方法。

• 補充：原根為何可以代替單位根？

考慮單位根的性質（也就是我們為什么使用單位根）：

1. \(n\) 個單位根互不相同。那么\(g^0,g^1,…,g^{p-2}\) 在模 \(p\) 意義下也不相同，成立。

2. \(z^n=1\)（\(z∈C\)）。根據費馬小定理：\(g^{p-1}\equiv 1\pmod p\)，成立。

3. 單位根對稱分布。其實這對於原根也成立。

證明：
  \(∵(g^{\frac{p-1}{2}})^2\equiv g^{p-1}\equiv 1 \pmod p\)
  又\(∵g^i \ne g^j(i\ne j \text{且}0≤i,j≤p-1)\)
  \(∴g^{\frac{p-1}{2}}= \ne 1,\text{則}g^{\frac{p-1}{2}}=-1\)
\(Q.E.D.\)

好了，我們真的可以使用 \(NTT\) 來完成 \(FFT\) 所能完成和所不能完成的事了。

可以嘗試過一下模板題：\(P3803\)【模板】多項式乘法(\(FFT\))

雖然題目沒有要求模數，但是仍然可以取一個比較大的模數從而完成。

這里取 \(p=998244353,g=3\)。

我們結合代碼來理解 \(NTT\) :

#include<bits/stdc++.h>
using namespace std;
const int N = (1 << 21) + 100, P = 998244353, G = 3;

inline int read() {
	static int x, f;
	static char c;
	while(!isdigit(c = getchar()) && c != '-');
	x = (f = c == '-') ? 0 : c - '0';
	while(isdigit(c = getchar())) 
		x = (x << 3) + (x << 1) + (c ^ 48);
	return f ? -x + 1 : x;
}

typedef long long ll;
inline ll power(ll a, ll k) {
	ll base = 1;
	for(; k; k >>= 1) {
		if(k & 1) base = (base * a) % P;
		a = (a * a) % P;
	}	
	return base;
}

int rev[1 << 21], limit = 1;
ll a[N], b[N];
inline void NTT(ll *f, int inv) {
	for(int i = 0; i < limit; i++)
		if(i < rev[i]) swap(f[i], f[rev[i]]);
	
	for(int len = 1; len < limit; len <<= 1) {
		ll temp = power(G, (P - 1) / (len << 1));
		if(inv) temp = power(temp, P - 2);
		//反NTT時，除法全部變為乘逆元 
		for(int l = 0; l < limit; l += len * 2) {
			ll w = 1;
			for(int k = 0; k < len; k++) {
				int x = f[l + k], y = w * f[l + k + len] % P;
				f[l + k] = (x + y) % P;
				f[l + k + len] = (x - y + P) % P;
				w = (w * temp) % P;
			}
		}
	} 
}

int main() {
	int n = read(), m = read(), L = -1;
	for(int i = 0; i <= n; i++) a[i] = (read() + P) % P;
	for(int i = 0; i <= m; i++) b[i] = (read() + P) % P;
	for( ; limit <= n + m; limit <<= 1, L++); 
	for(int i = 0; i < limit; i++)
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L);
	NTT(a, 0), NTT(b, 0);
	for(int i = 0; i < limit; i++)
		a[i] = (a[i] * b[i]) % P;
	NTT(a, 1);
	ll inv = power(limit, P - 2);
	for(int i = 0; i <= n + m; i++)
		printf("%lld ", (a[i] * inv) % P);
	
	return 0;
}

然后真的比\(FFT\)快好多呀……模板題評測記錄：

\(\text{遞歸}FFT\)：點我

\(\text{蝴蝶變換}FFT\)：點我

\(NTT\)：點我

---

Upd20200918:文章當中還有一些錯誤，近期本人正在矯正。