Concrete FFT polynomial multiplication example

Question

I have read a number of explanations of the steps involved in multiplying two polynomials using fast fourier transform and am not quite getting it in practice. I was wondering if I could get some help with a concrete example such as:

$$ p(x) = a_0 + a_2x^2 + a_4x^4 + a_6x^6 $$ $$ q(x) = b_0 + b_4x^4 + b_6x^6 + b_8x^8 $$

Since the $p(x)$'s degree is not a power of 2, do I pad it with zero coefficient terms and evaluate it as a degree 8 polynomial as well? I noticed that these polynomials don't have values for the odd terms; does this still multiply in O(nlog(n)) time complexity? I was unsure how to split it up.

When do I begin evaluating each term at the 8th roots of unity? i.e., $e^{\frac{i\pi k}{8}}$ for $k \in 2,4,6$.

Essentially, I seem to understand each component of component of the fft multiplication when I read it but I am yet to see a step by step concrete example of its process.

Can someone outline the steps for the multiplication of the above polynomials (or a similar simple multiplication) using fft? It would help me a lot.

Answer 1

Let us take your concrete example:- $$p(x) = a_0 + a_2x^2 + a_4x^4 + a_6x^6$$ $$q(x) = b_0 + b_4x^4 + b_6x^6 + b_8x^8$$

We can express the above two polynomials in vector format:- $$\mathbf{p}=[a_0,0,a_2,0,a_4,0,a_6]$$ $$\mathbf{q}=[b_0,0,0,0,b_4,0,b_6,0,b_8]$$

where the $k$th entry (for $k\in\{0,1,2,..\}$) represents the coefficient of $x^k$.

Once we have the vectors, we need to pad $\mathbf{p}$ with $8$ zeros (the length of vector $\mathbf{q}$ minus $1$), and we need to pad $\mathbf{q}$ with $6$ zeros (the length of vector $\mathbf{p}$ minus $1$), so we have:-

$$\mathbf{p'}=[a_0,0,a_2,0,a_4,0,a_6,\color{red}{0,0,0,0,0,0,0,0}]$$ $$\mathbf{q'}=[b_0,0,0,0,b_4,0,b_6,0,b_8,\color{red}{0,0,0,0,0,0}]$$
The reason why we need zero padding is to ensure that when we carry out the FFT on the vectors, multiply the FFTs element wise and carry out an Inverse FFT (IFFT), the result will correspond to a linear convolution.

It is worth noting that the length of $\mathbf{p'}$ and $\mathbf{q'}$ is $15$. As the FFT is designed to optimally process an integer power of $2$ number of samples, we can pad $\mathbf{p'}$ and $\mathbf{q'}$ by an extra $0$ element, so that we have $16$ samples to process.

$$\mathbf{p''}=[a_0,0,a_2,0,a_4,0,a_6,\color{red}{0,0,0,0,0,0,0,0,0}]$$ $$\mathbf{q''}=[b_0,0,0,0,b_4,0,b_6,0,b_8,\color{red}{0,0,0,0,0,0,0}]$$

Next we carry out the FFT on vectors $\mathbf{p''}$ and $\mathbf{q''}$, and we denote the resulting vectors by $\mathbf{P}$ and $\mathbf{Q}$ respectively (each element will be a complex number):-

$$\mathbf{P}=[A_0,A_1,A_2,A_3,A_4,A_5,A_6,A_7,A_8,A_9,A_{10},A_{11},A_{12},A_{13},A_{14},A_{15}]$$ $$\mathbf{Q}=[B_0,B_1,B_2,B_3,B_4,B_5,B_6,B_7,B_8,B_9,B_{10},B_{11},B_{12},B_{13},B_{14},B_{15}]$$ We then multiply each element of $\mathbf{P}$ with the corresponding element of $\mathbf{Q}$, resulting in vector $$\mathbf{R}=[C_0,C_1,C_2,...,C_{15}]$$ where $C_k=A_k\times B_k$.

The final step is to calculate the IFFT of vector $\mathbf{R}$, resulting in $$\mathbf{r}=[r_0,r_1,r_2,...,r_{15}]$$ where element $r_k$ is the coefficient of $x^k$, for $\in{0,1,..,15}$. Note that $r_{15}$ should be $0$, as we zero padded both vectors by an extra $0$ term.

Assuming that the coefficients of the polynomials $p(x)$ and $q(x)$ were real, then the elements of $\mathbf{r}$ should be real (any deviations, in terms of a small imaginary component, could arise from the finite precision nature of how the FFT is implemented on a processor).

Given the nature of the polynomials $p(x)$ and $q(x)$, the elements $r_m$ for odd $m$ will be zero (or a very small number, also due to the finite precision effects).

Answer 2

You should read "Introduction to algorithms" [CLRS] it explains every process of the algorithm and proof the necesary lemma/theorems to derive it.

for example it answers your question:

"How are the roots of unity for each polynomial used to reduce the problem size?"

The basic idea is that it uses a divide and conquer approach (nlogn part) based on the halving and cancellation lemmas. The roots of unity are used only to get the polynomials in their point-value representation (any different points could be used to get such representation) but these allow to apply conveniently the theorems above and get the asymptotic complexity desired.

"do I pad it with zero coefficient terms and evaluate it as a degree 8 polynomial as well? ".

Taking 2 polynomials in their point-value representation and leting $x_k$ be the points you are going to evaluate both $ 0 \leq k \leq n$ will yield $n+1$ point-value pairs, but you need $2n + 1$ pairs to interpolate a unique polynomial of degree $2n$ (the degree of the product polynomial of both of degree $n$). So you must extend the point-value representations of the polynomials to a $2n + 1$ point-value pairs each.

"When do I begin evaluating each term at the 8th roots of unity?" https://i.stack.imgur.com/2gYE2.jpg This picture answers you "when" (line 4). and is not only evaluated the 8th roots of unity it's also evaluated the 4th roots of unity, 2nd roots of unity. (divide and conquer... O(nlgn)) the size of vector $\vec{a}$ is reduced by a half as you can notice.

"I noticed that these polynomials don't have values for the odd terms; does this still multiply in O(nlog(n)) time complexity? "

Keep in mind that this notation only states an upper bound for the algorithm, this means that if your polynomials have "some" terms the upper bound will be the same as if it has "all" the terms.

If you truly grasp the algorithm and theory behind this (fairly easy but long enough for a simple answer in here IMO) you would be able to "see" the process and answer your doubts such as "when" and "how", "do I ..." and "complexity".

Good luck and keep learning.