1  # Discrete Fast Fourier Transform
  2  
  3  Available code files:
  4  
  5  * [fft2.sage](https://codeberg.org/darkrenaissance/darkfi/src/branch/master/script/research/zk/fft/fft2.sage):
  6    implementation using vandermonde matrices illustrating the theory section below.
  7  * [fft3.sage](https://codeberg.org/darkrenaissance/darkfi/src/branch/master/script/research/zk/fft/fft3.sage):
  8    simple example with $n = 4$ showing 3 steps of the algorithm.
  9  * [fft4.sage](https://codeberg.org/darkrenaissance/darkfi/src/branch/master/script/research/zk/fft/fft4.sage):
 10    illustrates the full working algorithm.
 11  
 12  ## Theory
 13  
 14  $$ f(x)g(x) ∈ 𝔽_{<2n}[x] $$
 15  $$ fg = \sum_{i + j < 2n - 2} a_i b_j x^{i + j} $$
 16  Complexity: $O(n^2)$
 17  
 18  Suppose $ω ∈ 𝔽$ is an nth root of unity.
 19  
 20  Recall: if $𝔽 = 𝔽_{p^k}$ then $∃N : 𝔽_{p^N}$
 21  contains all nth roots of unity.
 22  
 23  $$ \DFT_ω : 𝔽^n → 𝔽^n $$
 24  $$ \DFT_ω(f) = (f(ω^0), f(ω^1), …, f(ω^{n - 1})) $$
 25  
 26  $$ V_ω =
 27  \begin{pmatrix}
 28  1 & 1   & 1   & \cdots & 1 \\
 29  1 & ω^1 & ω^2 & \cdots & ω^{n - 1} \\
 30  1 & ω^2 & ω^4 & \cdots & ω^{2(n - 1)} \\
 31  \vdots \\
 32  1 & ω^{n - 1} & ω^{2(n - 1)} & \cdots & ω^{(n - 1)^2} \\
 33  \end{pmatrix}
 34  $$
 35  
 36  $$ \DFT_ω(f) = V_ω · f^T $$
 37  since vandermonde multiplication is simply evaluation of a polynomial.
 38  
 39  ### Lemma: $V_ω^{-1} = \frac{1}{n} V_{ω^{-1}}$
 40  
 41  Use $1 + ω + ⋯ + ω^{n - 1}$ and compute $V_ω V_{ω^{-1}}$
 42  
 43  Corollary: $\DFT_ω$ is invertible.
 44  
 45  ### Definitions
 46  
 47  1. Convolution $f * g = fg \mod (x^n - 1)$
 48  2. Pointwise product
 49     $$(a_0, …, a_{n - 1})·(b_0, …, b_{n - 1}) =
 50           (a_0 b_0, …, a_{n - 1} b_{n - 1}) ∈ 𝔽^n → 𝔽_{<n}[x]$$
 51  
 52  ### Theorem: $\DFT_ω(f*g) = \DFT_ω(f)·\DFT_ω(g)$
 53  
 54  $$ fg = q'(x^n - 1) + f*g  $$
 55  $$ ⇒ f*g = fg + q(x^n - 1) $$
 56  $$ \deg fg ≤ 2n - 2        $$
 57  
 58  $$
 59  \begin{align}
 60  (f*g)(ω^i) &= f(ω^i)g(ω^i) + q(ω^i)(ω^{in} - 1) \\
 61             &= f(ω^i)g(ω^i)
 62  \end{align}
 63  $$
 64  
 65  ### Result
 66  
 67  $$ f, g ∈ 𝔽_{<n/2}[x] $$
 68  $$ fg = f*g $$
 69  $$ \DFT_ω(f*g) = \DFT_ω(f)·\DFT_ω(g) $$
 70  $$ fg = \frac{1}{n} \DFT_{ω^{-1}} (\DFT_ω(f) · \DFT_ω(g)) $$
 71  
 72  ## Finite Field Extension Containing Nth Roots of Unity
 73  
 74  $$ μ_N = ⟨ω⟩, |𝔽_{p^N}^×| = p^N - 1 $$
 75  $$ \textrm{ord}(ω) = n | p^N - 1 $$
 76  but $𝔽_{p^N}^×$ is cyclic.
 77  
 78  For all $d | p^N - 1$, there exists $x ∈ 𝔽_{p^N}^×$
 79  with ord$(x) = d$.
 80  
 81  Finding $n | p^N - 1$ is sufficient for $ω ∈ 𝔽_{p^N}$
 82  $$ n | p^N - 1 ⇔ \textrm{ord}(p) = (ℤ / nℤ)^× $$
 83  
 84  ## FFT Algorithm Recursive Compute
 85  
 86  We recurse to a depth of $\log n$. Since each recursion uses $ω^i$, then
 87  in the final step $ω^i = 1$, and we simply return $f^T$.
 88  
 89  We only need to prove a single step of the algorithm produces the desired
 90  result, and then the correctness is inductively proven.
 91  
 92  $$
 93  \begin{align}
 94  f(X)    &= a_0 + a_1 X + a_2 X^2 + ⋯ + a_{n - 1} X^{n - 1} \\
 95          &= g(X) + X^{n/2} h(X)
 96  \end{align}
 97  $$
 98  
 99  ### Algorithm
100  
101  Implementation of this algorithm is available in [fft4.sage](https://codeberg.org/darkrenaissance/darkfi/src/branch/master/script/research/zk/fft/fft4.sage).
102  Particularly the function called `calc_dft()`.
103  
104  **function** DFT($n = 2^d, f(X)$)<br>
105  <span style="padding-left: 30px;">**if** $n = 1$ **then**</span><br>
106  <span style="padding-left: 60px;">    **return** $f(X)$</span><br>
107  <span style="padding-left: 30px;">**end**</span><br>
108  <span style="padding-left: 30px;">Write $f(X)$ as the sum of two polynomials with equal degree</span><br>
109  <span style="padding-left: 30px;">$f(X) = g(X) + X^{n/2} h(X)$</span><br>
110  <span style="padding-left: 30px;">Let $\mathbf{g}, \mathbf{h}$ be the vector representations of $g(X), h(X)$</span><br>
111  <span style="padding-left: 30px;"></span><br>
112  <span style="padding-left: 30px;">$\mathbf{r} = \mathbf{g} + \mathbf{h}$</span><br>
113  <span style="padding-left: 30px;">$\mathbf{s} = (\mathbf{g} - \mathbf{h})·(ω^0, …, ω^{n/2 - 1})$</span><br>
114  <span style="padding-left: 30px;">Let $r(X), s(X)$ be the polynomials represented by the vectors $\mathbf{r}, \mathbf{s}$</span><br>
115  <span style="padding-left: 30px;"></span><br>
116  <span style="padding-left: 30px;">Compute $(r(ω^0), …, r(ω^{n/2})) = \textrm{DFT}_{ω^2}(n/2, r(X))$</span><br>
117  <span style="padding-left: 30px;">Compute $(s(ω^0), …, s(ω^{n/2})) = \textrm{DFT}_{ω^2}(n/2, s(X))$</span><br>
118  <span style="padding-left: 30px;"></span><br>
119  <span style="padding-left: 30px;">**return** $(r(ω^0), s(ω^0), r(ω^2), s(ω^2), …, r(ω^{n/2}), s(ω^{n/2}))$</span><br>
120  end
121  
122  Sage code:
123  
124  ```python
125  def calc_dft(ω_powers, f):
126      m = len(f)
127      if m == 1:
128          return f
129      g, h = vector(f[:m/2]), vector(f[m/2:])
130  
131      r = g + h
132      s = dot(g - h, ω_powers)
133  
134      ω_powers = vector(ω_i for ω_i in ω_powers[::2])
135      rT = calc_dft(ω_powers, r)
136      sT = calc_dft(ω_powers, s)
137  
138      return list(alternate(rT, sT))
139  ```
140  
141  ### Even Values
142  
143  $$
144  \begin{align}
145  r(X)        &= g(X) + h(X) \\
146  \\
147  f(ω^{2i})   &= g(ω^{2i}) + (ω^{2i})^{n/2} h(ω^{2i}) \\
148              &= g(ω^{2i}) +                h(ω^{2i}) \\
149              &= (g + h)(ω^{2i}) \\
150  \end{align}
151  $$
152  
153  So then we can now compute $DFT_ω(f)_{k=2i} = DFT_{ω^2}(r)$
154  for the even powers of $f(ω^{2i})$.
155  
156  ### Odd Values
157  
158  For odd values $k = 2i + 1$
159  
160  $$
161  \begin{align}
162  s(X)        &= (g(X) - h(X))·(ω^0, …, ω^{n/2 - 1}) \\
163  \\
164  f(X)          &= a_0 + a_1 X + a_2 X^2 + ⋯ + a_{n - 1} X^{n - 1} \\
165                &= g(X) + X^{n/2} h(X) \\
166  f(ω^{2i + 1}) &= g(ω^{2i + 1}) + (ω^{2i + 1})^{n/2} h(ω^{2i + 1}) \\
167  \end{align}
168  $$
169  But observe that for any $n$th root of unity $ω^n = 1$ and $ω^{n/2} = -1$
170  $$ (ω^{2i + 1})^{n/2} = ω^{in} ω^{n/2} = ω^{n/2} = -1 $$
171  $$
172  \begin{align}
173  ⇒ f(ω^{2i + 1}) &= g(ω^{2i + 1}) - h(ω^{2i + 1}) \\
174                  &= (g - h)(ω^{2i + 1})
175  \end{align}
176  $$
177  
178  Let $\mathbf{s} = (\mathbf{g} - \mathbf{h})·(ω^0, …, ω^{n/2 - 1})$ be the
179  representation for $s(X)$. Then we can see that $s(ω^{2i}) = (g - h)(ω^{2i + 1})$
180  as desired.
181  
182  So then we can now compute $DFT_ω(f)_{k=2i + 1} = DFT_{ω^2}(s)$
183  for the odd powers of $f(ω^{2i + 1})$.
184  
185  ## Example
186  
187  Let $n = 8$
188  $$
189  \begin{align}
190  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) \\
191          &= (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) \\
192          &= g(X) + X^{n/2} h(X) \\
193  g(X)    &=  a_0 + a_1 X + a_2 X^2 + a_3 X^3 \\
194  h(X)    &=  a_4 + a_5 X + a_6 X^2 + a_7 X^3 \\
195  \end{align}
196  $$
197  Now vectorize $g(X), h(X)$
198  $$
199  \begin{align}
200  \mathbf{g} &= (a_0, a_1, a_2, a_3) \\
201  \mathbf{h} &= (a_4, a_5, a_6, a_7) \\
202  \end{align}
203  $$
204  Compute reduced polynomials in vector form
205  $$
206  \begin{align}
207  \mathbf{r} &=  \mathbf{g} + \mathbf{h} \\
208             &= (a_0 + a_4, a_1 + a_5, a_2 + a_6, a_3 + a_7) \\
209  \mathbf{s} &= (\mathbf{g} - \mathbf{h})·(1, ω, ω^2, ω^3) \\
210             &= (a_0 - a_4, a_1 - a_5, a_2 - a_6, a_3 - a_7)·(1, ω, ω^2, ω^3) \\
211             &= (a_0 - a_4, ω (a_1 - a_5), ω^2 (a_2 - a_6), ω^3 (a_3 - a_7)) \\
212  \end{align}
213  $$
214  Convert them to polynomials from the vectors. We also expand them out below
215  for completeness.
216  $$
217  \begin{align}
218  r(X)       &= r_0 + r_1 X + r_2 X^2 + r_3 X^3 \\
219             &= (a_0 + a_4) + (a_1 + a_5) X + (a_2 + a_6) X^2 + (a_3 + a_7) X^3 \\
220  s(X)       &= s_0 + s_1 X + s_2 X^2 + s_3 X^3 \\
221             &= (a_0 - a_4) + ω (a_1 - a_5) X + ω^2 (a_2 - a_6) X^2 + ω^3 (a_3 - a_7) X^3 \\
222  \end{align}
223  $$
224  Compute
225  $$ \textrm{DFT}_{ω^2}(4, r(X)), \textrm{DFT}_{ω^2}(4, s(X)) $$
226  The values returned will be
227  $$
228  (r(1), s(1), r(ω^2), s(ω^2), r(ω^4), s(ω^4), r(ω^6), s(ω^6))
229  =
230  (f(1), f(ω), f(ω^2), f(ω^3), f(ω^4), f(ω^5), f(ω^6), f(ω^7))
231  $$
232  Which is the output we return.
233  
234  ### Comparing Evaluations for $f(X)$ and $r(X), s(X)$
235  
236  We can see the evaluations are correct by substituting in $ω^i$.
237  
238  We expect that $s(X)$ on the domain $(1, ω^2, ω^4, ω^6)$ produces the values
239  $(f(1), f(ω^2), f(ω^4), f(ω^6))$, while $r(X)$ on the same domain produces
240  $(f(ω), f(ω^3), f(ω^5), f(ω^7))$.
241  
242  #### Even Values
243  
244  Let $k = 2i$, be an even number. Then note that $k$ is a multiple of 2, so
245  $4k$ is a multiple of $n ⇒ ω^{4k} = 1$,
246  $$
247  \begin{align}
248  r(X)       &= (a_0 + a_4) + (a_1 + a_5) X + (a_2 + a_6) X^2 + (a_3 + a_7) X^3 \\
249  r(ω^{2i})  &= (a_0 + a_4) + (a_1 + a_5) ω^{2i} + (a_2 + a_6) ω^{4i} + (a_3 + a_7) ω^{6i} \\
250  f(ω^k)     &= (a_0 + a_1 ω^k + a_2 ω^{2k} + a_3 ω^{3k}) + ω^{4k} (a_4     + a_5 ω^k   + a_6 ω^{2k} + a_7 ω^{3k}) \\
251             &= (a_0 + a_1 ω^k + a_2 ω^{2k} + a_3 ω^{3k}) +        (a_4     + a_5 ω^k   + a_6 ω^{2k} + a_7 ω^{3k}) \\
252             &= (a_0 + a_4) + (a_1 + a_5) ω^k + (a_2 + a_6) ω^{2k} + (a_3 + a_7) ω^{3k} \\
253             &= f(ω^{2i}) \\
254             &= (a_0 + a_4) + (a_1 + a_5) ω^{2i} + (a_2 + a_6) ω^{4i} + (a_3 + a_7) ω^{6i} \\
255             &= r(ω^{2i})
256  \end{align}
257  $$
258  
259  #### Odd Values
260  
261  For $k = 2i + 1$ odd, we have a similar relation where $4k = 8i + 4$, so
262  $ω^{4k} = ω^4$. But observe that $ω^4 = -1$.
263  $$
264  \begin{align}
265  s(X)       &= (a_0 - a_4) + ω (a_1 - a_5) X + ω^2 (a_2 - a_6) X^2 + ω^3 (a_3 - a_7) X^3 \\
266  s(ω^{2i})  &= (a_0 - a_4) + (a_1 - a_5) ω^{2i + 1} + (a_2 - a_6) ω^{4i + 2} + (a_3 - a_7) ω^{6i + 3} \\
267  f(ω^k)     &= (a_0 + a_1 ω^k + a_2 ω^{2k} + a_3 ω^{3k}) + ω^{4k} (a_4     + a_5 ω^k   + a_6 ω^{2k} + a_7 ω^{3k}) \\
268             &= (a_0 + a_1 ω^k + a_2 ω^{2k} + a_3 ω^{3k}) -        (a_4     + a_5 ω^k   + a_6 ω^{2k} + a_7 ω^{3k}) \\
269             &= f(ω^{2i + 1}) \\
270             &= (a_0 + a_1 ω^{2i + 1} + a_2 ω^{4i + 2} + a_3 ω^{6i + 3}) -  (a_4     + a_5 ω^{2i + 1}   + a_6 ω^{4i + 2} + a_7 ω^{6i + 3}) \\
271             &= 
272             (a_0  - a_4)
273             + (a_1 - a_5) ω^{2i + 1}
274             + (a_2 - a_6) ω^{4i + 2}
275             + (a_3 - a_7) ω^{6i + 3} \\
276             &= s(ω^{2i})
277  \end{align}
278  $$
279