Circle STARKs: Circle FFT

This is Part III of a four-part series on Circle STARKs. The original appeared on the ZK Security blog.

Recap#

In Part I of this series, we explored efficient field arithmetic over Mersenne primes and a key challenge in instantiating the STARK protocol using them, specifically the absence of a smooth multiplicative subgroup.

In Part II, we examined how this issue can be addressed using the circle group. We introduced the structure of its cosets, i.e., twin-coset and standard position coset, which can be used to instantiate STARK protocols. We also discussed the space of polynomials defined over the circle group, and the construction of vanishing polynomial over these cosets.

Introduction#

There are two common ways to represent a univariate polynomial $f(X)$ of degree less than $N$:

1. Monomial Basis:

   $$f(X) = \sum_{j = 0}^{N - 1} c_j \cdot b_j(X)$$

   where $c_j$ are the coefficients and $b_j(X)$ are the monomial basis polynomials such that:

   $$b(X) = [1, X, X^2, \cdots, X^{N - 1}]$$

2. Lagrange Basis:

   $$f(X) = \sum_{j = 0}^{N - 1} v_j \cdot \lambda_j(X)$$

   where $v_j = f(\omega_j)$ is the evaluation of $f$ at point $\omega_j$ and $\lambda_j(X)$ is the Lagrange basis polynomial satisfying:

   • $\lambda_j(\omega_j) = 1$
   • $\lambda_j(\omega_i) = 0$ for all $i \ne j$

   The Lagrange basis polynomial is defined as:

   $$\lambda_j(X) = \prod_{\substack{0 \le i < N, \\ i \ne j}} \frac{X - \omega_i}{\omega_j - \omega_i}$$

In general, the FFT algorithm performs a transformation from one polynomial basis to another. For example, it can convert a polynomial from the monomial basis to the Lagrange basis, and vice versa. This change of basis is useful because certain operations, such as polynomial multiplication, become significantly more efficient in the Lagrange basis.

In STARK protocols, the computation trace is usually represented as a set of polynomial evaluations at the points in a multiplicative subgroup, which means the coefficients of the polynomial can be computed using a Cooley-Tukey style FFT.

In Circle STARKs, the multiplicative subgroup is replaced by a twin-coset: not a group, but rather a union of cosets of the unit circle. The elements follow the “addition of angles” or “complex multiplication” group law, and each element corresponds to a point $(x, y)$ on the circle. In the Circle FFT, given evaluations over a twin-coset, we interpolate a bivariate polynomial.

In this post, we develop an intuitive understanding of the Circle FFT by drawing parallels to the classical Cooley–Tukey FFT. Throughout the post, we use FFT to refer to interpolation, moving from the Lagrange basis to a monomial basis, and inverse FFT to refer to evaluation, moving from the monomial basis to the Lagrange basis.

Cooley-Tukey FFT#

In this section, we present a generic version of the Cooley–Tukey FFT algorithm. This general version will eventually help us understand the Circle FFT as a concrete instance of the generic construct. To build intuition, we’ll walk through a concrete example alongside the general description.

Problem Statement#

Given the polynomial evaluations over a domain $D_n$ of size $N = k^n$, our goal is to compute the polynomial’s coefficients with respect to some basis polynomials such that:

Here, $c_j$ are the coefficients we wish to compute and $b_j^{(n)}(X)$ denotes the basis polynomials.

Algorithm#

The Cooley–Tukey FFT algorithm follows a classic divide-and-conquer strategy. It proceeds as follows:

• Step 1: Decompose the polynomial $f(X)$ into a collection of lower-degree polynomials and find the evaluations of those lower-degree polynomials over a smaller domain.

• Step 2: Recursively apply FFT to each of the lower-degree polynomials to compute their coefficients.

• Step 3: Combine the coefficients of the lower-degree polynomials to compute the coefficients of the original polynomial $f(X)$.

At each recursive layer $l$ of the FFT algorithm, we decompose a polynomial $f(X)$ using the following key components:

• $\psi$: The polynomial map $\psi$ is a $k$-to-$1$ map, meaning that it maps $k$ elements of the domain $D_l$ to a single element in the smaller domain $D_{l-1}$.

  $$\psi: D_l \rightarrow D_{l-1}$$

  Since $|D_l| = k^l$, this implies $|D_{l-1}| = k^{l-1}$. A common example of such a map over a multiplicative subgroup is: $\psi(X) = X^k$.

• $t_0, t_1, \cdots, t_{k-1}$: The twiddles $t_0, t_1, \cdots, t_{k-1}$ are polynomials which map elements from the domain $D_l$ to the field $\mathbb{F}$, i.e., for $j \in \{0, 1, \cdots, k-1\}$

  $$t_j: D_l \rightarrow \mathbb{F}$$

  For a fixed $z \in D_{l-1}$, the fiber of $z$ under the map $\psi$ consists of all elements in $D_l$ that map to $z$. Let $S_z$ denote the fiber of an element $z \in D_{l-1}$ under the polynomial map $\psi$, that is:

  $$S_z = \psi^{-1}(z) = \{x \in D_l \mid \psi(x) = z\} = \{x_0,x_1, \cdots, x_{k-1}\}$$

  Since $\psi$ is a $k$-to-1 map, the fiber of each $z \in D_{l-1}$ is a subset of $D_l$, containing exactly $k$ elements. The restriction of $f(X)$ to the fiber $S_z$ is a polynomial $f(X)|_{S_z}$ with the following evaluations:

  $$f|_{S_z} = \{f(x_0),f(x_1), \cdots, f(x_{k-1})\}$$

  Thus $f(X)|_{S_z}$ is of degree less than $k$ and it lies in a polynomial space of dimension at most $k$. The twiddles must be chosen to span this $k$-dimensional polynomial space. This ensures that the polynomial $f(X)$ can be correctly expressed in terms of the twiddles. Thus the twiddles must form a basis of a $k$-dimensional polynomial space over a fiber. This requirement will become clearer as we work through the algorithm and the example.

The polynomial map $\psi$ and the twiddles $t_0, t_1, \cdots, t_{k-1}$ need not be same across all layers of the FFT algorithm. We now describe each of the steps in more detail.

Step 1: Decompose $f(X)$ and compute evaluations of sub-polynomials#

We decompose the polynomial $f(X)$ into smaller-degree polynomials $f_i(X)$ as follows:

Each $f_i(X)$ is a polynomial of lower degree than $f(X)$, specifically:

Next, we compute the evaluations of each $f_i(X)$ over $D_{n-1} = \psi(D_n)$ using the evaluations of $f(X)$ over $D_n$. Let $S_z$ denote the fiber of an element $z \in D_{n-1}$ under the polynomial map $\psi$. Since it has $k$ elements, let us denote it as follows:

Such that for all $x_i \in S_z$, $\psi(x_i)=z$. Substituting all $x_i \in S_z$, into the decomposition of $f(X)$ gives:

Note that since for all $x_i \in S_z$, $\psi(x_i)=z$, the $f_0(z), f_1(z), \ldots, f_{k-1}(z)$ are same across all the equations. This gives us a system of $k$ linear equations in the $k$ unknowns $f_0(z), f_1(z), \ldots, f_{k-1}(z)$. Solving this system allows us to determine the values of each $f_i(z)$. Since the system has a constant size $k$, so solving it using e.g. Gaussian elimination takes constant time.

By repeating this process for every $z \in D_{n-1}$, we obtain the full set of evaluations for all the polynomials $f_0, f_1, \ldots, f_{k-1}$ over the domain $D_{n-1}$, enabling us to recursively apply the FFT on each of them.

Step 2: Recursively apply FFT to sub-polynomials#

Now that we have the evaluations of each $f_i$ over the domain $D_{n-1}$, the next step is to compute their coefficients. This brings us back to the same problem as before, but for polynomials of smaller degree i.e. $\deg(f_i) < k^{n-1}$ and over a smaller domain $D_{n-1}$. We apply the same divide-and-conquer strategy.

The polynomial map $\psi$ induces a sequence of progressively smaller domains:

where each domain $D_j = \psi(D_{j+1})$ for $j = 0, 1, \ldots, n - 1$, and $|D_j| = k^j$. We have assumed that the polynomial map $\psi$ is the same across all layers, but in general, each layer may use a different map.

We recursively apply FFT over smaller domains until we reach the base case: a constant polynomial, where all the evaluations over the domain are equal. In this case, interpolating gives a constant polynomial whose sole coefficient is the evaluation value itself. Thus, we simply return this evaluation value as the coefficient.

Once the base case is resolved, we work backward through the recursive calls, using the decomposition equation at each level to reconstruct the coefficients of the original $f_i$.

Step 3: Combine the coefficients#

Finally, we combine the coefficients of the polynomials $f_0, f_1, \ldots, f_{k-1}$ to compute the coefficients of the original polynomial $f(X)$ using the decomposition equation.

Sequence of Domains for Circle FFT#

In Circle FFT, the map $\psi$ is instantiated by a 2-to-1 map tailored for the circle group setting. The domain used here is the twin-coset:

where $G_{n-1}$ is a subgroup of the circle group of size $2^{n-1}$ and $Q \neq Q^{-1}$, hence the cosets $Q \cdot G_{n-1}$ and $Q^{-1} \cdot G_{n-1}$ are disjoint and $|D_n| = 2^n$. This section describes two specific 2-to-1 maps on $D_n$ that are central to the Circle FFT construction:

1. Quotient map $\phi_J$: This map sends each point in $D_n$ to an element of the quotient set $D_n / J$: $$\phi_J: D_n \rightarrow D_n / J.$$ where $J$ denotes the inversion map $J(P) = -P$. Quotient map $\phi_J$ maps two distinct points $P$ and $J(P)$ from $D_n$ to a single element in $D_n / J$. This can be interpreted as saying that two points are considered equivalent if they differ only by sign. This is same as the projection map $\pi_x$ that projects the points $(x,y)$ and $(x,-y)$ to the same $x$-coordinate: $$\pi_x((x, y)) = x$$ Since $\phi_J$ or $\pi_x$ map two points from $D_n$ to a single element in $D_n / J$, the size of $D_n / J$ will be half the size of $D_n$.
2. Squaring map $\pi$: This map was introduced in Part II. When applied to the twin-coset $D_n$ of size $2^n$, it produces a smaller twin-coset $D_{n-1}$ of size $2^{n-1}$. This is also a 2-to-1 map which halves the size of $D_{n}$. By recursively applying this map: $$D_j = \pi(D_{j+1}) \quad \text{for } j = n-1 \text{ down to } 1,$$ we generate a sequence of twin-cosets, until we reach $D_1$ (twin-coset of size 2). Each $D_j$ in this sequence corresponds to a twin-coset associated with the subgroup $G_{j-1}$, and has size $|D_j| = 2^j$.

In Circle FFT, we use both the projection map $\pi_x$ and the squaring map $\pi$ to construct the sequence of domains. A key property to note is that the maps $\pi_x$ and $\pi$ commute, meaning:

Intuitively, this means that projecting to the $x$-coordinate and then squaring is the same as squaring first and then projecting to the $x$-coordinate. This is further illustrated by the following two animations.

The following animation shows applying the projection map $\pi_x$ to a twin-coset $D_3$ followed by the squaring map $\pi$ to get the quotient domain $S_2$.

The following animation shows applying the squaring map $\pi$ to a twin-coset $D_3$ followed by the projection map $\pi_x$ to get the quotient domain $S_2$.

By applying these maps at each layer, we obtain a sequence of domains for the Circle FFT, similar to the domains in the Cooley–Tukey FFT. The key difference lies in the choice of maps: in Cooley–Tukey, we applied a single polynomial map $\psi$ at each layer, whereas in Circle FFT, we apply $\psi = \pi_x$ in the first layer, followed by $\psi = \pi$ in subsequent layers.

Let $S_j = D_j / J$ denote the quotient of the twin-coset $D_j$ under the involution $J$. The map $\pi: S_j \rightarrow S_{j-1}$ is a 2-to-1 map, defined by:

This is obtained using the doubling map and the equality $y^2 = 1 - x^2$ to compute the $x$-coordinate:

The sequence of domains for Circle FFT is shown as follows:

We will now describe the Circle FFT algorithm over these sequence of domains.

Circle FFT#

Let $\mathbb{F}_p$ be a Mersenne prime field with $p = 2^{31} - 1$, and let $\mathbb{F}$ denote its algebraic closure.

Given evaluations over a twin-coset $D_n$, the Circle FFT interpolates a bivariate polynomial $f(X, Y)$ from the space $L_N(\mathbb{F})$. This space consists of bivariate polynomials over the circle curve $X^2 + Y^2 = 1$ of total degree $\leq N/2$:

Let’s walk through the Circle FFT algorithm. As in the Cooley–Tukey case, we’ll also work through a concrete example alongside the algorithm.

Algorithm#

Circle FFT follows the same divide-and-conquer strategy as the Cooley–Tukey FFT. The bivariate polynomial is first decomposed into lower-degree polynomials. We then recursively compute the coefficients of these smaller polynomials and finally combine them to compute the coefficients of the original polynomial.

The key difference lies in how the polynomial is split. In Cooley–Tukey FFT, a single map $\psi$ is applied recursively at every level. In contrast, Circle FFT uses two distinct 2-to-1 maps: the projection map $\pi_x$ is applied in the first step, and the squaring map $\pi$ is used in all subsequent recursive steps.

Effectively, Circle FFT mirrors Cooley–Tukey, but with a split using $\psi = \pi_x$ initially, followed by recursive splits using $\psi = \pi$. We now walk through each step in detail.

Step 1: Decompose $f(X, Y)$ and compute evaluations of sub-polynomials#

In the first step, we split the bivariate polynomial $f(X, Y)$ over $D_n$ using the projection map $\pi_x$ and twiddles $t_0(X, Y) = 1$, $t_1(X, Y) = Y$:

Given the evaluations of $f(X, Y)$ over $D_n$, our goal is to compute the evaluations of $f_0(X)$ and $f_1(X)$ over the domain $S_n = \pi_x(D_n)$.

For any pair of points $P = (x, y)$ and its conjugate $J(P) = (x, -y)$ in $D_n$ that project to the same $x \in S_n$, substituting into the decomposition gives:

This system yields two equations in the unknowns $f_0(x)$ and $f_1(x)$, which can be solved directly.

Applying this to all pairs of points in $D_n$, gives the evaluations of $f_0(X)$ and $f_1(X)$ over the domain $S_n$.

Step 2: Recursively apply FFT to sub-polynomials#

Given the evaluations of $f_0(X)$ and $f_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 $f_0(X)$ using the squaring map $\pi$ and twiddles $t_{00}(X) = 1$, $t_{01}(X) = X$:

We compute the evaluations of $f_{00}(X)$ and $f_{01}(X)$ over $S_{n-1} = \pi(S_n)$ by solving a system of linear equations, as before.

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 $f_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 $f_1(X)$.

Step 3: Combine the coefficients#

Finally, we combine the coefficients of $f_0(X)$ and $f_1(X)$ to compute the coefficients of the original bivariate polynomial $f(X, Y)$, using the decomposition:

Implementation#

Here’s the Circle FFT implementation written in SageMath.

def cfft(evals, D):
    n = len(evals)
    assert (n & (n - 1)) == 0
    
    # Step 1: Split using map \pi_x
    S = pi_x(D)
    eval0 = []
    eval1 = []
    for e0, e1, P in zip(evals[:n//2], evals[n//2:], D):
        eval0.append((e0 + e1) / 2)
        eval1.append((e0 - e1) / (2 * get_y(P)))
    
    # Step 2: Recursively apply FFT using map \pi
    coeffs0 = cfft_pi(eval0, S)
    coeffs1 = cfft_pi(eval1, S)

    # Step 3: Recombine to compute coefficients
    coeffs = to_natural_order(coeffs0 + coeffs1)

    return coeffs

def cfft_pi(evals, S):
    n = len(evals)
    assert (n & (n - 1)) == 0

    if n == 1:
        return [evals[0]]
    
    S_next = pi(S)
    eval0 = []
    eval1 = []
    for e0, e1, x in zip(evals[:n//2], evals[n//2:], S):
        eval0.append((e0 + e1) / 2)
        eval1.append((e0 - e1) / (2 * x))
    
    coeffs0 = cfft_pi(eval0, S_next)
    coeffs1 = cfft_pi(eval1, S_next)
        
    return coeffs0 + coeffs1

Basis for FFT#

The FFT algorithm computes the coefficients with respect to some basis polynomials. For example, the Cooley-Tukey algorithm computes the coefficients $c_j$ with respect to basis polynomials $b_j^{(n)}(X)$ such that

The basis with respect to which the coefficients are computed depend on the choice of the polynomial map $\psi$ and the twiddles $t_0, t_1, \cdots t_{k-1}$ used at each layer of the FFT algorithm.

To see this in action, let us walk through a simple example. Consider the Cooley–Tukey FFT algorithm over a domain of size 4 to compute the coefficients of a polynomial of degree at most 3. First, we decompose the original polynomial $f(X)$:

Then, we decompose the two sub-polynomials $f_{0}(X)$ and $f_{1}(X)$:

Suppose $f_{00}(X)$, $f_{01}(X)$, $f_{10}(X)$ and $f_{11}(X)$ are constant polynomials. Then

Substitute the above values in decomposition of $f(X)$:

Thus the above algorithm computes the coefficients with respect to the basis: $t_0(X) \cdot t_0(\psi(X))$, $t_0(X) \cdot t_1(\psi(X))$, $t_1(X) \cdot t_0(\psi(X))$ and $t_1(X) \cdot t_1(\psi(X))$.

Thus depending on the choice of polynomial map $\psi$ and twiddles $t_0, t_1$, the FFT algorithm outputs the coefficients with respect to a particular basis.

• Choice 1: We use the polynomial map $\psi(X) = X^2$ and twiddles $t_0(X) = 1, t_1(X) = X$. This gives us coefficients with respect to the monomial basis as follows:

  $$b(X) = \{1, X, X^2, X^3\}$$

• Choice 2: We use the polynomial map $\psi(X) = X^2$ and twiddles $t_0(X) = X, t_1(X) = X+1$. This gives us coefficients with respect to the following basis:

  $$b(X) = \{X \cdot X^2,\; X^2 \cdot (X + 1),\; X \cdot (X^2 + 1),\; (X^2 + 1) \cdot (X + 1)\}$$

Similarly, circle FFT outputs coefficients $c_j$ with respect to some basis $b_j^{(n)}(X, Y)$ such that:

And the basis with respect to which the coefficients are computed depends on the choice of the polynomial map and the twiddles. In the concrete example of Circle FFT over twin-coset $D_3$, we saw that the algorithm computed the coefficient with respect to the following basis:

Using induction on $n$, we can show that the Circle FFT algorithm outputs coefficients with respect to the following basis [Theorem 2, Circle STARKs]:

where $0 \leq j \leq 2^n - 1$ and $(j_0, \ldots, j_{n-1}) \in \{0, 1\}^n$ is the binary representation of $j$, i.e.,