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)=a0+a1x+...+an1xn1p(x) = a_0 + a_1 x + ... + a_{n-1} x^{n-1} and let n=2kn = 2k.

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

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

It holds:

(p(1)p(ω)...p(ωn1))=(111...11ωω2...ωn1...............1ωn1ω2(n1)...ω(n1)(n1))(a0a1...an1)\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)p(x) as:

peven(x)=a0+a2x+...+an2xk1p_{even}(x) = a_0 + a_2 x + ... + a_{n-2} x^{k-1}

podd(x)=a1+a3x+...+an1xk1p_{odd}(x) = a_1 + a_3 x + ... + a_{n-1} x^{k-1}

p(x)=peven(x2)+xpodd(x2)p(x) = p_{even}(x^2) + x p_{odd}(x^2)

Note that:

(p(1)p(ω)...p(ωk1)p(ωk)p(ωk+1)...p(ωn1))=(peven(12)+1podd(12)peven(ω2)+ωpodd(ω2)...peven(ω2(k1))+ωk1podd(ω2(k1))peven(ω2k)+ωkpodd(ω2k)peven(ω2(k+1))+ωk+1podd(ω2(k+1))...peven(ω2(n1))+ωn1podd(ω2(n1)))\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 nn:

(p(1)p(ω)...p(ωk1)p(ωk)p(ωk+1)...p(ωn1))=(peven(12)+1podd(12)peven(ω2)+ωpodd(ω2)...peven(ω2(k1))+ωk1podd(ω2(k1))peven(1)+ωkpodd(1)peven(ω2)+ωk+1podd(ω2)...peven(ω2(k1))+ω2k1podd(ω2(k1)))\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 pevenp_{even} and poddp_{odd} as the second half. We only need to compute pevenp_{even} and poddp_{odd} for kk values.

This way, instead of evaluating the polynomial pp in 2k2k values, we evaluate two polynomials in kk 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 logn\log n steps.

And nlognn \cdot \log n is much less than n2n^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)=a7x7+a6x6+a5x5+a4x4+a3x3+a2x2+a1x+a0p(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 nn coefficients (its degree is n1n-1), because that’s what we need to describe values in the domain of size nn: {1,w,w2,...,wn1}.\{1, w, w^2, ..., w^{n-1}\}. In our case n=8n=8.

Let’s split pp into pevenandpoddp_{even} \text{ and } p_{odd}:

peven(x)=a6x3+a4x2+a2x+a0p_{even}(x) = a_6x^3 + a_4x^2 + a_2x + a_0

podd(x)=a7x3+a5x2+a3x+a1p_{odd}(x) = a_7x^3 + a_5x^2 + a_3x + a_1

p(x)=peven(x2)+xpodd(x2)p(x) = p_{even}(x^2) + x \cdot p_{odd}(x^2)

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

But we split further:

peveneven(x)=a4x+a0p_{eveneven}(x) = a_4x + a_0

pevenodd(x)=a6x+a2p_{evenodd}(x) = a_6x + a_2

peven(x2)=peveneven(x4)+x2pevenodd(x4)p_{even}(x^2) = p_{eveneven}(x^4) + x^2 \cdot p_{evenodd}(x^4)

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

And to compute p(x)p(x), we again just need to evaluate the function peven(x2)+xpodd(x2).p_{even}(x^2) + x \cdot p_{odd}(x^2). Note that we already have peven(x2)p_{even}(x^2) and podd(x2).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=(111...11ωω2...ωn1...............1ωn1ω2(n1)...ω(n1)(n1))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 FF for FFT actually, but it becomes helpful to see how to construct the inverse function.

Some quick calculation shows that the matrix

G=1n(111...11ω1ω2...ω(n1)...............1ω(n1)ω2(n1)...ω(n1)(n1))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 GG. So we can proceed similarly as above.

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