快速傅里葉變換

前置知識：複數。

本文將介紹一種算法，它支持在 \(O(n\log n)\) 的時間內計算兩個 \(n\) 度的多項式的乘法，比樸素的 \(O(n^2)\) 算法更高效。由於兩個整數的乘法也可以被當作多項式乘法，因此這個算法也可以用來加速大整數的乘法計算。

引入

我們現在引入兩個多項式 \(A\) 和 \(B\)：

\[ \begin{aligned} A ={}& 5x^2 + 3x + 7 \\ B ={}& 7x^2 + 2x + 1 \\ \end{aligned} \]

兩個多項式相乘的積 \(C = A \times B\)，我們可以在 \(O(n^2)\) 的時間複雜度中解得（這裏 \(n\) 為 \(A\) 或者 \(B\) 多項式的度）：

\[ \begin{aligned} C ={}& A \times B \\ ={}& 35x^4 + 31x^3 + 60x^2 + 17x + 7 \end{aligned} \]

很明顯，多項式 \(C\) 的係數 \(c_i\) 滿足 \(c_i = \sum_{j = 0}^i a_j b_{i - j}\)。而對於這種樸素算法而言，計算每一項的時間複雜度都為 \(O(n)\)，一共有 \(O(n)\) 項，那麼時間複雜度為 \(O(n^2)\)。

能否加速使得它的時間複雜度降低呢？如果使用快速傅里葉變換的話，那麼我們可以使得其複雜度降低到 \(O(n \log n)\)。

傅里葉變換

傅里葉變換（Fourier Transform）是一種分析信號的方法，它可分析信號的成分，也可用這些成分合成信號。許多波形可作為信號的成分，傅里葉變換用正弦波作為信號的成分。

設 \(f(t)\) 是關於時間 \(t\) 的函數，則傅里葉變換可以檢測頻率 \(\omega\) 的週期在 \(f(t)\) 出現的程度：

\[ F(\omega)=\mathbb{F}[f(t)]=\int_{-\infty}^{\infty}f(t)\mathrm{e}^{-\mathrm{i}{\omega}t}dt \]

它的逆變換是

\[ f(t)=\mathbb{F}^{-1}[F(\omega)]=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)\mathrm{e}^{\mathrm{i}{\omega}t}d\omega \]

逆變換的形式與正變換非常類似，分母 \(2\pi\) 恰好是指數函數的週期。

傅里葉變換相當於將時域的函數與週期為 \(2\pi\) 的復指數函數進行連續的內積。逆變換仍舊為一個內積。

傅里葉變換有相應的卷積定理，可以將時域的卷積轉化為頻域的乘積，也可以將頻域的卷積轉化為時域的乘積。

離散傅里葉變換

離散傅里葉變換（Discrete Fourier transform，DFT）是傅里葉變換在時域和頻域上都呈離散的形式，將信號的時域採樣變換為其 DTFT（discrete-time Fourier transform）的頻域採樣。

傅里葉變換是積分形式的連續的函數內積，離散傅里葉變換是求和形式的內積。

設 \(\{x_n\}_{n=0}^{N-1}\) 是某一滿足有限性條件的序列，它的離散傅里葉變換（DFT）為：

\[ X_k=\sum_{n=0}^{N-1}x_n\mathrm{e}^{-\mathrm{i}\frac{2\pi}{N}kn} \]

其中 \(\mathrm{e}\) 是自然對數的底數，\(i\) 是虛數單位。通常以符號 \(\mathcal {F}\) 表示這一變換，即

\[ \hat{x}=\mathcal{F}x \]

類似於積分形式，它的 逆離散傅里葉變換（IDFT）為：

\[ x_n=\frac{1}{N}\sum_{k=0}^{N-1}X_k\mathrm{e}^{\mathrm{i}\frac{2\pi}{N}kn} \]

可以記為：

\[ x=\mathcal{F}^{-1}\hat{x} \]

實際上，DFT 和 IDFT 變換式中和式前面的歸一化係數並不重要。在上面的定義中，DFT 和 IDFT 前的係數分別為 \(1\) 和 \(\frac {1}{N}\)。有時我們會將這兩個係數都改 \(\frac{1}{{\sqrt{N}}}\)。

離散傅里葉變換仍舊是時域到頻域的變換。由於求和形式的特殊性，可以有其他的解釋方法。

如果把序列 \(x_n\) 看作多項式 \(f(x)\) 的 \(x^n\) 項係數，則計算得到的 \(X_k\) 恰好是多項式 \(f(x)\) 代入單位根 \(\mathrm{e}^{\frac{-2\pi \mathrm{i}k}{N}}\) 的點值 \(f(\mathrm{e}^{\frac{-2\pi \mathrm{i}k}{N}})\)。

這便構成了卷積定理的另一種解釋辦法，即對多項式進行特殊的插值操作。離散傅里葉變換恰好是多項式在單位根處進行插值。

例如計算：

\[ \dbinom{n}{3}+\dbinom{n}{7}+\dbinom{n}{11}+\dbinom{n}{15}+\ldots \]

定義函數 \(f(x)\) 為：

\[ f(x)={(1+x)}^n=\dbinom{n}{0}x^0+\dbinom{n}{1}x^1+\dbinom{n}{2}x^2+\dbinom{n}{3}x^3+\ldots \]

然後可以發現，代入四次單位根 \(f(\mathrm{i})\) 得到這樣的序列：

\[ f(\mathrm{i})={(1+\mathrm{i})}^n=\dbinom{n}{0}+\dbinom{n}{1}\mathrm{i}-\dbinom{n}{2}-\dbinom{n}{3}\mathrm{i}+\ldots \]

於是下面的求和恰好可以把其餘各項消掉：

\[ f(1)+\mathrm{i}f(\mathrm{i})-f(-1)-\mathrm{i}f(-\mathrm{i})=4\dbinom{n}{3}+4\dbinom{n}{7}+4\dbinom{n}{11}+4\dbinom{n}{15}+\ldots \]

因此這道數學題的答案為：

\[ \dbinom{n}{3}+\dbinom{n}{7}+\dbinom{n}{11}+\dbinom{n}{15}+\ldots=\frac{2^n+\mathrm{i}(1+\mathrm{i})^n-\mathrm{i}(1-\mathrm{i})^n}{4} \]

這道數學題在單位根處插值，恰好構成離散傅里葉變換。

矩陣公式

由於離散傅立葉變換是一個 線性 算子，所以它可以用矩陣乘法來描述。在矩陣表示法中，離散傅立葉變換表示如下：

\[ \begin{bmatrix} X_{0} \\ X_{1} \\ X_{2} \\ \vdots \\ X_{N-1} \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \alpha & \alpha^{2} & \cdots & \alpha^{N-1} \\ 1 & \alpha^{2} & \alpha^{4} & \cdots & \alpha^{2(N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \alpha^{N-1} & \alpha^{2(N-1)} & \cdots & \alpha^{(N-1)(N-1)} \end{bmatrix} \begin{bmatrix} x_{0} \\ x_{1} \\ x_{2} \\ \vdots \\ x_{N-1} \end{bmatrix} \]

其中 \(\alpha = \mathrm{e}^{-\mathrm{i}\frac{2\pi}{N}}\)。

快速傅里葉變換

FFT 是一種高效實現 DFT 的算法，稱為快速傅立葉變換（Fast Fourier Transform，FFT）。它對傅里葉變換的理論並沒有新的發現，但是對於在計算機系統或者説數字系統中應用離散傅立葉變換，可以説是進了一大步。快速數論變換（NTT）是快速傅里葉變換（FFT）在數論基礎上的實現。

在 1965 年，Cooley 和 Tukey 發表了快速傅里葉變換算法。事實上 FFT 早在這之前就被發現過了，但是在當時現代計算機並未問世，人們沒有意識到 FFT 的重要性。一些調查者認為 FFT 是由 Runge 和 König 在 1924 年發現的。但事實上高斯早在 1805 年就發明了這個算法，但一直沒有發表。

分治法實現

FFT 算法的基本思想是分治。就 DFT 來説，它分治地來求當 \(x=\omega_n^k\) 的時候 \(f(x)\) 的值。基 - 2 FFT 的分治思想體現在將多項式分為奇次項和偶次項處理。

舉個例子，對於一共 \(8\) 項的多項式：

\[ f(x) = a_0 + a_1x + a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7 \]

按照次數的奇偶來分成兩組，然後右邊提出來一個 \(x\)：

\[ \begin{aligned} f(x) &= (a_0+a_2x^2+a_4x^4+a_6x^6) + (a_1x+a_3x^3+a_5x^5+a_7x^7)\\ &= (a_0+a_2x^2+a_4x^4+a_6x^6) + x(a_1+a_3x^2+a_5x^4+a_7x^6) \end{aligned} \]

分別用奇偶次次項數建立新的函數：

\[ \begin{aligned} G(x) &= a_0+a_2x+a_4x^2+a_6x^3\\ H(x) &= a_1+a_3x+a_5x^2+a_7x^3 \end{aligned} \]

那麼原來的 \(f(x)\) 用新函數表示為：

\[ f(x)=G\left(x^2\right) + x \times H\left(x^2\right) \]

利用偶數次單位根的性質 \(\omega^i_n = -\omega^{i + n/2}_n\)，和 \(G\left(x^2\right)\) 和 \(H\left(x^2\right)\) 是偶函數，我們知道在複平面上 \(\omega^i_n\) 和 \(\omega^{i+n/2}_n\) 的 \(G(x^2)\) 的 \(H(x^2)\) 對應的值相同。得到：

\[ \begin{aligned} f(\omega_n^k) &= G((\omega_n^k)^2) + \omega_n^k \times H((\omega_n^k)^2) \\ &= G(\omega_n^{2k}) + \omega_n^k \times H(\omega_n^{2k}) \\ &= G(\omega_{n/2}^k) + \omega_n^k \times H(\omega_{n/2}^k) \end{aligned} \]

和：

\[ \begin{aligned} f(\omega_n^{k+n/2}) &= G(\omega_n^{2k+n}) + \omega_n^{k+n/2} \times H(\omega_n^{2k+n}) \\ &= G(\omega_n^{2k}) - \omega_n^k \times H(\omega_n^{2k}) \\ &= G(\omega_{n/2}^k) - \omega_n^k \times H(\omega_{n/2}^k) \end{aligned} \]

因此我們求出了 \(G(\omega_{n/2}^k)\) 和 \(H(\omega_{n/2}^k)\) 後，就可以同時求出 \(f(\omega_n^k)\) 和 \(f(\omega_n^{k+n/2})\)。於是對 \(G\) 和 \(H\) 分別遞歸 DFT 即可。

考慮到分治 DFT 能處理的多項式長度只能是 \(2^m(m \in \mathbf{N}^ \ast )\)，否則在分治的時候左右不一樣長，右邊就取不到係數了。所以要在第一次 DFT 之前就把序列向上補成長度為 \(2^m(m \in \mathbf{N}^\ast )\)（高次係數補 \(0\)）、最高項次數為 \(2^m-1\) 的多項式。

在代入值的時候，因為要代入 \(n\) 個不同值，所以我們代入 \(\omega_n^0,\omega_n^1,\omega_n^2,\cdots, \omega_n^{n-1} (n=2^m(m \in \mathbf{N}^ \ast ))\) 一共 \(2^m\) 個不同值。

代碼實現方面，STL 提供了複數的模板，當然也可以手動實現。兩者區別在於，使用 STL 的 complex 可以調用 exp 函數求出 \(\omega_n\)。但事實上使用歐拉公式得到的虛數來求 \(\omega_n\) 也是等價的。

以上就是 FFT 算法中 DFT 的介紹，它將一個多項式從係數表示法變成了點值表示法。

值的注意的是，因為是單位復根，所以説我們需要令 \(n\) 項式的高位補為零，使得 \(n = 2 ^ k, k \in \mathbf{N}^ \ast\)。

遞歸版 FFT

#include <cmath>
#include <complex>

typedef std::complex<double> Comp;  // STL complex

const Comp I(0, 1);  // i
const int MAX_N = 1 << 20;

Comp tmp[MAX_N];

// rev=1,DFT; rev=-1,IDFT
void DFT(Comp* f, int n, int rev) {
  if (n == 1) return;
  for (int i = 0; i < n; ++i) tmp[i] = f[i];
  // 偶數放左邊，奇數放右邊
  for (int i = 0; i < n; ++i) {
    if (i & 1)
      f[n / 2 + i / 2] = tmp[i];
    else
      f[i / 2] = tmp[i];
  }
  Comp *g = f, *h = f + n / 2;
  // 遞歸 DFT
  DFT(g, n / 2, rev), DFT(h, n / 2, rev);
  // cur 是當前單位復根，對於 k = 0 而言，它對應的單位復根 omega^0_n = 1。
  // step 是兩個單位復根的差，即滿足 omega^k_n = step*omega^{k-1}*n，
  // 定義等價於 exp(I*(2*M_PI/n*rev))
  Comp cur(1, 0), step(cos(2 * M_PI / n), sin(2 * M_PI * rev / n));
  for (int k = 0; k < n / 2;
       ++k) {  // F(omega^k_n) = G(omega^k*{n/2}) + omega^k*n\*H(omega^k*{n/2})
    tmp[k] = g[k] + cur * h[k];
    // F(omega^{k+n/2}*n) = G(omega^k*{n/2}) - omega^k_n*H(omega^k\_{n/2})
    tmp[k + n / 2] = g[k] - cur * h[k];
    cur *= step;
  }
  for (int i = 0; i < n; ++i) f[i] = tmp[i];
}

時間複雜度 \(O(n\log n)\)。

倍增法實現

這個算法還可以從「分治」的角度繼續優化。對於基 - 2 FFT，我們每一次都會把整個多項式的奇數次項和偶數次項係數分開，一直分到只剩下一個係數。但是，這個遞歸的過程需要更多的內存。因此，我們可以先「模仿遞歸」把這些係數在原數組中「拆分」，然後再「倍增」地去合併這些算出來的值。

對於「拆分」，可以使用位逆序置換實現。

對於「合併」，使用蝶形運算優化可以做到只用 \(O(1)\) 的額外空間來完成。

位逆序置換

以 \(8\) 項多項式為例，模擬拆分的過程：

• 初始序列為 \(\{x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7\}\)
• 一次二分之後 \(\{x_0, x_2, x_4, x_6\},\{x_1, x_3, x_5, x_7 \}\)
• 兩次二分之後 \(\{x_0,x_4\} \{x_2, x_6\},\{x_1, x_5\},\{x_3, x_7 \}\)
• 三次二分之後 \(\{x_0\}\{x_4\}\{x_2\}\{x_6\}\{x_1\}\{x_5\}\{x_3\}\{x_7 \}\)

規律：其實就是原來的那個序列，每個數用二進制表示，然後把二進制翻轉對稱一下，就是最終那個位置的下標。比如 \(x_1\) 是 001，翻轉是 100，也就是 4，而且最後那個位置確實是 4。我們稱這個變換為位逆序置換（bit-reversal permutation），證明留給讀者自證。

根據它的定義，我們可以在 \(O(n\log n)\) 的時間內求出每個數變換後的結果：

位逆序置換實現（\(O(n\log n)\)）

/*
 * 進行 FFT 和 IFFT 前的反置變換
 * 位置 i 和 i 的二進制反轉後的位置互換
 * len 必須為 2 的冪
 */
void change(Complex y[], int len) {
  // 一開始 i 是 0...01，而 j 是 10...0，在二進制下相反對稱。
  // 之後 i 逐漸加一，而 j 依然維持着和 i 相反對稱，一直到 i = 1...11。
  for (int i = 1, j = len / 2, k; i < len - 1; i++) {
    // 交換互為小標反轉的元素，i < j 保證交換一次
    if (i < j) swap(y[i], y[j]);
    // i 做正常的 + 1，j 做反轉類型的 + 1，始終保持 i 和 j 是反轉的。
    // 這裏 k 代表了 0 出現的最高位。j 先減去高位的全為 1 的數字，直到遇到了
    // 0，之後再加上即可。
    k = len / 2;
    while (j >= k) {
      j = j - k;
      k = k / 2;
    }
    if (j < k) j += k;
  }
}

實際上，位逆序置換可以 \(O(n)\) 從小到大遞推實現，設 \(len=2^k\)，其中 \(k\) 表示二進制數的長度，設 \(R(x)\) 表示長度為 \(k\) 的二進制數 \(x\) 翻轉後的數（高位補 \(0\)）。我們要求的是 \(R(0),R(1),\cdots,R(n-1)\)。

首先 \(R(0)=0\)。

我們從小到大求 \(R(x)\)。因此在求 \(R(x)\) 時，\(R\left(\left\lfloor \dfrac{x}{2} \right\rfloor\right)\) 的值是已知的。因此我們把 \(x\) 右移一位（除以 \(2\)），然後翻轉，再右移一位，就得到了 \(x\) 除了（二進制）個位 之外其它位的翻轉結果。

考慮個位的翻轉結果：如果個位是 \(0\)，翻轉之後最高位就是 \(0\)。如果個位是 \(1\)，則翻轉後最高位是 \(1\)，因此還要加上 \(\dfrac{len}{2}=2^{k-1}\)。綜上

\[ R(x)=\left\lfloor \frac{R\left(\left\lfloor \frac{x}{2} \right\rfloor\right)}{2} \right\rfloor + (x\bmod 2)\times \frac{len}{2} \]

舉個例子：設 \(k=5\)，\(len=(100000)_2\)。為了翻轉 \((11001)_2\)：

1. 考慮 \((1100)_2\)，我們知道 \(R((1100)_2)=R((01100)_2)=(00110)_2\)，再右移一位就得到了 \((00011)_2\)。
2. 考慮個位，如果是 \(1\)，它就要翻轉到數的最高位，即翻轉數加上 \((10000)_2=2^{k-1}\)，如果是 \(0\) 則不用更改。

位逆序置換實現（\(O(n)\)）

// 同樣需要保證 len 是 2 的冪
// 記 rev[i] 為 i 翻轉後的值
void change(Complex y[], int len) {
  for (int i = 0; i < len; ++i) {
    rev[i] = rev[i >> 1] >> 1;
    if (i & 1) {  // 如果最後一位是 1，則翻轉成 len/2
      rev[i] |= len >> 1;
    }
  }
  for (int i = 0; i < len; ++i) {
    if (i < rev[i]) {  // 保證每對數只翻轉一次
      swap(y[i], y[rev[i]]);
    }
  }
  return;
}

蝶形運算優化

已知 \(G(\omega_{n/2}^k)\) 和 \(H(\omega_{n/2}^k)\) 後，需要使用下面兩個式子求出 \(f(\omega_n^k)\) 和 \(f(\omega_n^{k+n/2})\)：

\[ \begin{aligned} f(\omega_n^k) & = G(\omega_{n/2}^k) + \omega_n^k \times H(\omega_{n/2}^k) \\ f(\omega_n^{k+n/2}) & = G(\omega_{n/2}^k) - \omega_n^k \times H(\omega_{n/2}^k) \end{aligned} \]

使用位逆序置換後，對於給定的 \(n, k\)：

• \(G(\omega_{n/2}^k)\) 的值存儲在數組下標為 \(k\) 的位置，\(H(\omega_{n/2}^k)\) 的值存儲在數組下標為 \(k + \dfrac{n}{2}\) 的位置。
• \(f(\omega_n^k)\) 的值將存儲在數組下標為 \(k\) 的位置，\(f(\omega_n^{k+n/2})\) 的值將存儲在數組下標為 \(k + \dfrac{n}{2}\) 的位置。

因此可以直接在數組下標為 \(k\) 和 \(k + \frac{n}{2}\) 的位置進行覆寫，而不用開額外的數組保存值。此方法即稱為 蝶形運算，或更準確的，基 - 2 蝶形運算。

再詳細説明一下如何藉助蝶形運算完成所有段長度為 \(\frac{n}{2}\) 的合併操作：

1. 令段長度為 \(s = \frac{n}{2}\)；
2. 同時枚舉序列 \(\{G(\omega_{n/2}^k)\}\) 的左端點 \(l_g = 0, 2s, 4s, \cdots, N-2s\) 和序列 \(\{H(\omega_{n/2}^k)\}\) 的左端點 \(l_h = s, 3s, 5s, \cdots, N-s\)；
3. 合併兩個段時，枚舉 \(k = 0, 1, 2, \cdots, s-1\)，此時 \(G(\omega_{n/2}^k)\) 存儲在數組下標為 \(l_g + k\) 的位置，\(H(\omega_{n/2}^k)\) 存儲在數組下標為 \(l_h + k\) 的位置；
4. 使用蝶形運算求出 \(f(\omega_n^k)\) 和 \(f(\omega_n^{k+n/2})\)，然後直接在原位置覆寫。

快速傅里葉逆變換

傅里葉逆變換可以用傅里葉變換表示。對此我們有兩種理解方式。

線性代數角度

IDFT（傅里葉反變換）的作用，是把目標多項式的點值形式轉換成係數形式。而 DFT 本身是個線性變換，可以理解為將目標多項式當作向量，左乘一個矩陣得到變換後的向量，以模擬把單位復根代入多項式的過程：

\[ \begin{bmatrix}y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{bmatrix} = \begin{bmatrix}1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^1 & \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)^2} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{bmatrix} \]

現在我們已經得到最左邊的結果了，中間的 \(x\) 值在目標多項式的點值表示中也是一一對應的，所以，根據矩陣的基礎知識，我們只要在式子兩邊左乘中間那個大矩陣的逆矩陣就行了。

由於這個矩陣的元素非常特殊，它的逆矩陣也有特殊的性質，就是每一項 取倒數，再 除以變換的長度 \(n\)，就能得到它的逆矩陣。

注意：傅里葉變換的長度，並不是多項式的長度，變換的長度應比乘積多項式的長度長。待相乘的多項式不夠長，需要在高次項處補 \(0\)。

為了使計算的結果為原來的倒數，根據歐拉公式，可以得到

\[ \frac{1}{\omega_k}=\omega_k^{-1}=\mathrm{e}^{-\frac{2\pi \mathrm{i}}{k}}=\cos\left(\frac{2\pi}{k}\right)+\mathrm{i} \sin\left(-\frac{2\pi}{k}\right) \]

因此我們可以嘗試着把單位根 \(\omega_k\) 取成 \(\mathrm{e}^{-\frac{2\pi \mathrm{i}}{k}}\)，這樣我們的計算結果就會變成原來的倒數，之後唯一多的操作就只有再 除以它的長度 \(n\)，而其它的操作過程與 DFT 是完全相同的。我們可以定義一個函數，在裏面加一個參數 \(1\) 或者是 \(-1\)，然後把它乘到 \(\pi\) 上。傳入 \(1\) 就是 DFT，傳入 \(-1\) 就是 IDFT。

單位復根週期性

利用單位復根的週期性同樣可以理解 IDFT 與 DFT 之間的關係。

考慮原本的多項式是 \(f(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}=\sum_{i=0}^{n-1}a_ix^i\)。而 IDFT 就是把你的點值表示還原為係數表示。

考慮 構造法。我們已知 \(y_i=f\left( \omega_n^i \right),i\in\{0,1,\cdots,n-1\}\)，求 \(\{a_0,a_1,\cdots,a_{n-1}\}\)。構造多項式如下

\[ A(x)=\sum_{i=0}^{n-1}y_ix^i \]

相當於把 \(\{y_0,y_1,y_2,\cdots,y_{n-1}\}\) 當做多項式 \(A\) 的係數表示法。

這時我們有兩種推導方式，這對應了兩種實現方法。

方法一

設 \(b_i=\omega_n^{-i}\)，則多項式 \(A\) 在 \(x=b_0,b_1,\cdots,b_{n-1}\) 處的點值表示法為 \(\left\{ A(b_0),A(b_1),\cdots,A(b_{n-1}) \right\}\)。

對 \(A(x)\) 的定義式做一下變換，可以將 \(A(b_k)\) 表示為

\[ \begin{aligned} A(b_k)&=\sum_{i=0}^{n-1}f(\omega_n^i)\omega_n^{-ik}=\sum_{i=0}^{n-1}\omega_n^{-ik}\sum_{j=0}^{n-1}a_j(\omega_n^i)^{j}\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j\omega_n^{i(j-k)}=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\left(\omega_n^{j-k}\right)^i\\ \end{aligned} \]

記 \(S\left(\omega_n^a\right)=\sum_{i=0}^{n-1}\left(\omega_n^a\right)^i\)。

當 \(a=0 \pmod{n}\) 時，\(S\left(\omega_n^a\right)=n\)。

當 \(a\neq 0 \pmod{n}\) 時，我們錯位相減

\[ \begin{aligned} S\left(\omega_n^a\right)&=\sum_{i=0}^{n-1}\left(\omega_n^a\right)^i\\ \omega_n^a S\left(\omega_n^a\right)&=\sum_{i=1}^{n}\left(\omega_n^a\right)^i\\ S\left(\omega_n^a\right)&=\frac{\left(\omega_n^a\right)^n-\left(\omega_n^a\right)^0}{\omega_n^a-1}=0\\ \end{aligned} \]

也就是説

\[ S\left(\omega_n^a\right)= \begin{cases} n,&a=0\\ 0,&a\neq 0 \end{cases} \]

那麼代回原式

\[ A(b_k)=\sum_{j=0}^{n-1}a_jS\left(\omega_n^{j-k}\right)=a_k\cdot n \]

也就是説給定點 \(b_i=\omega_n^{-i}\)，則 \(A\) 的點值表示法為

\[ \begin{aligned} &\left\{ (b_0,A(b_0)),(b_1,A(b_1)),\cdots,(b_{n-1},A(b_{n-1})) \right\}\\ =&\left\{ (b_0,a_0\cdot n),(b_1,a_1\cdot n),\cdots,(b_{n-1},a_{n-1}\cdot n) \right\} \end{aligned} \]

綜上所述，我們取單位根為其倒數，對 \(\{y_0,y_1,y_2,\cdots,y_{n-1}\}\) 跑一遍 FFT，然後除以 \(n\) 即可得到 \(f(x)\) 的係數表示。

方法二

我們直接將 \(\omega_n^i\) 代入 \(A(x)\)。

推導的過程與方法一大同小異，最終我們得到 \(A(\omega_n^k) = \sum_{j=0}^{n-1}a_jS\left(\omega_n^{j+k}\right)\)。

當且僅當 \(j+k=0 \pmod{n}\) 時有 \(S\left(\omega_n^{j+k}\right) = n\)，否則為 \(0\)。因此 \(A(\omega_n^k) = a_{n-k}\cdot n\)。

這意味着我們將 \(\{y_0,y_1,y_2,\cdots,y_{n-1}\}\) 做 DFT 變換後除以 \(n\)，再反轉後 \(n - 1\) 個元素，同樣可以還原 \(f(x)\) 的係數表示。

代碼實現

所以我們 FFT 函數可以集 DFT 和 IDFT 於一身。代碼實現如下：

非遞歸版 FFT（對應方法一）

/*
 * 做 FFT
 * len 必須是 2^k 形式
 * on == 1 時是 DFT，on == -1 時是 IDFT
 */
void fft(Complex y[], int len, int on) {
  // 位逆序置換
  change(y, len);
  // 模擬合併過程，一開始，從長度為一合併到長度為二，一直合併到長度為 len。
  for (int h = 2; h <= len; h <<= 1) {
    // wn：當前單位復根的間隔：w^1_h
    Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
    // 合併，共 len / h 次。
    for (int j = 0; j < len; j += h) {
      // 計算當前單位復根，一開始是 1 = w^0_n，之後是以 wn 為間隔遞增： w^1_n
      // ...
      Complex w(1, 0);
      for (int k = j; k < j + h / 2; k++) {
        // 左側部分和右側是子問題的解
        Complex u = y[k];
        Complex t = w * y[k + h / 2];
        // 這就是把兩部分分治的結果加起來
        y[k] = u + t;
        y[k + h / 2] = u - t;
        // 後半個 「step」 中的ω一定和 「前半個」 中的成相反數
        // 「紅圈」上的點轉一整圈「轉回來」，轉半圈正好轉成相反數
        // 一個數相反數的平方與這個數自身的平方相等
        w = w * wn;
      }
    }
  }
  // 如果是 IDFT，它的逆矩陣的每一個元素不只是原元素取倒數，還要除以長度 len。
  if (on == -1) {
    for (int i = 0; i < len; i++) {
      y[i].x /= len;
    }
  }
}

非遞歸版 FFT（對應方法二）

/*
 * 做 FFT
 * len 必須是 2^k 形式
 * on == 1 時是 DFT，on == -1 時是 IDFT
 */
void fft(Complex y[], int len, int on) {
  change(y, len);
  for (int h = 2; h <= len; h <<= 1) {             // 模擬合併過程
    Complex wn(cos(2 * PI / h), sin(2 * PI / h));  // 計算當前單位復根
    for (int j = 0; j < len; j += h) {
      Complex w(1, 0);  // 計算當前單位復根
      for (int k = j; k < j + h / 2; k++) {
        Complex u = y[k];
        Complex t = w * y[k + h / 2];
        y[k] = u + t;  // 這就是把兩部分分治的結果加起來
        y[k + h / 2] = u - t;
        // 後半個 「step」 中的ω一定和 「前半個」 中的成相反數
        // 「紅圈」上的點轉一整圈「轉回來」，轉半圈正好轉成相反數
        // 一個數相反數的平方與這個數自身的平方相等
        w = w * wn;
      }
    }
  }
  if (on == -1) {
    reverse(y + 1, y + len);
    for (int i = 0; i < len; i++) {
      y[i].x /= len;
    }
  }
}

FFT 模板（HDU 1402 - A * B Problem Plus）

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>

const double PI = acos(-1.0);

struct Complex {
  double x, y;

  Complex(double _x = 0.0, double _y = 0.0) {
    x = _x;
    y = _y;
  }

  Complex operator-(const Complex &b) const {
    return Complex(x - b.x, y - b.y);
  }

  Complex operator+(const Complex &b) const {
    return Complex(x + b.x, y + b.y);
  }

  Complex operator*(const Complex &b) const {
    return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
  }
};

/*
 * 进行 FFT 和 IFFT 前的反置变换
 * 位置 i 和 i 的二进制反转后的位置互换
 *len 必须为 2 的幂
 */
void change(Complex y[], int len) {
  int i, j, k;

  for (int i = 1, j = len / 2; i < len - 1; i++) {
    if (i < j) std::swap(y[i], y[j]);

    // 交换互为小标反转的元素，i<j 保证交换一次
    // i 做正常的 + 1，j 做反转类型的 + 1，始终保持 i 和 j 是反转的
    k = len / 2;

    while (j >= k) {
      j = j - k;
      k = k / 2;
    }

    if (j < k) j += k;
  }
}

/*
 * 做 FFT
 *len 必须是 2^k 形式
 *on == 1 时是 DFT，on == -1 时是 IDFT
 */
void fft(Complex y[], int len, int on) {
  change(y, len);

  for (int h = 2; h <= len; h <<= 1) {
    Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));

    for (int j = 0; j < len; j += h) {
      Complex w(1, 0);

      for (int k = j; k < j + h / 2; k++) {
        Complex u = y[k];
        Complex t = w * y[k + h / 2];
        y[k] = u + t;
        y[k + h / 2] = u - t;
        w = w * wn;
      }
    }
  }

  if (on == -1) {
    for (int i = 0; i < len; i++) {
      y[i].x /= len;
    }
  }
}

const int MAXN = 200020;
Complex x1[MAXN], x2[MAXN];
char str1[MAXN / 2], str2[MAXN / 2];
int sum[MAXN];

int main() {
  while (scanf("%s%s", str1, str2) == 2) {
    int len1 = strlen(str1);
    int len2 = strlen(str2);
    int len = 1;

    while (len < len1 * 2 || len < len2 * 2) len <<= 1;

    for (int i = 0; i < len1; i++) x1[i] = Complex(str1[len1 - 1 - i] - '0', 0);

    for (int i = len1; i < len; i++) x1[i] = Complex(0, 0);

    for (int i = 0; i < len2; i++) x2[i] = Complex(str2[len2 - 1 - i] - '0', 0);

    for (int i = len2; i < len; i++) x2[i] = Complex(0, 0);

    fft(x1, len, 1);
    fft(x2, len, 1);

    for (int i = 0; i < len; i++) x1[i] = x1[i] * x2[i];

    fft(x1, len, -1);

    for (int i = 0; i < len; i++) sum[i] = int(x1[i].x + 0.5);

    for (int i = 0; i < len; i++) {
      sum[i + 1] += sum[i] / 10;
      sum[i] %= 10;
    }

    len = len1 + len2 - 1;

    while (sum[len] == 0 && len > 0) len--;

    for (int i = len; i >= 0; i--) printf("%c", sum[i] + '0');

    printf("\n");
  }

  return 0;
}

參考文獻

1. 桃醬的算法筆記.

---

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用