Algorithm

In this section, we will go through the Circle FFT algorithm specifically, to interpolate a bivariate polynomial given the evaluations over a circle domain. We will also go over a concrete example which will help us understand the algorithm.

Circle FFT follows a divide-and-conquer strategy, as in the classical Cooley–Tukey FFT. We recursively reduce the task of interpolating a polynomial over some domain to interpolating a lower degree polynomial over a smaller domain. Thus at each recursive layer, we have "smaller" polynomials and their evaluations over "smaller" domains. Let us first go over this sequence of domains for the Circle FFT algorithm.

Sequence of Domains for Circle FFT

In Circle FFT, we use a 2-to-1 map to halve the domain size at each recursive layer. The domain used here is the circle domain $D_{n}$ of size $∣ D_{n} ∣ = 2^{n}$ .

$D_{n} = q + ⟨ g_{n - 1} ⟩ \cup - q + ⟨ g_{n - 1} ⟩$

This section describes two specific 2-to-1 maps that are central to the Circle FFT construction:

Projection map $π_{x}$ : The projection map $π_{x}$ projects the points $(x, y)$ and $(x, - y)$ to the same $x$ -coordinate i.e. $π_{x} : D_{n} \to S_{n}, π_{x} ((x, y)) = x$ where $S_{n}$ is the set containing all the $x$ -coordinates from $D_{n}$ . This can be interpreted as saying that two points are considered equivalent if they differ only by sign. Since $π_{x}$ maps two points from $D_{n}$ to a single element in $S_{n}$ , the size of $S_{n}$ will be half the size of $D_{n}$ .
Squaring map $π$ : The squaring map $π$ is a 2-to-1 map defined by: $π : S_{n} \to S_{n - 1}, π (x) = 2 x^{2} - 1$ This is obtained using the doubling map and the equality $y^{2} = 1 - x^{2}$ to compute the $x$ -coordinate: $π (x, y) = (x, y) + (x, y) = (x^{2} - y^{2}, 2 x y) = (2 x^{2} - 1, 2 x y)$

In Circle FFT, we use both the projection map $π_{x}$ and the squaring map $π$ to construct the sequence of domains. The sequence of domains for Circle FFT is shown as follows: $D_{n} π_{x} S_{n} π S_{n - 1} π \dots π S_{1}$

Now that we know the sequence of domains, let us dive into the Circle FFT algorithm.

Circle FFT

The Circle FFT interpolates a bivariate polynomial $p (x, y)$ from the polynomial space $L_{N} (F)$ , given the evaluations over a circle domain $D_{n}$ of size $N = 2^{n}$ . The algorithm is a three step process described as follows.

Step 1: Decompose $p (x, y)$ into sub-polynomials

In the first step, we decompose the bivariate polynomial $p (x, y)$ over $D_{n}$ using the projection map $π_{x}$ into two univariate polynomials $p_{0} (x)$ and $p_{1} (x)$ as follows: $p (x, y) = p_{0} (π_{x} (x, y)) + y \cdot p_{1} (π_{x} (x, y)) = p_{0} (x) + y \cdot p_{1} (x)$

Now to continue with the FFT algorithm, we want to compute the evaluations of $p_{0} (x)$ and $p_{1} (x)$ over the smaller domain $S_{n} = π_{x} (D_{n})$ , given the evaluations of $p (x, y)$ over $D_{n}$ . This is done using the following equations:

$p_{0} (x) = \frac{p ( x , y ) + p ( x , - y )}{2}, p_{1} (x) = \frac{p ( x , y ) - p ( x , - y )}{2 \cdot y}$

Substituting all evaluations of $p (x, y)$ over $D_{n}$ in the above equations, gives the evaluations of $p_{0} (x)$ and $p_{1} (x)$ over the domain $S_{n}$ .

Note

To compute the evaluations of $p_{1} (x)$ over the domain $S_{n}$ , we subtract the evaluations i.e. compute $p (x, y) - p (x, - y)$ and then divide by $2 \cdot y$ . These values $y$ are the $y$ -coordinates of the points in the circle domain $D_{n}$ . They are also referred to as circle twiddles and they only depend on the circle domain $D_{n}$ . Therefore they can be precomputed before the start of the FFT algorithm. We will look into these in detail in the next section.

Example: To make all the calculations easier we will work on a concrete example over the Mersenne prime $p = 2^{5} - 1 = 31$ . Thus all calculations are over $F_{31}$ , i.e, modulo $31$ . We are given the evaluations $v = [13, 16, 9, 30, 29, 27, 13, 21]$ of $p (x, y)$ over the circle domain $D_{3} = [(7, 18), (13, 7), (24, 13), (18, 24), (7, 13), (13, 24), (24, 18), (18, 7)]$ and we want to compute the coefficients of $p (x, y)$ i.e. interpolate $p (x, y)$ given its evaluations $v$ over $D_{3}$ .

Step 1: . First, decompose the polynomial $p (x, y)$ using the map $π_{x} (x, y) = x$ :

$p (x, y) = p_{0} (π_{x} (x, y)) + y \cdot p_{1} (π_{x} (x, y)) = p_{0} (x) + y \cdot p_{1} (x)$

Given the evaluations $v$ of $p (x, y)$ over the circle domain $D_{3}$ we aim to compute the evaluations of $p_{0}$ and $p_{1}$ over the smaller domain $S_{3}$ : $S_{3} = π_{x} (D_{3}) = [7, 13, 24, 18]$ For $(x, y) = (7, 18)$ and $(x, - y) = (7, 13)$ in $D_{3}$ , substituting them into the following equations give: $p_{0} (x) = \frac{p ( x , y ) + p ( x , - y )}{2} = \frac{p ( 7 , 18 ) + p ( 7 , 13 )}{2} = \frac{13 + 29}{2} = 21$ $p_{1} (x) = \frac{p ( x , y ) - p ( x , - y )}{2 \cdot y} = \frac{p ( 7 , 18 ) - p ( 7 , 13 )}{2 \cdot 18} = \frac{13 - 29}{2 \cdot 18} = 3$

Repeating this process for other pairs $(x, y)$ and $(x, - y)$ in $D_{3}$ , we obtain: $v_{0} = [21, 6, 11, 10], v_{1} = [3, 28, 7, 6]$

Step 2: Recursively apply FFT to sub-polynomials

Given the evaluations of $p_{0} (x)$ and $p_{1} (x)$ over $S_{n}$ , we now compute their coefficients.

This step mirrors the recursive structure of the Cooley–Tukey FFT. Each polynomial is recursively split using the squaring map $π (x) = 2 x^{2} - 1$ , and the process continues over successively smaller domains.

We begin by decomposing $p_{0} (x)$ using the squaring map $π$ as follows: $p_{0} (x) = p_{00} (π (x)) + x \cdot p_{01} (π (x))$

We compute the evaluations of $p_{00} (x)$ and $p_{01} (x)$ over $S_{n - 1} = π (S_{n})$ using the following equations.

$p_{00} (π (x)) = \frac{p _{0} ( x ) + p _{0} ( - x )}{2}, p_{01} (π (x)) = \frac{p _{0} ( x ) - p _{0} ( - x )}{2 \cdot x}$

This recursive process continues until we reach the base case: all evaluations of the polynomial over the domain are same. At that level, the coefficient is simply the evaluation of the constant polynomial.

Finally, we reconstruct the coefficients of $p_{0} (x)$ by working backward through the recursive calls, using the decomposition equation at each level. The same process applies to compute the coefficients of $p_{1} (x)$ .

Note

Similar to circle twiddles, to compute the evaluations of $p_{01} (x)$ over $S_{n - 1}$ we divide by the values $x$ which are the values from the domain $S_{n}$ . These values $x$ are referred to as line twiddles. For the next recursive layer, we divide by values from $S_{n - 1}$ i.e. $π (x)$ and for the next recursive layer we divide by $π^{2} (x)$ and so on. Thus the line twiddles is a vector of values $x, π (x), π^{2} (x), \dots$ and so on. These only depend on the initial domain $D_{n}$ and thus can be precomputed before the start of the algorithm.

Step 2: Given the evaluations of $p_{0}$ and $p_{1}$ over $S_{3} = [7, 13, 24, 18]$ : $v_{0} = [21, 6, 11, 10], v_{1} = [3, 28, 7, 6]$ we recursively apply the FFT algorithm to compute the coefficients of $p_{0}$ and $p_{1}$ .

At each layer, we decompose the polynomial using the polynomial map $π (x)$ . For example, the decomposition of $p_{0} (x)$ is:

$p_{0} (x) = p_{00} (π (x)) + x \cdot p_{01} (π (x))$

Omitting the intermediate steps of the recursive calls, we eventually obtain the coefficients of $p_{0} (x)$ and $p_{1} (x)$ as follows: $p_{0} (x) = 12 + 26 \cdot x + π (x) + 28 \cdot x \cdot π (x)$ $p_{1} (x) = 11 + 26 \cdot x + 14 \cdot π (x) + 20 \cdot x \cdot π (x)$

Step 3: Combine the coefficients

Finally, we combine the coefficients of $p_{0} (x)$ and $p_{1} (x)$ to compute the coefficients of the original bivariate polynomial $p (x, y)$ , using the decomposition: $p (x, y) = p_{0} (x) + y \cdot p_{1} (x)$

Step 3: Given the coefficients of $p_{0} (x)$ and $p_{1} (x)$ : $p_{0} (x) = 12 + 26 \cdot x + π (x) + 28 \cdot x \cdot π (x)$ $p_{1} (x) = 11 + 26 \cdot x + 14 \cdot π (x) + 20 \cdot x \cdot π (x)$ we reconstruct the original polynomial $p (x, y)$ using the decomposition: $p (x, y) = p_{0} (x) + y \cdot p_{1} (x)$

Substituting the expressions for $p_{0}$ and $p_{1}$ , we get: $p (x, y) = 12 + 26 \cdot x + π (x) + 28 \cdot x \cdot π (x) + y \cdot (11 + 26 \cdot x + 14 \cdot π (x) + 20 \cdot x \cdot π (x))$ $p (x, y) = 12 + 26 \cdot x + π (x) + 28 \cdot x \cdot π (x) +$ $11 \cdot y + 26 \cdot x \cdot y + 14 \cdot y \cdot π (x) + 20 \cdot x \cdot y \cdot π (x)$

This completes an overview of the interpolation algorithm using Circle FFT. In the next section, we will see how the twiddle values are computed and stored for Circle FFT.