​

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 + \langle g_{n-1} \rangle \cup -q + \langle g_{n-1} \rangle$ This section describes two specific 2-to-1 maps that are central to the Circle FFT construction:

1. Projection map $\pi_x$: The projection map $\pi_x$ projects the points $(x,y)$ and $(x,-y)$ to the same $x$-coordinate i.e. $\pi_x: D_n \rightarrow S_n, \quad \pi_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 $\pi_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$.
2. Squaring map $\pi$: The squaring map $\pi$ is a 2-to-1 map defined by: $\pi: S_n \rightarrow S_{n-1}, \quad \pi(x) = 2x^2 - 1$ This is obtained using the doubling map and the equality $y^2 = 1 - x^2$ to compute the $x$-coordinate: $\pi(x, y) = (x, y) + (x, y) = (x^2 - y^2, 2xy) = (2x^2-1, 2xy)$

In Circle FFT, we use both the projection map $\pi_x$ and the squaring map $\pi$ to construct the sequence of domains. The sequence of domains for Circle FFT is shown as follows: $D_n \xrightarrow{\pi_x} S_n \xrightarrow{\pi} S_{n-1} \xrightarrow{\pi} \cdots \xrightarrow{\pi} 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 $\pi_x$ into two univariate polynomials $p_0(x)$ and $p_1(x)$ as follows: $p(x, y) = p_0(\pi_x(x, y)) + y \cdot p_1(\pi_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 = \pi_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}, \quad 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$.

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 $\mathbb{F}_{31}$, i.e, modulo $31$. We are given the evaluations $\vec{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 $\vec{v}$ over $D_3$. Step 1: . First, decompose the polynomial $p(x, y)$ using the map $\pi_x(x, y) = x$: $p(x, y) = p_0(\pi_x(x, y)) + y \cdot p_1(\pi_x(x, y)) = p_0(x) + y \cdot p_1(x)$ Given the evaluations $\vec{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 = \pi_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: $\vec{v_0} = [21, 6, 11, 10], \quad \vec{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 $\pi(x) = 2x^2 - 1$, and the process continues over successively smaller domains. We begin by decomposing $p_0(x)$ using the squaring map $\pi$ as follows: $p_0(x) = p_{00}(\pi(x)) + x \cdot p_{01}(\pi(x))$ We compute the evaluations of $p_{00}(x)$ and $p_{01}(x)$ over $S_{n-1} = \pi(S_n)$ using the following equations. $$p_{00}(\pi(x)) = \frac{p_0(x) + p_0(-x)}{2}, \quad p_{01}(\pi(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)$.

Step 2: Given the evaluations of $p_0$ and $p_1$ over $S_3 = [7, 13, 24, 18]$: $\vec{v_0} = [21, 6, 11, 10], \quad \vec{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 $\pi(x)$. For example, the decomposition of $p_0(x)$ is: $p_0(x) = p_{00}(\pi(x)) + x \cdot p_{01}(\pi(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 + \pi(x) + 28 \cdot x \cdot \pi(x)$ $p_1(x) = 11 + 26 \cdot x + 14 \cdot \pi(x) + 20 \cdot x \cdot \pi(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 + \pi(x) + 28 \cdot x \cdot \pi(x)$ $p_1(x) = 11 + 26 \cdot x + 14 \cdot \pi(x) + 20 \cdot x \cdot \pi(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 + \pi(x) + 28 \cdot x \cdot \pi(x) + y \cdot (11 + 26 \cdot x + 14 \cdot \pi(x) + 20 \cdot x \cdot \pi(x))$ $p(x, y) = 12 + 26 \cdot x + \pi(x) + 28 \cdot x \cdot \pi(x) + \quad \quad \quad \quad \quad \quad$ $\quad \quad \quad \quad \quad \quad 11 \cdot y + 26 \cdot x \cdot y + 14 \cdot y \cdot \pi(x) + 20 \cdot x \cdot y \cdot \pi(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.