Discrete Fast Fourier Transform

Available code files:

• fft2.sage: implementation using vandermonde matrices illustrating the theory section below.
• fft3.sage: simple example with $n=4$ showing 3 steps of the algorithm.
• fft4.sage: illustrates the full working algorithm.

Theory

$$f(x)g(x)\in {\mathbb{F}}_{<2n}[x]$$ $$fg=\sum _{i+j<2n-2}{a}_i{b}_j{x}^{i+j}$$ Complexity: $O(n^2)$

Suppose $\omega \in \mathbb{F}$ is an nth root of unity.

Recall: if $\mathbb{F}={\mathbb{F}}_{p^k}$ then $\exists N:{\mathbb{F}}_{p^N}$ contains all nth roots of unity.

$${\text{DFT}}_{\omega }:{\mathbb{F}}^n\to {\mathbb{F}}^n$$ $${\text{DFT}}_{\omega }(f)=(f({\omega }^0),f({\omega }^1),\dots ,f({\omega }^{n-1}))$$

$$V_{\omega }=\left(\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& {\omega }^1& {\omega }^2& \cdots & {\omega }^{n-1}\\ 1& {\omega }^2& {\omega }^4& \cdots & {\omega }^{2(n-1)}\\ ⋮& & & & \\ 1& {\omega }^{n-1}& {\omega }^{2(n-1)}& \cdots & {\omega }^{(n-1{)}^2}\end{array}\right)$$

$${\text{DFT}}_{\omega }(f)=V_{\omega }\cdot f^T$$ since vandermonde multiplication is simply evaluation of a polynomial.

Lemma: $V_{\omega }^{-1}=\frac{1}{n}{V}_{{\omega }^{-1}}$

Use $1+\omega +\cdots +{\omega }^{n-1}$ and compute $V_{\omega }{V}_{{\omega }^{-1}}$

Corollary: ${\text{DFT}}_{\omega }$ is invertible.

Definitions

1. Convolution $f\ast g=fg\mathrm{mod}(x^n-1)$
2. Pointwise product $$(a_0,\dots ,a_{n-1})\cdot (b_0,\dots ,b_{n-1})=(a_0{b}_0,\dots ,a_{n-1}{b}_{n-1})\in {\mathbb{F}}^n\to {\mathbb{F}}_{<n}[x]$$

Theorem: ${\text{DFT}}_{\omega }(f\ast g)={\text{DFT}}_{\omega }(f)\cdot {\text{DFT}}_{\omega }(g)$

$$fg=q^{\prime }(x^n-1)+f\ast g$$ $$\Rightarrow f\ast g=fg+q(x^n-1)$$ $$\mathrm{deg}fg\le 2n-2$$

$$\begin{array}{rl}(f\ast g)({\omega }^i)& =f({\omega }^i)g({\omega }^i)+q({\omega }^i)({\omega }^{in}-1)\\ & =f({\omega }^i)g({\omega }^i)\end{array}$$

Result

$$f,g\in {\mathbb{F}}_{<n/2}[x]$$ $$fg=f\ast g$$ $${\text{DFT}}_{\omega }(f\ast g)={\text{DFT}}_{\omega }(f)\cdot {\text{DFT}}_{\omega }(g)$$ $$fg=\frac{1}{n}{\text{DFT}}_{{\omega }^{-1}}({\text{DFT}}_{\omega }(f)\cdot {\text{DFT}}_{\omega }(g))$$

Finite Field Extension Containing Nth Roots of Unity

$${\mu }_N=⟨\omega ⟩,|{\mathbb{F}}_{p^N}^{\times }|=p^N-1$$ $$\text{ord}(\omega )=n|p^N-1$$ but ${\mathbb{F}}_{p^N}^{\times }$ is cyclic.

For all $d|p^N-1$, there exists $x\in {\mathbb{F}}_{p^N}^{\times }$ with ord$(x)=d$.

Finding $n|p^N-1$ is sufficient for $\omega \in {\mathbb{F}}_{p^N}$ $$n|p^N-1\iff \text{ord}(p)=(Z/nZ{)}^{\times }$$

FFT Algorithm Recursive Compute

We recurse to a depth of $\log n$. Since each recursion uses ${\omega }^i$, then in the final step ${\omega }^i=1$, and we simply return $f^T$.

We only need to prove a single step of the algorithm produces the desired result, and then the correctness is inductively proven.

$$\begin{array}{rl}f(X)& =a_0+a_{1}X+a_2{X}^2+\cdots +a_{n-1}{X}^{n-1}\\ & =g(X)+X^{n/2}h(X)\end{array}$$

Algorithm

Implementation of this algorithm is available in fft4.sage. Particularly the function called calc_dft().

function DFT($n=2^d,f(X)$)
if $n=1$ then
return $f(X)$
end
Write $f(X)$ as the sum of two polynomials with equal degree
$f(X)=g(X)+X^{n/2}h(X)$
Let $\mathbf{g},\mathbf{h}$ be the vector representations of $g(X),h(X)$

$\mathbf{r}=\mathbf{g}+\mathbf{h}$
$\mathbf{s}=(\mathbf{g}-\mathbf{h})\cdot ({\omega }^0,\dots ,{\omega }^{n/2-1})$
Let $r(X),s(X)$ be the polynomials represented by the vectors $\mathbf{r},\mathbf{s}$

Compute $(r({\omega }^0),\dots ,r({\omega }^{n/2}))={\text{DFT}}_{{\omega }^2}(n/2,r(X))$
Compute $(s({\omega }^0),\dots ,s({\omega }^{n/2}))={\text{DFT}}_{{\omega }^2}(n/2,s(X))$

return $(r({\omega }^0),s({\omega }^0),r({\omega }^2),s({\omega }^2),\dots ,r({\omega }^{n/2}),s({\omega }^{n/2}))$
end

Sage code:

def calc_dft(ω_powers, f):
    m = len(f)
    if m == 1:
        return f
    g, h = vector(f[:m/2]), vector(f[m/2:])

    r = g + h
    s = dot(g - h, ω_powers)

    ω_powers = vector(ω_i for ω_i in ω_powers[::2])
    rT = calc_dft(ω_powers, r)
    sT = calc_dft(ω_powers, s)

    return list(alternate(rT, sT))

Even Values

$$\begin{array}{rl}r(X)& =g(X)+h(X)\\ & \\ f({\omega }^{2i})& =g({\omega }^{2i})+({\omega }^{2i}{)}^{n/2}h({\omega }^{2i})\\ & =g({\omega }^{2i})+h({\omega }^{2i})\\ & =(g+h)({\omega }^{2i})\end{array}$$

So then we can now compute $DF{T}_{\omega }(f{)}_{k=2i}=DF{T}_{{\omega }^2}(r)$ for the even powers of $f({\omega }^{2i})$.

Odd Values

For odd values $k=2i+1$

$$\begin{array}{rl}s(X)& =(g(X)-h(X))\cdot ({\omega }^0,\dots ,{\omega }^{n/2-1})\\ & \\ f(X)& =a_0+a_{1}X+a_2{X}^2+\cdots +a_{n-1}{X}^{n-1}\\ & =g(X)+X^{n/2}h(X)\\ f({\omega }^{2i+1})& =g({\omega }^{2i+1})+({\omega }^{2i+1}{)}^{n/2}h({\omega }^{2i+1})\end{array}$$ But observe that for any $n$th root of unity ${\omega }^n=1$ and ${\omega }^{n/2}=-1$ $$({\omega }^{2i+1}{)}^{n/2}={\omega }^{in}{\omega }^{n/2}={\omega }^{n/2}=-1$$ $$\begin{array}{rl}\Rightarrow f({\omega }^{2i+1})& =g({\omega }^{2i+1})-h({\omega }^{2i+1})\\ & =(g-h)({\omega }^{2i+1})\end{array}$$

Let $\mathbf{s}=(\mathbf{g}-\mathbf{h})\cdot ({\omega }^0,\dots ,{\omega }^{n/2-1})$ be the representation for $s(X)$. Then we can see that $s({\omega }^{2i})=(g-h)({\omega }^{2i+1})$ as desired.

So then we can now compute $DF{T}_{\omega }(f{)}_{k=2i+1}=DF{T}_{{\omega }^2}(s)$ for the odd powers of $f({\omega }^{2i+1})$.

Example

Let $n=8$ $$\begin{array}{rl}f(X)& =(a_0+a_{1}X+a_2{X}^2+a_3{X}^3)+(a_4{X}^4+a_5{X}^5+a_6{X}^6+a_7{X}^7)\\ & =(a_0+a_{1}X+a_2{X}^2+a_3{X}^3)+X^4(a_4+a_{5}X+a_6{X}^2+a_7{X}^3)\\ & =g(X)+X^{n/2}h(X)\\ g(X)& =a_0+a_{1}X+a_2{X}^2+a_3{X}^3\\ h(X)& =a_4+a_{5}X+a_6{X}^2+a_7{X}^3\end{array}$$ Now vectorize $g(X),h(X)$ $$\begin{array}{rl}\mathbf{g}& =(a_0,a_1,a_2,a_3)\\ \mathbf{h}& =(a_4,a_5,a_6,a_7)\end{array}$$ Compute reduced polynomials in vector form $$\begin{array}{rl}\mathbf{r}& =\mathbf{g}+\mathbf{h}\\ & =(a_0+a_4,a_1+a_5,a_2+a_6,a_3+a_7)\\ \mathbf{s}& =(\mathbf{g}-\mathbf{h})\cdot (1,\omega ,{\omega }^2,{\omega }^3)\\ & =(a_0-a_4,a_1-a_5,a_2-a_6,a_3-a_7)\cdot (1,\omega ,{\omega }^2,{\omega }^3)\\ & =(a_0-a_4,\omega (a_1-a_5),{\omega }^2(a_2-a_6),{\omega }^3(a_3-a_7))\end{array}$$ Convert them to polynomials from the vectors. We also expand them out below for completeness. $$\begin{array}{rl}r(X)& =r_0+r_{1}X+r_2{X}^2+r_3{X}^3\\ & =(a_0+a_4)+(a_1+a_5)X+(a_2+a_6)X^2+(a_3+a_7)X^3\\ s(X)& =s_0+s_{1}X+s_2{X}^2+s_3{X}^3\\ & =(a_0-a_4)+\omega (a_1-a_5)X+{\omega }^2(a_2-a_6)X^2+{\omega }^3(a_3-a_7)X^3\end{array}$$ Compute $${\text{DFT}}_{{\omega }^2}(4,r(X)),{\text{DFT}}_{{\omega }^2}(4,s(X))$$ The values returned will be $$(r(1),s(1),r({\omega }^2),s({\omega }^2),r({\omega }^4),s({\omega }^4),r({\omega }^6),s({\omega }^6))=(f(1),f(\omega ),f({\omega }^2),f({\omega }^3),f({\omega }^4),f({\omega }^5),f({\omega }^6),f({\omega }^7))$$ Which is the output we return.

Comparing Evaluations for $f(X)$ and $r(X),s(X)$

We can see the evaluations are correct by substituting in ${\omega }^i$.

We expect that $s(X)$ on the domain $(1,{\omega }^2,{\omega }^4,{\omega }^6)$ produces the values $(f(1),f({\omega }^2),f({\omega }^4),f({\omega }^6))$, while $r(X)$ on the same domain produces $(f(\omega ),f({\omega }^3),f({\omega }^5),f({\omega }^7))$.

Even Values

Let $k=2i$, be an even number. Then note that $k$ is a multiple of 2, so $4k$ is a multiple of $n\Rightarrow {\omega }^{4k}=1$, $$\begin{array}{rl}r(X)& =(a_0+a_4)+(a_1+a_5)X+(a_2+a_6)X^2+(a_3+a_7)X^3\\ r({\omega }^{2i})& =(a_0+a_4)+(a_1+a_5){\omega }^{2i}+(a_2+a_6){\omega }^{4i}+(a_3+a_7){\omega }^{6i}\\ f({\omega }^k)& =(a_0+a_1{\omega }^k+a_2{\omega }^{2k}+a_3{\omega }^{3k})+{\omega }^{4k}(a_4+a_5{\omega }^k+a_6{\omega }^{2k}+a_7{\omega }^{3k})\\ & =(a_0+a_1{\omega }^k+a_2{\omega }^{2k}+a_3{\omega }^{3k})+(a_4+a_5{\omega }^k+a_6{\omega }^{2k}+a_7{\omega }^{3k})\\ & =(a_0+a_4)+(a_1+a_5){\omega }^k+(a_2+a_6){\omega }^{2k}+(a_3+a_7){\omega }^{3k}\\ & =f({\omega }^{2i})\\ & =(a_0+a_4)+(a_1+a_5){\omega }^{2i}+(a_2+a_6){\omega }^{4i}+(a_3+a_7){\omega }^{6i}\\ & =r({\omega }^{2i})\end{array}$$

Odd Values

For $k=2i+1$ odd, we have a similar relation where $4k=8i+4$, so ${\omega }^{4k}={\omega }^4$. But observe that ${\omega }^4=-1$. $$\begin{array}{rl}s(X)& =(a_0-a_4)+\omega (a_1-a_5)X+{\omega }^2(a_2-a_6)X^2+{\omega }^3(a_3-a_7)X^3\\ s({\omega }^{2i})& =(a_0-a_4)+(a_1-a_5){\omega }^{2i+1}+(a_2-a_6){\omega }^{4i+2}+(a_3-a_7){\omega }^{6i+3}\\ f({\omega }^k)& =(a_0+a_1{\omega }^k+a_2{\omega }^{2k}+a_3{\omega }^{3k})+{\omega }^{4k}(a_4+a_5{\omega }^k+a_6{\omega }^{2k}+a_7{\omega }^{3k})\\ & =(a_0+a_1{\omega }^k+a_2{\omega }^{2k}+a_3{\omega }^{3k})-(a_4+a_5{\omega }^k+a_6{\omega }^{2k}+a_7{\omega }^{3k})\\ & =f({\omega }^{2i+1})\\ & =(a_0+a_1{\omega }^{2i+1}+a_2{\omega }^{4i+2}+a_3{\omega }^{6i+3})-(a_4+a_5{\omega }^{2i+1}+a_6{\omega }^{4i+2}+a_7{\omega }^{6i+3})\\ & =(a_0-a_4)+(a_1-a_5){\omega }^{2i+1}+(a_2-a_6){\omega }^{4i+2}+(a_3-a_7){\omega }^{6i+3}\\ & =s({\omega }^{2i})\end{array}$$