Circle FFT and its implementation

Abstract: This paper primarily examines the work of Ulrich Haböck et al. titled "Circle STARKs," interpreting Circle FFT and its implementation aspects. The analysis focuses on three main areas: 1) the foundation for accelerating basic operations with CFFT; 2) the implementation of CFFT calculations; and 3) the coding implementation of CFFT.

1. Introduction

In traditional STARKs, Fast Fourier Transform (FFT) is used to efficiently perform interpolation and write adjacent row constraints. To utilize FFT for fundamental polynomial operations, the finite field

${\mathbb{F}}_p$ must possess a smooth-order root of unity. Besides the FFT itself, the choice of finite field

$p$ directly influences the complexity of arithmetic operations, thereby impacting the efficiency of STARKs. The most effective field for arithmetic operations is the Mersenne prime field, specifically

$p=2^e-1$. Notably,

$p=2^{31}-1$ allows for highly efficient arithmetic operations on 32-bit computers. For these reasons, traditional FFT is inadequate in meeting the efficiency demands of STARKs. In reference [1], the authors continue the ECFFT [2][3] approach, constructing Circle STARKs based on the Mersenne prime

$p=2^{31}-1$ on the circle defined by the equation

$x^2+y^2=1$. The core innovation lies in the introduction of Circle FFT (CFFT) within STARKs.

2. Circle FFT

2.1 Fundamentals of CFFT

The choice of field significantly affects the efficiency of FFT acceleration and the feasibility of using FFT for this purpose. A prime

$p$ that is friendly to CFFT possesses the property

$p\equiv 3mod4$. The set of

$p+1$ points on the curve

$C=C\left(\mathbb{F}p\right)$ forms a group, defined by the group operation

$(x0,y0)\cdot (x1,y1):=(x0x1-y0y1,x0y1+x1y0)$, where the circular group

$C\left(\mathbb{F}p\right)$ is a cyclic group. Additionally, we define the rotation operation

$T_P(x,y):=P\cdot (x,y)$, the squaring map

$\pi (x,y):=(x,y)\cdot (x,y)$, and the group inverse

$J(x,y):=(x,-y)$.

When using FFT acceleration, the length of the coefficients to be computed must be

$2^n(n\ge 1)$. For CFFT, we define a double coset of size

$N=2^n$ as

$D=Q\cdot Gn-1\cup Q^{-1}\cdot Gn-1$. Here,

$Gn-1$ is a subgroup of size

$2^{n-1}$ of

$C\left(\mathbb{F}p\right)$, and

$Q\in C\left(\mathbb{F}p\right)$. When the prime

$p\equiv 3mod4$, we have

$D=Q\cdot Gn=Q\cdot Gn-1\cup Q^{-1}\cdot Gn-1$. When

$D$ is a coset of the subgroup

$Gn$, it is the standard position coset of size

$N$. Furthermore, under the squaring map, the set

$\pi \left(D\right)$ of size

$N/2$ shares the same properties as

$D$, meaning it is also a double coset or standard coset. The standard position coset is illustrated in Figure 1, where the elements are evenly distributed along the circle, corresponding to the binary size of the set. For

$n\ge 1$，

${\pi }^{n-1}(D)=\{(x_D,y_D)\}$ contains only two elements.

Figure 1: Diagram of the three minimal standard position cosets in the affine plane over

${\mathbb{F}}_p$

In Circle STARKs, prior to CFFT acceleration, we first define the space of bivariate polynomials whose coefficients belong to the extension field

$F$:

$$\begin{array}{}\text{(1)}& {\mathcal{L}}_N\left(F\right)=\{p(x,y)\in F[x,y]/(x^2+y^2-1):\mathrm{deg}p\le \frac{N}{2}\}\end{array}$$

where

$N$ is an even number. The space

$\mathcal{L}N\left(F\right)$ has rotational invariance and good separability. These key properties are necessary for efficient encoding. By repeatedly substituting

$y^2=1-x^2$ into the polynomial

$p(x,y)$ derived from

$\mathcal{L}N\left(F\right)$, we obtain

$$\begin{array}{}\text{(2)}& p(x,y)=p_0\left(x\right)+y\cdot p_1\left(x\right)\end{array}$$，
where

$p0\left(x\right)\in F{\left[x\right]}^{\le N/2}$ and

$p1\left(x\right)\in F{\left[x\right]}^{\le N/2-1}$. Equation

$(2)$ is referred to as the canonical form of the polynomial derived from

${\mathcal{L}}_N\left(F\right)$. The canonical form indicates that the set of monomials

$$\begin{array}{}\text{(3)}& 1,x,\dots ,x^{\frac{N}{2}},y,y\cdot x,\dots ,y\cdot x^{\frac{N}{2}-1}\end{array}$$

spans the space

$\mathcal{L}N\left(F\right)$. Based on the dimension of

$\mathcal{L}N\left(F\right)$, the canonical form must be a basis, specifically a monomial basis. Additionally, for details on the cyclic code, refer to reference [1].

For the double coset

$D=Q\cdot Gn-1\cup Q^{-1}\cdot Gn-1$, with

$Q\in C\left(\mathbb{F}p\right)\mathrm{\setminus }Gn$, CFFT completes the interpolation of functions derived from

$F^D$ by calculating the coefficients corresponding to the basis of Circle FFT

${\mathcal{B}}_n$. The basis for Circle FFT consists of N-dimensional polynomial bases that depend solely on the size of the field. We define the

$n$-th order FFT basis as

$$\begin{array}{}\text{(4)}& b_j^{(n)}(x,y):=y^{j_0}\cdot v_1{(x)}^{j_1}\cdot \dots v_{n-1}^{j_{n-1}}(x),0\le j\le 2^n-1\end{array}$$

where

$0\le j\le 2^n-1$ and

$$j=j_0+j_1\cdot 2+\dots +j_{n-1}\cdot 2^{n-1}((j_0,\dots ,j_{n-1})\in {\{0,1\}}^n)$$;

$vk(x)$ is the vanishing polynomial of the standard position coset of size

$2^k$ [1]. The family formed by

$b_j^{(n)}$ constitutes the set

$\mathcal{B}n\in \mathcal{L}N(F)$.

According to the definitions above, the double coset

$D$ exhibits invariance under the wrapping operation

$J$, and each

$J$-orbit in

$D$ contains only two points. Consequently, the quotient mapping

$$\begin{array}{}\text{(5)}& {\varphi }_J:D\to D/J,P\mapsto \{P,J\left(P\right)\}\end{array}$$

is a two-to-one mapping. The double coset obtained through recursion is given by

$Dj=\pi \left(Dj+1\right)$, for

$j=n-1,\dots ,1$, corresponding to the descending chain of subgroups

$Gn-1\supset Gn-2\supset \dots \supset G_0$. Since the operations

$J$ and

$\pi$ are commutative, we can obtain the commutative diagram:

Figure 2: Commutative diagram of the J and
𝜋 mappings: each mapping is two-to-one, and the coset size is halved.

Alternatively, the quotient

$Sj=Dj/J$ can be considered as a subset of the x-axis, and

$\varphi J=\pi x$ serves as the projection onto the x-axis. Thus, the commutative diagram can be transformed into

Figure 3: Transformed commutative diagram.

In the diagram, the squaring homomorphism mapping

$\pi :Sj\to Sj-1$ is a two-to-one mapping given by

$x\mapsto 2\cdot x^2-1$.

2.2 Mathematical Description of CFFT

Given a double coset of size

$2^n$, FFT is a divide-and-conquer algorithm that iteratively simplifies the interpolation problem for

$f\in F^D$ along the projection chain:

$$\begin{array}{}\text{(6)}& D=D_n\stackrel{{\varphi }_J={\pi }_x}{\to }{S}_n\stackrel{\pi }{\to }{S}_{n-1}\stackrel{\pi }{\to }\dots \stackrel{\pi }{\to }{S}_1\end{array}$$
From the projection chain, it is evident that the entire chain is divided into two stages.

In the first stage, we use

$t0(x,y)=y$ as the "reference odd function" to decompose

$f\in F^D$ into "odd" and "even" parts according to the wrapping function

$J$, resulting in two unique functions

$f0$ and

$f_1\in F^{S_n}$ defined on the univariate domain

$S_n$, specifically:

$$\begin{array}{}\text{(7)}& f_0\left(x\right)=\frac{f(x,y)+f(x,-y)}{2}\end{array}$$

$$\begin{array}{}\text{(8)}& f_1\left(x\right)=\frac{f(x,y)-f(x,-y)}{2y}\end{array}$$
and

$$\begin{array}{}\text{(9)}& f(x,y)=f_0\left(x\right)+y\cdot f_1\left(x\right)\end{array}$$

If we factor out the coefficient

$\frac{1}{2}$ from equations (7) and (8) and consider

$y^{\prime }=y^{-1}$ as the rotation factor, this decomposition is quite similar to the butterfly operation in traditional FFT.

In the second stage, we select

$t1(x,y)=y$ as the "reference odd function." On the projection domain

$Sj-1=\pi \left(S_j\right)$, the iterative decomposition yields:

$$\begin{array}{}\text{(10)}& f_0\left(\pi \left(x\right)\right)=\frac{f\left(x\right)+f(-x)}{2}\end{array}$$

$$\begin{array}{}\text{(11)}& f_1\left(\pi \left(x\right)\right)=\frac{f\left(x\right)-f(-x)}{2\cdot x}\end{array}$$
and

$$\begin{array}{}\text{(12)}& f\left(x\right)=f_0\left(\pi \left(x\right)\right)+x\cdot f_1\left(\pi \left(x\right)\right)\end{array}$$

where

$\pi \left(x\right)=2\cdot x^2-1$. This decomposition process continues until the projection domain reduces to the single point

$S_1$.

The output of the CFFT algorithm is a constant

$ck=fk0,\dots ,kn-j\in F$, where

$0\le k\le 2^n-1$ and

$k0,\dots ,kn-1$ are the binary representation bits of

$k$. Thus, for a given basis space

$\mathcal{B}{j}^{(0)}\subseteq \mathcal{B}j$:

$$\begin{array}{}\text{(13)}& {\mathcal{B}}_j^{(0)}=\{b_{2k}^{\left(j\right)}:0\le k\le 2^{j-1}-1\},1\le j\le n\end{array}$$

We have:

$$\begin{array}{}\text{(14)}& f\left(x\right)=\sum _{k^{\prime }=0}^{2^j-1}{c}_{k^{\prime }}\cdot b_{2\cdot k^{\prime }}^{(j+1)}\left(x\right)\end{array}$$

When

$j=n$,

$c_{k^{\prime }}$ represents the output of the CFFT.

3. CFFT Implementation

In Algorithm [1], the authors provide the algorithm's pseudocode.

From the pseudocode, it is evident that the entire algorithm is divided into two stages: 1) operations at the first level; and 2) operations in the subsequent layers. The primary difference lies in the calculation of the rotation factors.

3.1 Calculation of Rotation Factors

As seen in the pseudocode, although the implementation of CFFT differs between the first layer and the subsequent layers, the inconsistency mainly arises from the calculation of rotation factors, while the butterfly operation remains consistent. In the RUST implementation, the calculation of rotation factors is

From lines 132 to 142 of the code, it can be observed that the calculation of rotation factors distinguishes between the first layer and the remaining layers. It is necessary to differentiate between the first step and the second step, which are completed in init_domain and working_domain, respectively. This represents a significant difference from the calculation of rotation factors in traditional FFT. The generation of the initial domain is accomplished by calling the cfft_domain function.

3.2 Butterfly Operator

In the implementation code, the butterfly operators for the forward and inverse transformations are as follows:

and

It can be seen that in the implementation, the butterfly operator of CFFT is essentially consistent with that of traditional FFT.

3.3 Basic Implementation of CFFT

After the rotation factors have been precomputed, CFFT can be implemented directly based on the pseudocode.

In the implementation, the butterfly operation is realized using direct code. Alternatively, it can be called directly using fn butterfly_cfft().

---

References:
[1] Ulrich Haböck and David Levit and Shahar Papini. Circle STARKs. In Cryptology ePrint Archive, Paper 2024/278, 2024. https://eprint.iacr.org/2024/278
[2] Eli Ben-Sasson, Dan Carmon, Swastik Kopparty, and David Levit. Elliptic Curve Fast Fourier Transform (ECFFT) Part I: Fast polynomial algorithms over all finite fields. In Electronic Colloquium on Computational Complexity, volume TR21-103, 2021. https: //eccc.weizmann.ac.il/report/2021/103/.
[3] Eli Ben-Sasson, Dan Carmon, Swastik Kopparty, and David Levit. Scalable and transparent proofs over all large fields, via elliptic curves (ECFFT part II). In IACR preprint archive, 2022. https://eprint.iacr.org/2022/1542.
[4] GitHub - Plonky3/Plonky3 at Circle-Fast-Fourier-Transform. https://github.com/Plonky3/Plonky3