Cyclotomic fast Fourier transform

The cyclotomic fast Fourier transform is a type of fast Fourier transform algorithm over finite fields.^[1] This algorithm first decomposes a discrete Fourier transform into several circular convolutions. This then derives the discrete Fourier transform results from the circular convolution results. When applied to a discrete Fourier transform over ${\displaystyle \mathrm {GF} (2^{m})}$, this algorithm has a very low multiplicative complexity. In practice, since there usually exist efficient algorithms for circular convolutions with specific lengths, this algorithm is very efficient.^[2]

Background

The discrete Fourier transform over finite fields finds widespread application in the decoding of error-correcting codes such as BCH codes and Reed–Solomon codes. Generalized from the complex field, a discrete Fourier transform of a sequence ${\displaystyle \{f_{i}\}_{0}^{N-1}}$ over a finite field ${\displaystyle \mathrm {GF} (p^{m})}$ is defined as $${\displaystyle F_{j}=\sum _{i=0}^{N-1}f_{i}\alpha ^{ij},0\leq j\leq N-1,}$$ where ${\displaystyle \alpha }$ is the N-th primitive root of 1 in ${\displaystyle \mathrm {GF} (p^{m})}$. If the polynomial representation of ${\displaystyle \{f_{i}\}_{0}^{N-1}}$ can defined as $${\displaystyle f(x)=f_{0}+f_{1}x+f_{2}x^{2}+\cdots +f_{N-1}x^{N-1}=\sum _{0}^{N-1}f_{i}x^{i},}$$ it is easy to see that ${\displaystyle F_{j}}$ is simply ${\displaystyle f(\alpha ^{j})}$. That is, the discrete Fourier transform of a sequence converts it to a polynomial evaluation problem.

Written in matrix format, $${\displaystyle \mathbf {F} =\left[{\begin{matrix}F_{0}\\F_{1}\\\vdots \\F_{N-1}\end{matrix}}\right]=\left[{\begin{matrix}\alpha ^{0}&\alpha ^{0}&\cdots &\alpha ^{0}\\\alpha ^{0}&\alpha ^{1}&\cdots &\alpha ^{N-1}\\\vdots &\vdots &\ddots &\vdots \\\alpha ^{0}&\alpha ^{N-1}&\cdots &\alpha ^{(N-1)(N-1)}\end{matrix}}\right]\left[{\begin{matrix}f_{0}\\f_{1}\\\vdots \\f_{N-1}\end{matrix}}\right]={\mathcal {F}}\mathbf {f} .}$$

Direct evaluation of the discrete Fourier transform has an ${\displaystyle O(N^{2})}$ complexity. Fast Fourier transforms are just efficient algorithms evaluating the above matrix-vector product.

Algorithm

First, define a linearized polynomial over ${\displaystyle \mathrm {GF} (p^{m})}$ as $${\displaystyle L(x)=\sum _{i}l_{i}x^{p^{i}},l_{i}\in \mathrm {GF} (p^{m}).}$$ Here ${\displaystyle L(x)}$ is called linearized, because ${\displaystyle L(x_{1}+x_{2})=L(x_{1})+L(x_{2})}$, which comes from the fact that for elements ${\displaystyle x_{1},x_{2}\in \mathrm {GF} (p^{m})}$ and ${\displaystyle (x_{1}+x_{2})^{p}=x_{1}^{p}+x_{2}^{p}}$.

Notice that ${\displaystyle p}$ is invertible modulo ${\displaystyle N}$ because ${\displaystyle N}$ must divide the order ${\displaystyle p^{m}-1}$ of the multiplicative group of the field ${\displaystyle \mathrm {GF} (p^{m})}$. So, the elements ${\displaystyle \{0,1,2,\ldots ,N-1\}}$ can be partitioned into ${\displaystyle l+1}$ cyclotomic cosets modulo ${\displaystyle N}$: $${\displaystyle {\begin{aligned}&\{0\},\\&\{k_{1},pk_{1},p^{2}k_{1},\ldots ,p^{m_{1}-1}k_{1}\},\\&\ldots ,\\&\{k_{l},pk_{l},p^{2}k_{l},\ldots ,p^{m_{l}-1}k_{l}\},\end{aligned}}}$$ where ${\displaystyle k_{i}=p^{m_{i}}k_{i}{\pmod {N}}}$. Therefore, the input to the Fourier transform can be rewritten as $${\displaystyle f(x)=\sum _{i=0}^{l}L_{i}(x^{k_{i}}),\quad L_{i}(y)=\sum _{t=0}^{m_{i}-1}y^{p^{t}}f_{p^{t}k_{i}{\bmod {N}}}.}$$

In this way, the polynomial representation is decomposed into a sum of linear polynomials. Hence, ${\displaystyle F_{j}}$ is given by $${\displaystyle F_{j}=f(\alpha ^{j})=\sum _{i=0}^{l}L_{i}(\alpha ^{jk_{i}})}$$. Expanding ${\displaystyle \alpha ^{jk_{i}}\in \mathrm {GF} (p^{m_{i}})}$ with the proper basis ${\displaystyle \{\beta _{i,0},\beta _{i,1},\ldots ,\beta _{i,m_{i}-1}\}}$ yields ${\displaystyle \alpha ^{jk_{i}}=\sum _{s=0}^{m_{i}-1}a_{ijs}\beta _{i,s}}$, where ${\displaystyle a_{ijs}\in \mathrm {GF} (p)}$. By the property of the linearized polynomial ${\displaystyle L_{i}(x)}$, this yields $${\displaystyle F_{j}=\sum _{i=0}^{l}\sum _{s=0}^{m_{i}-1}a_{ijs}\left(\sum _{t=0}^{m_{i}-1}\beta _{i,s}^{p^{t}}f_{p^{t}k_{i}{\bmod {N}}}\right).}$$

This equation can be rewritten in matrix form as ${\displaystyle \mathbf {F} =\mathbf {AL\Pi f} }$, where ${\displaystyle \mathbf {A} }$ is an ${\displaystyle N\times N}$ matrix over ${\displaystyle \mathrm {GF} (p)}$ that contains the elements ${\displaystyle a_{ijs}}$, ${\displaystyle \mathbf {L} }$ is a block diagonal matrix, and ${\displaystyle \mathbf {\Pi } }$ is a permutation matrix regrouping the elements in ${\displaystyle \mathbf {f} }$ according to the cyclotomic coset index.

Note that if the normal basis ${\displaystyle \{\gamma _{i}^{p^{0}},\gamma _{i}^{p^{1}},\cdots ,\gamma _{i}^{p^{m_{i}-1}}\}}$ is used to expand the field elements of ${\displaystyle \mathrm {GF} (p^{m_{i}})}$, the i-th block of ${\displaystyle \mathbf {L} }$ is given by: $${\displaystyle \mathbf {L} _{i}={\begin{bmatrix}\gamma _{i}^{p^{0}}&\gamma _{i}^{p^{1}}&\cdots &\gamma _{i}^{p^{m_{i}-1}}\\\gamma _{i}^{p^{1}}&\gamma _{i}^{p^{2}}&\cdots &\gamma _{i}^{p^{0}}\\\vdots &\vdots &\ddots &\vdots \\\gamma _{i}^{p^{m_{i}-1}}&\gamma _{i}^{p^{0}}&\cdots &\gamma _{i}^{p^{m_{i}-2}}\\\end{bmatrix}},}$$ which is a circulant matrix. It is well known that a circulant matrix-vector product can be efficiently computed by convolutions. Hence, it is successful to reduce the discrete Fourier transform into short convolutions.

Complexity

When applied to a characteristic-2 field ${\displaystyle \mathrm {GF} (2^{m})}$, the matrix ${\displaystyle \mathbf {A} }$ is just a binary matrix. Only addition is used when calculating the matrix-vector product of ${\displaystyle \mathrm {A} }$ and ${\displaystyle \mathrm {L\Pi f} }$. It has been shown that the multiplicative complexity of the cyclotomic algorithm is given by ${\displaystyle O(n(\log _{2}n)^{\log _{2}{\frac {3}{2}}})}$, and the additive complexity is given by ${\displaystyle O(n^{2}/(\log _{2}n)^{\log _{2}{\frac {8}{3}}})}$.^[2]