The Fast Fourier Transform (FFT) developed by Cooley and Tukey in 1965 has its origins in the work of Gauss. The FFT, its variants and extensions to finite fields, are a fundamental algorithmic tool and a beautiful example of interplay between algebra and combinatorics. There are many great resources on FFT, see ingopedia’s curated list.
Given a field $F$, a polynomial of degree at most $d-1$ is defined as $P(X)=a_0+a_1 X+ \dots + a_{d-1} X^{d-1}$, where each $a_i \in F$. $P(X)$ can be viewed as a vector $\vec{a}=\langle a_0,\dots,a_{d-1} \rangle$ of $d$ elements in $F$.
Polynomials can also be evaluated. Given $d$ unique points $x_0,\dots,x_{d-1}$ in $F$, we can represent $P(X)$ as the vector $\vec{y}=\langle y_0,\dots,y_{d-1} \rangle$, where $y_i=P(x_i)$ for $i=0,\ldots,d-1$. We use zero-based indexing throughout: $a_i$ denotes a coefficient, $x_i$ an evaluation point, and $y_i$ the corresponding evaluation.
In a previous post we showed the isomorphism $\vec{a} \leftrightarrow \vec{y}$ between degree at most $d-1$ polynomials and their evaluations at $d$ unique points. Here we focus on the computational cost of this mapping: the number of field operations.
Given $\vec{a}$ for a polynomial of degree at most $d-1$, one can compute each $y_i=P(x_i)$ using Horner’s method in $O(d)$ field operations, by writing $P(X)=a_0+X(a_1+X(a_2+\cdots+Xa_{d-1}))$. Computing all $d$ evaluations therefore takes $O(d^2)$ operations. Can we compute the evaluation $\vec{a}\to\vec{y}$ in better than quadratic time?
Conversely, given evaluations $y_0,\ldots,y_{d-1}$ at $d$ distinct points $x_0,\ldots,x_{d-1}$, Lagrange interpolation recovers
\[P(X)=\sum_{i=0}^{d-1}y_iL_i(X),\qquad L_i(X)=\frac{\prod_{j\ne i}(X-x_j)}{\prod_{j\ne i}(x_i-x_j)}.\]Evaluating a single $L_i$ at one point takes $O(d)$ operations using this product formula. To recover the entire coefficient vector of $P$ in $O(d^2)$ operations, we share work across the basis polynomials. First construct the coefficients of $Q(X)=\prod_{j=0}^{d-1}(X-x_j)$ by multiplying its linear factors one at a time. Multiplying a degree-$k$ polynomial by one linear factor takes $O(k)$ operations, so constructing $Q$ takes $O(d^2)$ operations in total. For each $i$, division by $X-x_i$ takes $O(d)$ operations, and computing the nonzero normalization factor $\prod_{j\ne i}(x_i-x_j)$ and scaling the quotient also take $O(d)$ operations. Thus all $d$ basis polynomials can be constructed in $O(d^2)$ operations. Summing their coefficient vectors with weights $y_i$ takes another $O(d^2)$ operations. Can we compute the interpolation $\vec{y}\to\vec{a}$ in better than quadratic time?
The Fast Fourier Transform allows interpolating and evaluating a polynomial on special points, roots of unity, with $O( d \log d)$ operations.
FFT Theorem:
- Given $d=2^D$ and a field $F$ of characteristic $p=qd+1$, let $\omega$ be a primitive $d$-th root of unity. There exists an algorithm that in $O(d \log d )$ operations:
- Computes $P(\omega^0),\dots,P(\omega^{d-1})$ given the coefficients of a degree at most $d-1$ polynomial $P(X)$.
- Computes the coefficients of the degree at most $d-1$ polynomial $P(X)$ given $d$ evaluations $P(\omega^0),\dots,P(\omega^{d-1})$.
This result uses two powerful ingredients: divide-and-conquer, and roots of unity. We start with the first direction: computing the evaluations given the coefficients.
The first ingredient: decompose a polynomial into its odd and even coefficients each with half the degree
Assume that $P(X)$ has degree at most $d-1$ and that $d=2^D$ is a power of two. A classic divide-and-conquer algorithm design principle when applied to polynomials is to “split” a degree at most $d-1$ polynomial into two degree at most $d/2-1$ polynomials by dividing the polynomial to its even and odd coefficients:
\[P(X)= P_e(X^2) + XP_o(X^2) \label{\star} \tag{*}\]Where $P_e(X),P_o(X)$ are polynomials induced by the even and odd coefficients of $P(X)$ as follows:
\[P_e(X)=a_0 + a_2 X+ a_4 X^2 + \dots+ a_{d-2}X^{d/2 -1}\] \[P_o(X)=a_1 + a_3 X+ a_5 X^2 + \dots+ a_{d-1}X^{d/2 -1}\]Here is a pictorial representation of the odd and even indexes:
Why is this such a powerful divide-and-conquer decomposition?
- Given $d$ evaluations $P_e(x_0^2),\ldots,P_e(x_{d-1}^2)$ and $d$ evaluations $P_o(x_0^2),\ldots,P_o(x_{d-1}^2)$, one can compute the desired $d$ evaluations $P(x_0),\ldots,P(x_{d-1})$ in $2d$ operations, simply because $P(x_i)=P_e(x_i^2)+x_iP_o(x_i^2)$ via $(\ref{\star})$.
- The degree of $P_e(X)$ and $P_o(X)$ is at most $d/2-1$. So if we manage to repeat this recursion we can get $O(\log d)$ depth!
- However, to continue recursively there seems to be a challenge:
- $P_e(X)$ and $P_o(X)$ are of degree less than $d/2$, so each recursive call gives us only $d/2$ evaluations at points $z_0,\ldots,z_{d/2-1}$.
- To compute all $d$ evaluations of $P(X)$ via $(\ref{\star})$, we need evaluations of both smaller polynomials at the $d$ squares $x_0^2,\ldots,x_{d-1}^2$.
- How can we generate these $d$ evaluations from just $d/2$ evaluations using $O(d)$ operations?
The second ingredient: use roots of unity
For $g\in F$, if $g^d=1$ then we call $g$ a $d$-th root of unity. Moreover, if $g^0,g^1,\ldots,g^{d-1}$ are $d$ distinct elements, then we call $g$ a primitive $d$-th root of unity.
Suppose $g$ is such a primitive root, with $d=2^D$. Its powers $1,g,\ldots,g^{d-1}$ form a cyclic group of order $d$.
Now square these roots. For $d=8$, we get:
\[\begin{array}{c|cccc|cccc} \text{root}&1&g&g^2&g^3&g^4&g^5&g^6&g^7\\ \hline \text{square}&1&g^2&g^4&g^6&1&g^2&g^4&g^6 \end{array}\]The second half repeats the first! We asked for evaluations at eight squares, but there are only four distinct points.
This happens for every power-of-two $d\geq2$: for $j=0,\ldots,d/2-1$,
\[(g^{j+d/2})^2=g^{2j+d}=g^{2j}=(g^j)^2.\]Thus each root in the first half has the same square as the corresponding root in the second half. The $d$ squares are just the $d/2$ distinct points
\[1,g^2,g^4,\ldots,g^{d-2},\]repeated twice. These points form the cyclic group of order $d/2$ generated by $g^2$.
So we can solve the challenge above: choose $x_j=g^j$ for $j=0,\ldots,d-1$, and recursively evaluate $P_e$ and $P_o$ at the $d/2$ points $(g^2)^j$, for $j=0,\ldots,d/2-1$. Repeating each resulting evaluation vector twice gives their evaluations at all $d$ squares $x_j^2$.
The $d/2$ evaluations at the $d/2$-th roots, repeated twice, give all $d$ evaluations at the squares of the $d$-th roots. This was the missing piece to complete the divide-and-conquer!
As a visual example, fix $d=8$, every square of a power of $\omega_8$, say $(\omega_8^i)^2$ equals a power of $\omega_4$ (in particular $\omega_4^{i \bmod 4}$). For a specific example consider the squares of $\omega_8^3$ and $\omega_8^7$. Observe that $(\omega_8^3)^2 = \omega_8^6 = \omega_4^3$ and similarly $(\omega_8^7)^2 = \omega_8^{14} = \omega_8^{6} = \omega_4^3$:
Finding a primitive root of unity
If $F$ is a prime Field of characteristic $p=dq+1$ then we are guaranteed that for any nonzero $x \in F$, $g=x^q$ is a $d$-th root of unity because $g^d = x^{qd}=x^{p-1}=1$ (using Fermat’s little theorem). It is easy to check that the $d$-th roots of unity form a multiplicative group, and hence a cyclic group.
Choosing $x$ uniformly at random from $F\setminus {0}$ we get that $x^q$ is a uniformly random $d$-th roots of unity.
To see why this is true, denote $f(x)=x^q$ and see that $|x \text{ s.t. } f(x)=y|=q$ for all $y$. Indeed $f(x)$ maps a set of size $dq$ to a set of size at most $d$ (since $f(x)^d=1$ is a root of $x^d -1$ which can have at most $d$ roots). Finally, for any $y$ there can be at most $q$ values such that $x^q=y$. So $f(x)$ must map exactly $q$ elements to each of the $d$ options. Here we used twice the fundamental fact that nontrivial polynomials cannot have more roots than their degree.
Since $d=2^D$, the cyclic group of order $d$ has $d/2$ primitive roots (all the odd powers are generators). Moreover, if $g=x^q$ is not a primitive root, then it’s an even power, so $g^{d/2} = 1$.
That implies a simple Las Vegas algorithm: Randomly choose a nonzero $x \in F$, set $g=x^q$, and if $g^{d/2} \neq 1$ then output $g$, otherwise repeat. Since the probability of success is $1/2$, the expected number of rounds is 2.
Computational cost of the FFT
Each recursive call halves the number of coefficients, so the depth is $\log_2 d = D$. At layer $i$, there are $2^i$ subproblems, each involving a polynomial of degree at most $d/2^i-1$. Combining the evaluations from the two children of each non-leaf subproblem takes $2d/2^i$ field operations. Thus each non-leaf layer takes $2^i \times (2d/2^i)=2d$ field operations.
So overall $2 d \log_2 d$ field operations.
Interpolation: inverse FFT
Given $d$ evaluations $y_j=P(\omega^j)$ for $j=0,\ldots,d-1$, can we interpolate and compute the coefficients $\vec{a}$ of $P(X)$ in $O(d\log d)$ operations? Yes, using the same two ingredients!
Set $m=d/2$ and recall
\[P(X)=P_e(X^2)+XP_o(X^2).\]If we had the $m$ evaluations
\[e_j=P_e((\omega^2)^j),\qquad o_j=P_o((\omega^2)^j), \qquad j=0,\ldots,m-1,\]then we could recursively obtain the coefficient vectors $\vec{b},\vec{c}$ of $P_e$ and $P_o$, respectively. By $(\ref{\star})$, we recover the original coefficients as $a_{2i}=b_i$ and $a_{2i+1}=c_i$, for $i=0,\ldots,m-1$.
So the only challenge is to compute $\vec{e}$ and $\vec{o}$ from $\vec{y}$. For each $j=0,\ldots,m-1$, evaluate $(\ref{\star})$ at $\omega^j$ and $\omega^{j+m}$:
\[y_j=P(\omega^j)=P_e(\omega^{2j})+\omega^jP_o(\omega^{2j}) =e_j+\omega^jo_j,\] \[y_{j+m}=P(\omega^{j+m})=P_e(\omega^{2j})-\omega^jP_o(\omega^{2j}) =e_j-\omega^jo_j.\]Here we use $(\omega^{j+m})^2=\omega^{2j}$ and $\omega^{j+m}=-\omega^j$. The latter follows from $\omega^m=-1$: since $(\omega^m)^2=1$, we have $(\omega^m-1)(\omega^m+1)=0$, and $\omega^m\ne1$ because $\omega$ has order $d=2m$.
These two evaluations contain the same even contribution and opposite odd contributions. Adding them cancels the odd contribution; subtracting them cancels the even contribution. We can therefore recover both:
\[e_j=\frac{y_j+y_{j+m}}{2},\qquad o_j=\frac{y_j-y_{j+m}}{2\omega^j}.\]Thus we can compute all $m$ pairs $e_j,o_j$ in $O(d)$ operations. Again there are $\log_2 d$ layers of recursion, each requiring $O(d)$ operations, for a total of $O(d\log d)$ operations.
Inverse FFT via Inversion
The FFT is a linear transformation of the Vandermonde matrix where $V_{i,j}=\omega^{ij}$ for $i,j=0,\ldots,d-1$. The inverse matrix is $U_{i,j}= \omega^{-ij}/d$. This can be seen as the $i,j$ value of the product $VU$ is $(VU){i,j} = \frac{1}{d} \sum{k=0}^{d-1} \omega^{k(i-j)}$ which equals 1 for $i=j$ and 0 otherwise (because $(\omega -1 )(1+\omega+\dots+\omega^{d-1}) = \omega^d -1 = 0$). Hence, running the regular FFT with $\omega^{-1}$ and dividing by $d$ is the same as multiplying by the inverse $U$. This obtains inverse FFT as well.
Do we need fields of characteristic $q 2^D+1$?
The requirement to have a field of characteristic $q 2^D+1$ is limiting. What if we are using a field that does not have a large power of two as a root of unity, but we still want to work with large, degree at most $d-1$ polynomials? Recent work of Ben-Sasson, Carmon, Kopparty, and Levit shows how to do FFT for any field in $O(d \log d)$ operations!
Acknowledgments
Many thanks to Gilad Stern, Valeria Nikolaenko, and Dan Boneh for insightful discussions and comments.
Your comments and thoughts on x.