FFT in finite fields

The finite field fast Fourier transform (opens new window) is incredibly useful in cryptography, in particular in modern ZKP (opens new window) schemes like PLONK.

Let’s have the polynomial $p(x) = a_0 + a_1 x + ... + a_{n-1} x^{n-1}$ and let $n = 2k$.

Also, let’s have the field $F_p$ and in it, we have the multiplicative subgroup $H = \{\omega^1, \omega^2,..., \omega^n = 1\}.$

We would like to have an efficient evaluation of the polynomial $p$ at all points of $H$ (this is for example needed in PLONK). The FFT enables us to do exactly this.

It holds:

$\begin{pmatrix} p(1) \\ p(\omega) \\ ... \\ p(\omega^{n-1}) \end{pmatrix} = \begin{pmatrix}1 & 1 & 1 &... & 1 \\1 & \omega & \omega^2 & ... & \omega^{n-1} \\... & ... & ... & ...& ... \\1 & \omega^{n-1} & \omega^{2(n-1)} & ... & \omega^{(n-1)(n-1)} \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ ... \\ a_{n-1} \end{pmatrix}$

The FFT trick is to express $p(x)$ as:

$p_{even}(x) = a_0 + a_2 x + ... + a_{n-2} x^{k-1}$

$p_{odd}(x) = a_1 + a_3 x + ... + a_{n-1} x^{k-1}$

$p(x) = p_{even}(x^2) + x p_{odd}(x^2)$

Note that:

$\begin{pmatrix} p(1) \\ p(\omega) \\ ... \\ p(\omega^{k-1}) \\ p(\omega^k) \\ p(\omega^{k+1}) \\ ... \\ p(\omega^{n-1}) \end{pmatrix} = \begin{pmatrix} p_{even}(1^2) + 1 \cdot p_{odd}(1^2) \\ p_{even}(\omega^2) + \omega \cdot p_{odd}(\omega^2) \\ ... \\ p_{even}(\omega^{2(k-1)}) + \omega^{k-1} \cdot p_{odd}(\omega^{2(k-1)}) \\ p_{even}(\omega^{2k}) + \omega^k \cdot p_{odd}(\omega^{2k}) \\ p_{even}(\omega^{2(k+1)}) + \omega^{k+1} \cdot p_{odd}(\omega^{2(k+1)}) \\ ... \\ p_{even}(\omega^{2(n-1)}) + \omega^{n-1} \cdot p_{odd}(\omega^{2(n-1)}) \end{pmatrix}$

And we can simplify it to have exponents smaller than $n$:

$\begin{pmatrix} p(1) \\ p(\omega) \\ ... \\ p(\omega^{k-1}) \\ p(\omega^k) \\ p(\omega^{k+1}) \\ ... \\ p(\omega^{n-1}) \end{pmatrix} = \begin{pmatrix} p_{even}(1^2) + 1 \cdot p_{odd}(1^2) \\ p_{even}(\omega^2) + \omega \cdot p_{odd}(\omega^2) \\ ... \\ p_{even}(\omega^{2(k-1)}) + \omega^{k-1} \cdot p_{odd}(\omega^{2(k-1)}) \\ p_{even}(1) + \omega^k \cdot p_{odd}(1) \\ p_{even}(\omega^{2}) + \omega^{k+1} \cdot p_{odd}(\omega^{2}) \\ ... \\ p_{even}(\omega^{2(k-1)}) + \omega^{2k-1} \cdot p_{odd}(\omega^{2(k-1)}) \end{pmatrix}$

The first half of the vector uses the same arguments for $p_{even}$ and $p_{odd}$ as the second half. We only need to compute $p_{even}$ and $p_{odd}$ for $k$ values.

This way, instead of evaluating the polynomial $p$ in $2k$ values, we evaluate two polynomials in $k$ values. Hey, but isn’t this the same number of operations?

Yeah, it’s even more operations, because we need to do the multiplications with the powers of $\omega$ and some additions too. However, if we recursively split the polynomials in this way, we will end up with $\log n$ steps.

And $n \cdot \log n$ is much less than $n^2$, which is the number of operations needed if we evaluate the polynomial the naive way (like in the first matrix expression at the top of the page).

Let’s observe the example:

$p(x) = a_7x^7 + a_6x^6 + a_5x^5 + a_4x^4 + a_3x^3 + a_2x^2 + a_1x + a_0$

The polynomial has $n$ coefficients (its degree is $n-1$), because that’s what we need to describe values in the domain of size $n$: $\{1, w, w^2, ..., w^{n-1}\}.$ In our case $n=8$.

Let’s split $p$ into $p_{even} \text{ and } p_{odd}$:

$p_{even}(x) = a_6x^3 + a_4x^2 + a_2x + a_0$

$p_{odd}(x) = a_7x^3 + a_5x^2 + a_3x + a_1$

$p(x) = p_{even}(x^2) + x \cdot p_{odd}(x^2)$

We would need to evaluate both polynomials at $\{1, w^2, w^4, w^6\}.$

But we split further:

$p_{eveneven}(x) = a_4x + a_0$

$p_{evenodd}(x) = a_6x + a_2$

$p_{even}(x^2) = p_{eveneven}(x^4) + x^2 \cdot p_{evenodd}(x^4)$

To compute $p_{even}(x)$ we just need to have $x^2$ (which is precomputed because $x$ is some power of $\omega$) and evaluate two functions, $p_{eveneven}$ and $p_{evenodd}$. But we need to evaluate these two functions only at $\{1, w^4\}$. We do the same for $p_{odd}$.

And to compute $p(x)$, we again just need to evaluate the function $p_{even}(x^2) + x \cdot p_{odd}(x^2).$ Note that we already have $p_{even}(x^2)$ and $p_{odd}(x^2).$

Now to the inverse operation – the inverse transformation in our case is a polynomial interpolation. That is, we have a set of evaluation points and need to find the coefficients of the polynomial.

Let’s observe the matrix we had above:

$F = \begin{pmatrix}1 & 1 & 1 &... & 1 \\1 & \omega & \omega^2 & ... & \omega^{n-1} \\... & ... & ... & ...& ... \\1 & \omega^{n-1} & \omega^{2(n-1)} & ... & \omega^{(n-1)(n-1)} \end{pmatrix}$

We didn’t need the matrix $F$ for FFT actually, but it becomes helpful to see how to construct the inverse function.

Some quick calculation shows that the matrix

$G = \frac{1}{n} \begin{pmatrix}1 & 1 & 1 &... & 1 \\1 & \omega^{-1} & \omega^{-2} & ... & \omega^{-(n-1)} \\... & ... & ... & ...& ... \\1 & \omega^{-(n-1)} & \omega^{-2(n-1)} & ... & \omega^{-(n-1)(n-1)} \end{pmatrix}$

is the inverse of $G$. So we can proceed similarly as above.

Next time, I will write about how the FFT helps speed up schemes like PLONK.