[多項式全家桶]快速傅立葉變換 - FFT

阿新 • • 發佈：2021-06-15

#1.0 多項式

#1.1 環和域

這一節沒啥卵用...

#1.1.1 環

如果一個非空集合 \(R\) 上定義了兩個二元運算 \(+,\cdot\)（分別稱為加法和乘法），滿足：

\((R,+)\) 構成 \(\text{Abel}\) 群（滿足交換律的群）；
乘法結合律：\((a\cdot b)\cdot c=a\cdot(b\cdot c),\forall\ a,b,c\in R;\)
分配率：\((a+b)\cdot c=a\cdot c+b\cdot c,\ c\cdot(a+b)=c\cdot a+c\cdot b,\forall\ a,b,c\in R,\)

則稱 \(R\)

關於運算 \(+,\cdot\) 構成一個環，記為 \((R;+,\cdot)\)，或簡記為 \(R.\)

若環 \(R\) 中成立

乘法交換律：\(a\cdot b=b\cdot a,\forall\ a,b\in R,\)

則稱 \(R\) 為交換環。

環 \(R\) 中若存在乘法么元，即存在 \(e\in R\)，使得對任意的 \(a\in R\)，恆有

\[e\cdot a=a\cdot e=a, \]

則稱 \(R\) 為么環。

設 \(R\) 為么環，\(a\in R\)，如果存在 \(b\in R\) 使得 \(b\cdot a=e\)，則稱 \(b\) 為 \(a\) 的一個左逆元

。類似的，如果存在 \(c\in R\) 使得 \(a\cdot c=e\)，則稱 \(c\) 是 \(a\) 的一個右逆元。如果 \(b\in R\) 使得 \(ba=ab=e\)，則稱 \(b\) 為 \(a\) 的逆元，記為 \(a^{-1}\)。同時稱 \(a\) 為可逆元或單位元。

#1.1.2 域

設 \(D\) 是含有至少兩個元素的么環。如果 \(D\) 的每個非零元素都可逆，則稱 \(D\) 是一個體。具有乘法交換律的體稱為域。

設 \(p\) 是一個質數，那麼模 \(p\) 意義下的集合 \({0,1,\cdots,p−1}\) 是一個域，記為 \(F_p\)。

下文的多項式都在域上。

#1.2 多項式

#1.2.1 定義

設 \(R\) 是一個環，\(a_0,a_1,\cdots,a_n\) 為 \(R\) 中的元素且 \(a_n\ne0,x_0\) 始終表示 \(1\)，則

\[f(x)=\sum\limits_{k=0}^na_kx^k \]

稱為 \(R\) 上次數為 \(n\) 的多項式。

#1.2.2 卷積

對於陣列 \(a,b\)，令

\[c_k=\sum\limits_{i+j=k}a_ib_j=\sum\limits_{i=0}^ka_ib_{k-i}=\sum\limits_{j=0}^ka_{k-j}b_j, \]

則稱 \(c\) 為 \(a\) 和 \(b\) 的卷積。

#1.2.3 多項式乘法

設 \(f(x)=\sum\limits_{k=0}^na_kx^k\) 和 \(g(x)=\sum\limits_{k=0}^mb_kx^k\)，利用分配率可得

\[f(x)g(x)=\sum\limits_{k=0}^{n+m}c_kx^k. \]

其中 \(c\) 是 \(a\) 和 \(b\) 的卷積。

#1.2.4 多項式的點值表示

設 \(f(x)=\sum\limits_{k=0}^na_kx^k\) 為次數不大於 \(n\) 的多項式，對於 \(n+1\) 個點 \(x_0,x_1,\cdots,x_n\)，有

\[\begin{pmatrix}f(x_0)\\f(x_1)\\f(x_1)\\\vdots\\f(x_n)\end{pmatrix}=\begin{pmatrix}1&x_0&x_0^2&\cdots&x_0^n\\1&x_1&x_1^2&\cdots&x_1^n\\1&x_2&x_2^2&\cdots&x_2^n\\\vdots&\vdots&\vdots&\ddots&\vdots\\1&x_n&x_n^2&\cdots&x_n^n\end{pmatrix}\begin{pmatrix}a_0\\a_1\\a_2\\\vdots\\a_n\end{pmatrix} \]

若 \(x_0,x_1,\cdots,x_n\) 互不相同，則等式右邊的矩陣可逆，那麼

\[\begin{pmatrix}a_0\\a_1\\a_2\\\vdots\\a_n\end{pmatrix}=\begin{pmatrix}1&x_0&x_0^2&\cdots&x_0^n\\1&x_1&x_1^2&\cdots&x_1^n\\1&x_2&x_2^2&\cdots&x_2^n\\\vdots&\vdots&\vdots&\ddots&\vdots\\1&x_n&x_n^2&\cdots&x_n^n\end{pmatrix}^{-1}\begin{pmatrix}f(x_0)\\f(x_1)\\f(x_1)\\\vdots\\f(x_n)\end{pmatrix} \]

即 \(n+1\) 個點值對 \((x_k,f(x_k))\) 可以唯一確定 \(f(x).\)

那麼，我們來看如何用點值來處理多項式乘法。

設 \(f(x)=\sum\limits_{i=0}^na_ix^i,g(x)=\sum\limits_{i=0}^mb_ix^i\)，設 \(h(x)\) 為 \(f(x)\) 和 \(g(x)\) 的乘積，於是有

\[h(x)=f(x)g(x)=\sum\limits_{i=0}^{n+m}c_ix^i, \]

其中 \(c\) 是 \(a\) 和 \(b\) 的卷積。

那麼，我們只需要 \(n+m+1\) 個互不相同的點 \(x_0,x_1,x_2,\cdots,x_{n+m}\)，若已知所有的 \(f(x_k)\) 和 \(g(x_k)\)，只需計算所有的 \(h(x_k)\)，即 \(f(x_k)g(x_k)\) 的值便可唯一確定 \(h(x)\)，即 \(f(x)g(x).\)

當然，對於 \(f(x)+g(x)\) 也是類似的。

#1.2.5 多項式的根

設 \(F,k\) 是域且 \(F\subseteq K\)，\(f(x)\) 是 \(F\) 上的多項式， \(\alpha\) 是 \(K\) 中的元素。若 \(f(\alpha)=0\)，則稱 \(\alpha\) 為 \(f(x)\) 的一個根。

\(f(x)\) 的根不一定在 \(F\) 中。

#1.2.6 單位根

多項式 \(x^n-1\) 的根稱為 \(n\) 次單位根。

設 \(\omega_n=e^{\frac{2\pi i}{n}}\) 為複數，由尤拉公式

\[e^{ix}=\cos x+i\sin x \]

容易驗證 \(\omega_n^n=1\)；

對於質數 \(p=kn+1\)，設 \(\omega_n=g^{\frac{p-1}n}\) 為 \(F_p\) 中的元素，其中 \(g\) 為模 \(p\) 的原根，那麼顯然 \(\omega_n^n=1.\)

在上述兩種情況下，\(1,\omega_n,\omega_n^2,\cdots,\omega_n^{n-1}\) 為 \(n\) 個互不相同的 \(n\) 次單位根。

#1.2.7 單位根の性質

\(\omega_{2n}^{2k}=\omega_{n}^k\)

對應的點/向量是相同的。

\[(\cos\dfrac{4k\pi}{2n},\sin\dfrac{4k\pi}{2n})=(\cos\dfrac{2k\pi}{n},\sin\dfrac{2k\pi}{n}) \]

#2.0 快速傅立葉變換（FFT）

#2.1 離散傅立葉變換（DFT）

#2.1.1 定義

設 \(a\) 為大小為 \(n\) 的陣列，\(f(x)=\sum\limits_{k=0}^{n-1}a_kx^k\)，定義

\[\text{DFT}(a)=(f(1),f(\omega_n),f(\omega_n^2),\cdots,f(\omega_n^{n-1})), \]

稱為 \(a\) 的離散傅立葉變換，或 \(a\) 的 \(\text{DFT}.\)

\(\text{DFT}\) 將 \(f(x)\) 的係數陣列 \(a\) 變換成 \(f(x)\) 在每個 \(\omega_n^k\) 處的值，從而得到了 \(f(x)\) 的一個點值表示。

#2.1.2 IDFT

\(\text{DFT}\) 的逆變換稱為 \(\text{IDFT}\)，也記作 \(\text{DFT}^{−1}\)。\(\text{IDFT}\) 將 \(f(x)\) 在 \(n\) 個 \(n\) 次單位根處的點值表示變換為其係數表示。

只要將 \(\text{DFT}\) 中的 \(\omega_n\) 換成 \(\omega_n^{-1}\) 並將結果除以 \(n\)，就可以得到 \(\text{IDFT}\)。

更形式化地說，設 \(g(x)=\sum\limits_{k=0}^{n-1}f(\omega_n^k)x^k\)，則有

\[\begin{aligned} a=&\text{DFT}^{-1}(\text{DFT(a)})\\ =&\dfrac{1}{n}(g(1),g(\omega_n^{-1}),g(\omega_n^{-2}),\cdots,g(\omega_n^{-(n-1)})). \end{aligned} \]

#2.2 卷積

#2.2.1 迴圈卷積

對於大小為 \(n\) 的陣列 \(a,b\)，令

\[c_k=\sum\limits_{(i+j)\bmod\ n=k}a_ib_j=\sum\limits_{i+j=k}a_ib_j+\sum\limits_{i+j=n+k}a_ib_j \]

則稱 \(c\) 為 \(a\) 和 \(b\) 的迴圈卷積，用 \(a\otimes b\) 表示。

#2.2.2 卷積定理

設 \(f(x)=\sum\limits_{k=0}^{n-1}a_nx^k,g(x)=\sum\limits_{k=0}^{n-1}b_kx^k,h(x)=\sum\limits_{k=0}^{n-1}c_kx^k\)。對於任意正數 \(s\) 由於 \((\omega_n^s)^n=1\)，有

\[\begin{aligned} f(\omega_n^s)g(\omega_n^s)=&\sum\limits_{k=0}^{n-1}\left(\sum\limits_{i+j=k}a_ib_j+\sum\limits_{i+j=n+k}a_ib_j\right)(\omega_n^s)^k\\ =&h(\omega_n^s) \end{aligned} \]

注意這裡 \(c=a\otimes b.\)

由於 \(h(x)\) 是次數不大於 \(n-1\) 的多項式，由點值表示的唯一性得

\[\text{DFT}^{-1}(\text{DFT}(a)\circ\text{DFT}(b))=a\otimes b. \]

其中 \(\circ\) 為逐元素相乘的運算。

狹義的 \(\text{DFT}\) 僅指代複數域上的 \(\text{DFT}\)，因此常把 \(F_p\) 上的 \(\text{DFT}\) 稱為數論變換，或 \(\text{NTT}\)。當然，這兩種形式的 \(\text{DFT}\) 的原理是相同的。

#2.3 FFT 綜述

快速傅立葉變換，或 \(\text{FFT}\)，是快速計算 \(\text{DFT}\) 的方法，其時間複雜度僅為 \(O(n\log n)\)。

在 OI 中，通常只實現要求 \(n\) 為 \(2\) 的冪的 \(\text{FFT}\)。如果要計算 \(a\)和 \(b\) 的卷積 \(c\)，可以取 \(n\) 為不小於 \(c\) 的大小的最小的 \(2\) 的冪，將\(a\) 和 \(b\) 的大小均擴充為 \(n\) 後再計算。如果要計算迴圈卷積，可以先計算普通卷積，再將第 \(k+n\) 項加到第 \(k\) 項上。

#2.4 FFT の實現

設 \(n\) 為 \(2\) 的冪。對於 \(f(x)=\sum\limits_{k=0}^{n-1}a_kx^k\)，設

\[\begin{aligned} f_0(x)=&\sum\limits_{k=0}^{n/2-1}a_{2k}x^k\\ f_1(x)=&\sum\limits_{k=0}^{n/2-1}a_{2k+1}x^k \end{aligned} \]

那麼有 \(f_(x)=f_0(x^2)+xf_1(x^2).\)

當 \(0\leq s<\dfrac{n}{2}\) 時，有

\[\begin{aligned} f(\omega_n^s)&=f_0(\omega_n^{2s})+\omega_n^sf_1(\omega_n^{2s})\\ &=f_0(\omega_{\frac{n}{2}}^{s})+\omega_n^sf_1(\omega_{\frac{n}{2}}^{s}) \end{aligned} \]

顯然有 \(\omega_n^{\frac{n}{2}}=-1\)，於是當 \(\dfrac{n}{2}\leq s<n\) 時，有

\[\begin{aligned} f(\omega_n^s)=&f_0((\omega_n^s)^2)+\omega_n^sf_1((\omega_n^s)^2)\\ =&f_0((\omega_n^{s-\frac{n}{2}})^2)-\omega_n^{s-\frac{n}{2}}f_1((\omega_n^{s-\frac n 2})^2)\\ =&f_0(\omega_n^{2s-n})-\omega_n^{s-\frac{n}{2}}f_1(\omega_n^{2s-n})\\ =&f_0(\omega_{\frac n 2}^{s-\frac n 2})-\omega_n^{s-\frac n 2}f_1(\omega_\frac n 2^{s-\frac n 2}) \end{aligned} \]

也就是說，對於 \(0\leq s<\dfrac{n}{2}\)，有

\[\begin{aligned} f(\omega_n^s)&=f_0(\omega_{\frac{n}{2}}^{s})+\omega_n^sf_1(\omega_{\frac{n}{2}}^{s})\\ f(\omega_n^{s+\frac n 2})&=f_0(\omega_\frac n 2^s)-\omega_n^sf_1(\omega_\frac n 2^s) \end{aligned} \]

那麼只要遞迴地算出 \(f_0(x)\) 和 \(f_1(x)\) 的 \(\text{DFT}\)，就可以以 \(O(n)\) 的時間複雜度計算 \(f(x)\) 的 \(\text{DFT}\)，總時間複雜度為 \(O(n\log n)\)。

對於 \(n=2^k\)，可以將所有下標視為 \(k\) 位二進位制數。將下標根據是否小於 \(\dfrac{n}{2}\) 劃分成兩部分是根據最高位進行劃分，而根據奇偶性劃分成兩部分是根據最低位進行劃分。對於前者，可以自然地將遞迴樹展開。而對於後者，只要將所有下標的二進位制位翻轉，就可以轉化為前者的情況。

展開遞迴樹不僅會減少遞迴的開銷，還會去除冗餘的操作，並且優化了記憶體訪問效率。

#3.0 程式碼實現

#include <iostream>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <complex>
#define ll long long
#define mset(l,x) memset(l,x,sizeof(l))
using namespace std;

const int N = 10000010;
const int INF = 0x3fffffff;
const double PI = acos(-1);

typedef complex <double> cp;

int n,m,t = 1,res[N],ta[N],tb[N],lim,r[N];
cp ca[N],cb[N];

inline void fft(cp *a,int type){
    for (int i = 0;i < t;i ++)
      if (i < r[i]) swap(a[i],a[r[i]]);
    for (int mid = 1;mid < t;mid <<= 1){
        cp OMG(cos(PI / mid),type *sin(PI / mid));
        for (int R = mid << 1,j = 0;j < t;j += R){
            cp omg(1,0);
            for (int k = 0;k < mid;k ++,omg = omg * OMG){
                cp tmp1 = a[j + k],tmp2 = omg * a[j + mid + k];
                a[j + k] = tmp1 + tmp2;
                a[j + mid + k] = tmp1 - tmp2;
            }
        }
    }
}

int main(){
    scanf("%d%d",&n,&m);
    while (t <= n + m) t <<= 1,lim ++;
    for (int i = 0;i <= n;i ++){
        scanf("%d",&ta[i]);
        ca[i].real(ta[i]);
    }
    for (int i = 0;i <= m;i ++){
        scanf("%d",&tb[i]);
        cb[i].real(tb[i]);
    }
    for (int i = 0;i < t;i ++)
      r[i] = (r[i >> 1] >> 1) | ((i & 1) << (lim - 1));
    fft(ca,1);fft(cb,1);
    for (int i = 0;i < t;i ++)
      ca[i] *= cb[i];
    fft(ca,-1);
    for (int i = 0;i < t;i ++)
      res[i] = floor(ca[i].real() / t + 0.5);
    for (int i = 0;i <= n + m;i ++)
      cout << res[i] << " ";
    return 0;
}