Show how to do FFT by hand

Question

Say you have two polynomials: $3 + x$ and $2x^2 + 2$.

I'm trying to understand how FFT helps us multiply these two polynomials. However, I can't find any worked out examples. Can someone show me how FFT algorithm would multiply these two polynomials. (Note: there is nothing special about these polynomials, but I wanted to keep it simple to make it easier to follow.)

I've looked at the algorithms in pseudocode, but all of them seem to be have problems (don't specify what the input should be, undefined variables). And surprisingly, I can't find where anyone has actually walked through (by hand) an example of multiplying polynomials using FFT.

Answer 1

Suppose we use fourth roots of unity, which corresponds to substituting $1,i,-1,-i$ for $x$. We also use decimation-in-time rather than decimation-in-frequency in the FFT algorithm. (We also apply a bit-reversal operation seamlessly.)

In order to compute the transform of the first polynomial, we start by writing the coefficients: $$ 3,1,0,0. $$ The Fourier transform of the even coefficients $3,0$ is $3,3$, and of the odd coefficients $1,0$ is $1,1$. (This transform is just $a,b \mapsto a+b,a-b$.) Therefore the transform of the first polynomial is $$ 4,3+i,2,3-i. $$ This is obtained using $X_{0,2} = E_0 \pm O_0$, $X_{1,3} = E_1 \mp i O_1$. ( From twiddle factor calculation ).

Let's do the same for the second polynomial. The coefficients are $$2,0,2,0.$$ The even coefficients $2,2$ transform to $4,0$, and the odd coefficients $0,0$ transform to $0,0$. Therefore the transform of the second polynomial is $$ 4,0,4,0. $$

We obtain the Fourier transform of the product polynomial by multiplying the two Fourier transforms pointwise: $$ 16, 0, 8, 0. $$ It remains to compute the inverse Fourier transform. The even coefficients $16,8$ inverse-transform to $12,4$, and the odd coefficients $0,0$ inverse-transform to $0,0$. (The inverse transform is $x,y \mapsto (x+y)/2,(x-y)/2$.) Therefore the transform of the product polynomial is $$6,2,6,2.$$ This is obtained using $X_{0,2} = (E_0 \pm O_0)/2$, $X_{1,3} = (E_1 \mp i O_1)/2$. We have obtained the desired answer $$ (3 + x)(2 + 2x^2) = 6+2x+6x^2+2x^3. $$

Answer 2

Define the polynomials, where deg(A) = q and deg(B) = p. The deg(C) = q + p.

In this case, deg(C) = 1 + 2 = 3.

$$ A = 3 + x \\ B = 2x^2 + 2 \\ C = A*B = ? $$

We can easily find C in $O(n^2)$ time by brute-force multiplication of coefficients. By applying FFT (and inverse FFT), we could achieve this in $O(n\log(n))$ time. Explicitly:

1. Convert the coefficient representation of A and B to its value representation. This process is called evaluation. Performing Divide-and-Conquer (D&C) for this would take $O(n\log(n))$ time.
2. Multiply component-wise the polynomials in their value representation. This returns the value representation of C = A*B. This take $O(n)$ time.
3. Invert C using inverse FFT to get C in its coefficient representation. This process is called interpolation and it also takes $O(n\log(n))$ time.

Continuing along, we represent each polynomial as a vector whose value are its coefficients. We pad the vector with 0's up to the smallest power of two, $n = 2^k, n \geq deg(C) $. Thus $n = 4$. Choosing a power of two provides us a way to recursively apply our divide-and-conquer algorithm.

$$ \begin{align} &A = 3+x+0x^2+0x^3 \Rightarrow &\vec{a} = [3, 1, 0, 0]\\ &B = 2+0x+2x^+0x^3 \Rightarrow & \vec{b} = [2, 0, 2, 0] \end{align} $$

Let $A', B'$ be the value representation of A and B, respectively. Notice that FFT (Fast Fourier Transform) is a linear transformation (linear map) and can be represented as a matrix, $M$. Thus

$$ A' = M \overrightarrow{a} \\ B' = M \overrightarrow{b} $$

We define $M = M_n(\omega)$ where $\omega$ is complex roots $n^{th}$ complex roots of unity. Notice n = 4, in this example. Also notice that the entry in the $j^{th}$ row and $k^{th}$ column is $\omega_n^{jk}$ . See more about the DFT matrix here

$$ M_4(w) = \begin{bmatrix} 1 & 1 & 1 & ... & 1\\ 1 & \omega^1 & \omega^2 & ... & \omega^{n-1}\\ 1 & \omega^2 & \omega^{4} & ... & ... \\ ... & ... & ... & \omega^{jk} & ...\\ 1 & \omega^{n-1} & \omega^{2(n-1)} & ... & \omega^{(n-1)(n-1)} \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 &1 \\ 1 &\omega & \omega^2 & \omega^3\\ 1 &\omega^2 & \omega^4 & \omega^6\\ 1 &\omega^3 & \omega^6 & \omega^9 \end{bmatrix} $$

Given the $\omega_4 = 4^{th}$ roots of unity, we have the ordered set equality:

$$\{\omega^{0},\omega^1, \omega^2, \omega^3, \omega^4, \omega^5, ... \} = \{1, i, -1, -i, 1, i, ...\}$$

This can be visualized as iterating thru roots of the unit circle in the counter-clockwise direction.

Also, notice the mod n nature, i.e. $\omega^6 = \omega^{6 \mod n} = \omega^{2} = -1 $ and $-i = \omega^{3} = \omega^{3+n}$

To complete step 1 (evaluation) we find $A', B'$ by performing

$$ A' = M * \vec{a} = \begin{bmatrix} 1 & 1 & 1 &1 \\ 1 &\omega & \omega^2 & \omega^3\\ 1 &\omega^2 & \omega^4 & \omega^6\\ 1 &\omega^3 & \omega^6 & \omega^9 \end{bmatrix} \begin{bmatrix} 3 \\ 1 \\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix} 3 + 1 \\ 3 + 1 \omega \\ 3 + \omega^2 \\ 3 + \omega^3 \end{bmatrix} = \begin{bmatrix} 4 \\ 3 + i \\ 2 \\ 3 - i \end{bmatrix} \\ B' = M * \vec{b} = \begin{bmatrix} 1 & 1 & 1 &1 \\ 1 &\omega & \omega^2 & \omega^3\\ 1 &\omega^2 & \omega^4 & \omega^6\\ 1 &\omega^3 & \omega^6 & \omega^9 \end{bmatrix} \begin{bmatrix} 2 \\ 0 \\ 2 \\ 0 \end{bmatrix} = \begin{bmatrix} 2 + 2 \\ 2 + 2 \omega^2 \\ 2 + 2\omega^4 \\ 2 + 2\omega^6 \end{bmatrix} = \begin{bmatrix} 4 \\ 0 \\ 4 \\ 0 \end{bmatrix} $$

This step can be achieved using D&C algorithms (beyond the scope of this answer).

Multiply $A' * B'$ component-wise (step 2)

$$ A' * B' = \begin{bmatrix} 4 \\ 3 + i \\ 2 \\ 3 - i \end{bmatrix} \begin{bmatrix} 4 \\ 0 \\ 4 \\ 0 \end{bmatrix} = \begin{bmatrix} 16 \\ 0\\ 8 \\ 0 \end{bmatrix} = C'\\ $$

Finally, the last step is to represent C' into coefficients. Notice

$$ C' = M \vec{c} \\ \Rightarrow M^{-1}C' = M^{-1} M \vec{c} \\ \Rightarrow \vec{c} = M^{-1}C' $$

Notice $M_n^{-1} = \frac{1}{n} M_n(\omega^{-1})$¹ and $\omega^j = -\omega^{n/2 + j}$.

$$ M_n^{-1} = \frac{1}{4} \begin{bmatrix} 1 & 1 & 1 &1 \\ 1 &\omega^{-1} & \omega^{-2} & \omega^{-3}\\ 1 &\omega^{-2} & \omega^{-4} & \omega^{-6}\\ 1 &\omega^{-3} & \omega^{-6} & \omega^{-9} \end{bmatrix} = \frac{1}{4} \begin{bmatrix} 1 & 1 & 1 & 1\\ 1 & -i & -1 & i \\ 1 & -1 & 1 & -1\\ 1 & i & -1 & -i \end{bmatrix} $$

$\omega^{-j}$ can be visualized as iterating thru roots of the unit circle in the clockwise direction.

$$\{\omega^{0},\omega^{-1}, \omega^{-2}, \omega^{-3}, \omega^{-4}, \omega^{-5}, ... \} = \{1, -i, -1, i, 1, -i, ...\}$$

Also, it is true that, given the $n^{th}$ root of unity, the equality $\omega^{-j} = \omega^{n-j}$ holds. (Do you see why?)

Then, $$ \vec{c} = M^{-1}C' = \frac{1}{n} M_n(w^{-1}) = \frac{1}{4} \begin{bmatrix} 1 & 1 & 1 & 1\\ 1 & -i & -1 & i \\ 1 & -1 & 1 & -1\\ 1 & i & -1 & -i \end{bmatrix} \begin{bmatrix} 16 \\ 0\\ 8 \\ 0 \end{bmatrix} = \begin{bmatrix} (16 + 8) /4 \\ (16 - 8) /4 \\ (16 + 8) /4 \\ (16 - 8) /4 \end{bmatrix} = \begin{bmatrix} 6 \\ 2\\ 6 \\ 2 \end{bmatrix} $$

Thus, we get the polynomial $$ C = A * B = 6 + 2x + 6x^2 + 2x^3$$ ¹: Inversion Formula pg 73, Algorithms by Dasgupta et. al. (C) 2006