Cyclotomic fast Fourier transform
# Cyclotomic fast Fourier transform

The cyclotomic fast Fourier transform is a type of fast Fourier transform algorithm over finite fields. 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 $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.

## 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 $\{f_{i}\}_{0}^{N-1}$ over a finite field $GF(p^m)$ is defined as

    $F_{j}=\sum _{i=0}^{N-1}f_{i}\alpha ^{ij},0\leq j\leq N-1,$

where $\alpha$ is the N-th primitive root of 1 in $GF(p^m).$ If we define the polynomial representation of $\{f_{i}\}_{0}^{N-1}$ as

    $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 $F_{j}$ is simply $f(\alpha ^{j})$. That is, the discrete Fourier transform of a sequence converts it to a polynomial evaluation problem.

Written in matrix format,

    $\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 $O(N^{2})$ complexity. Fast Fourier transforms are just efficient algorithms evaluating the above matrix-vector product.

## Algorithm

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

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

$L(x)$ is called linearized because $L(x_{1}+x_{2})=L(x_{1})+L(x_{2})$, which comes from the fact that for elements $x_{1},x_{2}\in \mathrm {GF} (p^{m}),$$(x_{1}+x_{2})^{p}=x_{1}^{p}+x_{2}^{p}.$

Notice that $p$ is invertible modulo $N$ because $N$ must divide the order $p^{m}-1$ of the multiplicative group of the field $\mathrm {GF} (p^{m})$. So, the elements $\{0,1,2,\ldots ,N-1\}$ can be partitioned into $l+1$ cyclotomic cosets modulo $N$:

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

where $k_{i}=p^{m_{i}}k_{i}{\pmod {N}}$. Therefore, the input to the Fourier transform can be rewritten as

    $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 $F_{j}$ is given by

    $F_{j}=f(\alpha ^{j})=\sum _{i=0}^{l}L_{i}(\alpha ^{jk_{i}})$.

Expanding $\alpha ^{jk_{i}}\in \mathrm {GF} (p^{m_{i}})$ with the proper basis $\{\beta _{i,0},\beta _{i,1},\ldots ,\beta _{i,m_{i}-1}\}$, we have $\alpha ^{jk_{i}}=\sum _{s=0}^{m_{i}-1}a_{ijs}\beta _{i,s}$ where $a_{ijs}\in \mathrm {GF} (p)$, and by the property of the linearized polynomial $L_{i}(x)$, we have

    $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 $\mathbf {F} =\mathbf {AL\Pi f}$, where $\mathbf {A}$ is an $N\times N$ matrix over GF(p) that contains the elements $a_{ijs}$, $\mathbf {L}$ is a block diagonal matrix, and $\mathbf {\Pi }$ is a permutation matrix regrouping the elements in $\mathbf {f}$ according to the cyclotomic coset index.

Note that if the normal basis $\{\gamma _{i}^{p^{0}},\gamma _{i}^{p^{1}},\cdots ,\gamma _{i}^{p^{m_{i}-1}}\}$ is used to expand the field elements of $\mathrm {GF} (p^{m_{i}})$, the i-th block of $\mathbf {L}$ is given by:

    $\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.

## Complexity

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

