This section describes how the Cooley-Tukey fast Fourier transform works. As we learned in the previous section, the key is to select evaluation points that yield an efficient FFT algorithm.
Specifically, say we have such that , and for any .
Put another way, all the values are distinct and .
Put yet another way, the group generated by inside (written ) has size .
We call such an a primitive -th root of unity.
Suppose we have an which is a primitive th root of unity and let .
The FFT algorithm will let us compute for this set.
Actually, it is easier to see how it will let us compute the algorithm efficiently.
We will describe an algorithm that takes as input
- a primitive th root of unity
- in dense coefficients form (i.e., as a vector of coefficients of length ).
and outputs the vector of evaluations and does it in time (which is to say, if ).
Notice that naively, computing each evaluation using the coefficients of would require time , and so computing all of them would require time .
The algorithm can be defined recursively as follows.
If , then is a primitive st root of unity, and is a polynomial of degree . That means and also is a constant . So, we can immediately output the array of evaluations .
If , then we will split into two polynomials, recursively call on them, and reconstruct the result from the recursive calls.
To that end, define to be the polynomial whose coefficients are all the even-index coefficients of and the polynomial whose coefficients are all the odd-index coefficients of . In terms of the array representation, this just means splitting out every other entry into two arrays. So that can be done in time .
Write , so that and . Then
Now, notice that if is a th root of unity, then is a th root of unity. Thus we can recurse with and similarly for . Let
By assumption . So, for any we have
Now, since may be larger than , we need to reduce it mod , relying on the fact that if is an th root of unity then since . Thus, and so we have
We can compute the array in time (since each entry is the previous entry times ). Then we can compute each entry of the output in as
There are such entries, so this takes time .
This concludes the recursive definition of the algorithm .
- the coefficients of polynomial
- the even coefficients of corresponding to
- the odd coefficients of corresponding to
Now let’s analyze the time complexity. Let be the complexity on an instance of size (that is, for ).
Looking back at what we have done, we have done
- for computing and
- two recursive calls, each of size
- for computing the powers of
- for combining the results of the recursive calls
In total, this is . Solving this recurrence yields . Basically, there are recursions before we hit the base case, and each step takes time .
Now, in practice there are ways to describe this algorithm non-recursively that have better concrete performance, but that’s out of scope for this document. Read the code if you are interested.
So far we have a fast way to compute all at once, where is the set of powers of a th root of unity . For convenience let .
Now we want to go the other way and compute a polynomial given an array of evaluations. Specifically, evaluations uniquely define a degree polynomial. This can be written as a system of equations
which can be rewritten as a matrix vector product. This matrix is a Vandermonde matrix and it just so happens that square Vandermonde matrices are invertible, iff the are unique. Since we purposely selected our to be the powers of , a primitive -th root of unity, by definition are unique.
Therefore, to compute the polynomial given the corresponding array of evaluations (i.e. interpolation) we can solve for the polynomial’s coefficients using the inverse of the matrix. All we need now is the inverse of this matrix, which is slightly complicated to compute. I’m going to skip it for now, but if you have the details please make a pull request.
Substituting in the inverse matrix we obtain the equation for interpolation. Observe that this equation is nearly identical to the original equation for evaluation, except with the following substitution. Consequently and perhaps surprisingly, we can reuse the FFT algorithm in order to compute the inverse– .
So, suppose we have an array of field elements (which you can think of as a function ) and we want to compute the coefficients of a polynomial with .
To this end, define a polynomial by . That is, the polynomial whose coefficients are the evaluations in our array that we’re hoping to interpolate.
Now, let .
That is, we’re going to feed into the FFT algorithm defined above with as the th root of unity. It is not hard to check that if is an n-th root of unity, so is . Remember: the resulting values are the evaluations of on the powers of , so .
Now, let . That is, re-interpret the values returned by the FFT as the coefficients of a polynomial. I claim that is almost the polynomial we are looking for. Let’s calculate what values takes on at the powers of .
Now, let’s examine the quantity . We claim that if , then , and if , then . The first claim is clear since
For the second claim, we will prove that . This implies that . So either or . The former cannot be the case since it implies which in turn implies which is impossible since we are in the case . Thus we have as desired.
So let’s show that is invariant under multiplication by . Basically, it will come down to the fact that .
So now we know that
So if we define , then for every as desired. Thus we have our interpolation algorithm, sometimes called an inverse FFT or IFFT:
Input: the points we want to interpolate and a th root of unity.
Interpret the input array as the coefficients of a polynomial .
Output the polynomial . I.e., in terms of the dense-coefficients form, output the vector .
Note that this algorithm also takes time
Polynomials can be represented as a list of coefficients or a list of evaluations on a set
If the set is the set of powers of a root of unity, there are time algorithms for converting back and forth between those two representations
In evaluations form, polynomials can be added and multiplied in time
- TODO: caveat about hitting degree
Evaluations<F: FftField>that wrap a
Vec<F>and implement the FFT algorithms described above for converting between them
Familiarize yourself with the types and functions provided by