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 DFT into several circular convolutions, and then derives the DFT results from the circular convolution results. When applied to a DFT over ${\displaystyle 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]

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 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 GF(p^m). If we define the polynomial representation of ${\displaystyle \{f_{i}\}_{0}^{N-1}}$ 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 DFT has an ${\displaystyle O(N^{2})}$ complexity. Fast Fourier transforms are just efficient algorithms evaluating the above matrix-vector product.

First, we define a linearized polynomial over GF(p^m) as

${\displaystyle L(x)=\sum _{i}l_{i}x^{p^{i}},l_{i}\in \mathrm {GF} (p^{m}).}$

${\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}),}$${\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 \{0\},}$

${\displaystyle \{k_{1},pk_{1},p^{2}k_{1},\ldots ,p^{m_{1}-1}k_{1}\},}$

${\displaystyle \ldots ,}$

${\displaystyle \{k_{l},pk_{l},p^{2}k_{l},\ldots ,p^{m_{l}-1}k_{l}\},}$

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, and 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}\}}$, we have ${\displaystyle \alpha ^{jk_{i}}=\sum _{s=0}^{m_{i}-1}a_{ijs}\beta _{i,s}}$ where ${\displaystyle a_{ijs}\in \mathrm {GF} (p)}$, and by the property of the linearized polynomial ${\displaystyle L_{i}(x)}$, we have

${\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 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 we successfully reduce the discrete Fourier transform into short convolutions.

When applied to a characteristic-2 field 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]