Convolution
Consider two polynomials $\gdef\mydef{\gdef}\mydef\s#1{\left(#1\right)}\mydef\m#1{\left[#1\right]}\mydef\b#1{\left\{#1\right\}}\mydef\abs#1{\left\lvert#1\right\rvert} F\s x = \sum_kf_kx^k$, $G\s x = \sum_kg_kx^k$ and their product $H\s x = \s{F \cdot G}\s x = \sum_kh_kx^k$, where
$$ h_k = \sum_{i + j = k}f_ig_j. \tag*{} $$
Toeplitz Matrix
We can describe convolution using the language of linear algebra.
Embed the Toeplitz matrix into a circulant matrix $C$:
$$ C_{ij} = t_{i - j} = \begin{bmatrix} t_0 & t_{-1} & t_{-2} & \cdots & \cdots & t_{-\s{n - 1}} \\ t_1 & t_0 & t_{-1} & \ddots & & \vdots \\ t_2 & t_1 & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & t_{-1} & t_{-2} \\ \vdots & & \ddots & t_1 & t_0 & t_{-1} \\ t_{n - 1} & \cdots & \cdots & t_2 & t_1 & t_0 \end{bmatrix}, \tag*{} $$
where $t_k = t_{k \bmod n}$, so that $C$ is circulant. Then consider $\mathbf y = \mathbf t * \mathbf x$:
$$ y_i = \sum_kt_kx_{i - k} = \sum_jt_{i - j}x_j = \sum_jC_{ij}x_j \tag*{}. $$
Therefore $\mathbf y = \mathbf t * \mathbf x = C\mathbf x$.
We want to find a way to compute $C\mathbf x$ efficiently, so we want to diagonalize $C$.
Diagonalization
Define cyclic shift operator $S$:
$$ S = \begin{bmatrix} 0 & 1 & 0 & \cdots \\ 0 & 0 & 1 & \cdots \\ \vdots & \vdots & \vdots & \ddots \\ 1 & 0 & 0 & \cdots \end{bmatrix}. \tag*{} $$
Then
$$ C = t_0I + t_1S + t_2S^2 + \cdots + t_{n - 1}S^{n - 1}. \tag*{} $$
Since each circulant matrix is a polynomial in $S$, diagonalizing $S$ simultaneously diagonalizes all circulant matrices. We want to find a basis $\b{\mathbf v_k}$ such that
$$ S\mathbf v_k = \lambda_k\mathbf v_k, \tag*{} $$
then
$$ C\mathbf v_k = \s{\sum_jt_j\lambda_k^j}\mathbf v_k. \tag*{} $$
The Fourier basis consists of eigenvectors of $S$:
$$ \mathbf v_k = \s{1, \omega^k_n, \omega^{2k}_n, \ldots, \omega^{\s{n - 1}k}_n}^{\intercal}, \tag*{} $$
Where $\omega_n = e^{2\pi i/n}$.
Discrete Fourier Transform
Define operator
$$ C_{\mathbf a}\s{\mathbf b} := \mathbf a * \mathbf b. \tag*{} $$
By assembling the Fourier basis vectors into the discrete Fourier transform matrix, we obtain a matrix that diagonalizes $S$.
$$ F = \left[\omega_n^{jk}\right]_{0 \leq j, k < n} = \begin{bmatrix} \omega_n^0 & \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0 \\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n - 1}\\ \omega_n^0 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2\s{n - 1}} \\ \omega_n^0 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3\s{n - 1}} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ \omega_n^0 & \omega_n^{n - 1} & \omega_n^{2\s{n - 1}} & \omega_n^{3\s{n - 1}} & \cdots & \omega_n^{\s{n - 1}^2} \end{bmatrix}. \tag*{} $$
We have
$$ F^{-1}SF = \operatorname{diag}\s{\omega_n^0, \omega_n^1, \ldots, \omega_n^{n - 1}}, \tag*{} $$
so
$$ F^{-1}C_{\mathbf a}F = \sum_ja_j\s{F^{-1}SF}^j = \operatorname{diag}\s{\sum_ja_j\omega_n^{jk}}_{k = 0}^{n - 1}. \tag*{} $$
Let
$$ \hat a_k = \sum_{j = 0}^{n - 1}a_j\omega_n^{jk}. \tag*{} $$
Let $A\s x = \sum_{j = 0}^{n - 1}a_jx^j$.
$$ F^{-1}C_{\mathbf a}F = \operatorname{diag}\s{\hat{\mathbf a}} = \operatorname{diag}\s{A\s{\omega^k_n}}_{k = 0}^{n - 1}. \tag*{} $$
Let $\mathcal F\b{\mathbf a} = F\mathbf a = \hat{\mathbf a}$.
Back to our original problem: polynomial multiplication.
Since $\s{A \cdot B}\s x = A\s x \cdot B\s x$, by multiplying each element of one vector by the corresponding element of the other vector:
$$ \mathcal F\b{\mathbf a * \mathbf b} = \mathcal F\b{\mathbf a} \cdot \mathcal F\b{\mathbf b}. \tag*{} $$
Finally, applying the inverse DFT, we obtain:
$$ A \cdot B = \mathcal F^{-1}\b{\mathcal F\b{\mathbf a} \cdot \mathcal F\b{\mathbf b}}. \tag*{} $$
At this point, we already know how to use DFT to compute polynomial multiplication. The remaining problem is to apply DFT efficiently.
Fast Fourier Transform
Assume $n$ is a power of $2$.
$$ F_n = \underbrace{ \begin{bmatrix} I & D \\ I & -D \end{bmatrix} }_{\text{butterfly}} \cdot \underbrace{ \begin{bmatrix} F_{n / 2} & 0 \\ 0 & F_{n / 2} \end{bmatrix}}_{\text{recursion}} \cdot \underbrace{P}_{\text{rearrange}}, \tag*{} $$
where
$$ D = \operatorname{diag}\s{1, \omega_n^1, \omega_n^2, \ldots, \omega_n^{n/2 - 1}}, \tag*{} $$
and $P$ is the permutation matrix that maps
$$ \s{a_0, a_1, \ldots, a_{n - 1}} \mapsto \s{a_0, a_2, a_4, \ldots, a_1, a_3, a_5, \ldots}. \tag*{} $$
Let
$$ A\s x = \sum_{k = 0}^{n - 1}a_kx^k. \tag*{} $$
Split it into two polynomials:
$$ A_0\s x = \sum_{k = 0}^{n/2 - 1}a_{2k}x^k = a_0 + a_2x + a_4x^2 + \cdots + a_{n - 2}x^{n/2 - 1}, \tag*{} $$
$$ A_1\s x = \sum_{k = 0}^{n/2 - 1}a_{2k + 1}x^k = a_1 + a_3x + a_5x^2 + \cdots + a_{n - 1}x^{n/2 - 1}, \tag*{} $$
so that
$$ A\s x = A_0\s{x^2} + xA_1\s{x^2}. \tag*{} $$
Pick $x = \omega_n^k$ in $A\s x$:
$$ \begin{aligned} A\s{\omega_n^k} &= A_0\s{\s{\omega_n^k}^2} + \omega_n^kA_1\s{\s{\omega_n^k}^2} \\ &= A_0\s{\omega_{n/2}^k} + \omega_n^kA_1\s{\omega_{n/2}^k}. \end{aligned} \tag{$\dagger$} $$
And
$$ \begin{aligned} A\s{\omega_n^{k + n/2}} = A\s{-\omega_n^k} &= A_0\s{\s{-\omega_n^k}^2} - \omega_n^kA_1\s{\s{-\omega_n^k}^2} \\ &= A_0\s{\omega_{n/2}^k} - \omega_n^kA_1\s{\omega_{n/2}^k}. \end{aligned} \tag{$\ddagger$} $$
Therefore, after applying DFT to $A_0$ and $A_1$, use $\dagger$ to compute $A\s{\omega_n^k}$ for $0 \leq k < n/2$ and use $\ddagger$ for $n/2 \leq k < n$.
We have
$$ T_{\text{DFT}}\s n = 2T_{\text{DFT}}\s{n/2} + O\s n. \tag*{} $$
Through master theorem, $T_{\text{DFT}}\s n = O\s{n \log n}$.
Inverse DFT
At this point, we have
$$ \hat{\mathbf a} = \mathcal F\b{\mathbf a} = F\mathbf a. \tag*{} $$
Hence
$$ \mathbf a = \mathcal F^{-1}\b{\hat{\mathbf a}} = F^{-1}\hat{\mathbf a}. \tag*{} $$
A quick check can verify that the inverse of the matrix $F$ has the following form:
$$ F^{-1} = \frac1n\left[\omega_n^{-jk}\right]_{0 \leq j, k < n} = \frac1n \begin{bmatrix} \omega_n^0 & \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0 \\ \omega_n^0 & \omega_n^{-1} & \omega_n^{-2} & \omega_n^{-3} & \cdots & \omega_n^{-\s{n - 1}}\\ \omega_n^0 & \omega_n^{-2} & \omega_n^{-4} & \omega_n^{-6} & \cdots & \omega_n^{-2\s{n - 1}} \\ \omega_n^0 & \omega_n^{-3} & \omega_n^{-6} & \omega_n^{-9} & \cdots & \omega_n^{-3\s{n - 1}} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ \omega_n^0 & \omega_n^{-\s{n - 1}} & \omega_n^{-2\s{n - 1}} & \omega_n^{-3\s{n - 1}} & \cdots & \omega_n^{-\s{n - 1}^2} \end{bmatrix}. \tag*{} $$
Note that these problems are almost the same as DFT, so the coefficients $a_k$ can be found by the same divide and conquer algorithm, as well as the direct FFT, only instead of $\omega_n^k$, we have to use $\omega_n^{-k}$, and at the end we need to divide the resulting coefficients by $n$.
- Categories: 学习笔记
- Tag: Polynomial FFT