fft algo

doc/book.toml

@@ -18,6 +18,7 @@ runnable = false
 [output.katex]
 
 [preprocessor.katex]
+macros = "latex-macros.txt"
 
 [preprocessor.toc]
 command = "mdbook-toc"

doc/latex-macros.txt

@@ -0,0 +1 @@
+\DFT:{\textrm{DFT}}

doc/src/SUMMARY.md

@@ -37,6 +37,8 @@
   	- [Network Protocol](misc/hashchain/network_protocol.md)
   - [darkwiki](misc/darkwiki.md)
   - [dnetview](misc/dnetview.md)
+- [Crypto](crypto/crypto.md)
+  - [FFT](crypto/fft.md)
 - [Testnet Guide](testnet/guide.md)
 - [Learn](learn/learn.md)
   - [Research](learn/research.md)

doc/src/crypto/crypto.md

@@ -0,0 +1 @@
+# Crypto

doc/src/crypto/fft.md

@@ -0,0 +1,279 @@
+# Discrete Fast Fourier Transform
+
+Available code files:
+
+* [fft2.sage](https://github.com/darkrenaissance/darkfi/blob/master/script/research/zk/fft/fft2.sage):
+  implementation using vandermonde matrices illustrating the theory section below.
+* [fft3.sage](https://github.com/darkrenaissance/darkfi/blob/master/script/research/zk/fft/fft3.sage):
+  simple example with $n = 4$ showing 3 steps of the algorithm.
+* [fft4.sage](https://github.com/darkrenaissance/darkfi/blob/master/script/research/zk/fft/fft4.sage):
+  illustrates the full working algorithm.
+
+## Theory
+
+$$ f(x)g(x) ∈ 𝔽_{<2n}[x] $$
+$$ fg = \sum_{i + j < 2n - 2} a_i b_j x^{i + j} $$
+Complexity: $O(n^2)$
+
+Suppose $ω ∈ 𝔽$ is an nth root of unity.
+
+Recall: if $𝔽 = 𝔽_{p^k}$ then $∃N : 𝔽_{p^N}$
+contains all nth roots of unity.
+
+$$ \DFT_ω : 𝔽^n → 𝔽^n $$
+$$ \DFT_ω(f) = (f(ω^0), f(ω^1), …, f(ω^{n - 1})) $$
+
+$$ V_ω =
+\begin{pmatrix}
+1 & 1   & 1   & \cdots & 1 \\
+1 & ω^1 & ω^2 & \cdots & ω^{n - 1} \\
+1 & ω^2 & ω^4 & \cdots & ω^{2(n - 1)} \\
+\vdots \\
+1 & ω^{n - 1} & ω^{2(n - 1)} & \cdots & ω^{(n - 1)^2} \\
+\end{pmatrix}
+$$
+
+$$ \DFT_ω(f) = V_ω · f^T $$
+since vandermonde multiplication is simply evaluation of a polynomial.
+
+### Lemma: $V_ω^{-1} = \frac{1}{n} V_{ω^{-1}}$
+
+Use $1 + ω + ⋯ + ω^{n - 1}$ and compute $V_ω V_{ω^{-1}}$
+
+Corollary: $\DFT_ω$ is invertible.
+
+### Definitions
+
+1. Convolution $f * g = fg \mod (x^n - 1)$
+2. Pointwise product
+   $$(a_0, …, a_{n - 1})·(b_0, …, b_{n - 1}) =
+         (a_0 b_0, …, a_{n - 1} b_{n - 1}) ∈ 𝔽^n → 𝔽_{<n}[x]$$
+
+### Theorem: $\DFT_ω(f*g) = \DFT_ω(f)·\DFT_ω(g)$
+
+$$ fg = q'(x^n - 1) + f*g  $$
+$$ ⇒ f*g = fg + q(x^n - 1) $$
+$$ \deg fg ≤ 2n - 2        $$
+
+$$
+\begin{align}
+(f*g)(ω^i) &= f(ω^i)g(ω^i) + q(ω^i)(ω^{in} - 1) \\
+           &= f(ω^i)g(ω^i)
+\end{align}
+$$
+
+### Result
+
+$$ f, g ∈ 𝔽_{<n/2}[x] $$
+$$ fg = f*g $$
+$$ \DFT_ω(f*g) = \DFT_ω(f)·\DFT_ω(g) $$
+$$ fg = \frac{1}{n} \DFT_{ω^{-1}} (\DFT_ω(f) · \DFT_ω(g)) $$
+
+## Finite Field Extension Containing Nth Roots of Unity
+
+$$ μ_N = ⟨ω⟩, |𝔽_{p^N}^×| = p^N - 1 $$
+$$ \textrm{ord}(ω) = n | p^N - 1 $$
+but $𝔽_{p^N}^×$ is cyclic.
+
+For all $d | p^N - 1$, there exists $x ∈ 𝔽_{p^N}^×$
+with ord$(x) = d$.
+
+Finding $n | p^N - 1$ is sufficient for $ω ∈ 𝔽_{p^N}$
+$$ n | p^N - 1 ⇔ \textrm{ord}(p) = (ℤ / nℤ)^× $$
+
+## FFT Algorithm Recursive Compute
+
+We recurse to a depth of $\log n$. Since each recursion uses $ω^i$, then
+in the final step $ω^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{align}
+f(X)    &= a_0 + a_1 X + a_2 X^2 + ⋯ + a_{n - 1} X^{n - 1} \\
+        &= g(X) + X^{n/2} h(X)
+\end{align}
+$$
+
+### Algorithm
+
+Implementation of this algorithm is available in [fft4.sage](https://github.com/darkrenaissance/darkfi/blob/master/script/research/zk/fft/fft4.sage).
+Particularly the function called `calc_dft()`.
+
+**function** DFT($n = 2^d, f(X)$)<br>
+<span style="padding-left: 30px;">**if** $n = 1$ **then**</span><br>
+<span style="padding-left: 60px;">    **return** $f(X)$</span><br>
+<span style="padding-left: 30px;">**end**</span><br>
+<span style="padding-left: 30px;">Write $f(X)$ as the sum of two polynomials with equal degree</span><br>
+<span style="padding-left: 30px;">$f(X) = g(X) + X^{n/2} h(X)$</span><br>
+<span style="padding-left: 30px;">Let $\mathbf{g}, \mathbf{h}$ be the vector representations of $g(X), h(X)$</span><br>
+<span style="padding-left: 30px;"></span><br>
+<span style="padding-left: 30px;">$\mathbf{r} = \mathbf{g} + \mathbf{h}$</span><br>
+<span style="padding-left: 30px;">$\mathbf{s} = (\mathbf{g} - \mathbf{h})·(ω^0, …, ω^{n/2 - 1})$</span><br>
+<span style="padding-left: 30px;">Let $r(X), s(X)$ be the polynomials represented by the vectors $\mathbf{r}, \mathbf{s}$</span><br>
+<span style="padding-left: 30px;"></span><br>
+<span style="padding-left: 30px;">Compute $(r(ω^0), …, r(ω^{n/2})) = \textrm{DFT}_{ω^2}(n/2, r(X))$</span><br>
+<span style="padding-left: 30px;">Compute $(s(ω^0), …, s(ω^{n/2})) = \textrm{DFT}_{ω^2}(n/2, s(X))$</span><br>
+<span style="padding-left: 30px;"></span><br>
+<span style="padding-left: 30px;">**return** $(r(ω^0), s(ω^0), r(ω^2), s(ω^2), …, r(ω^{n/2}), s(ω^{n/2}))$</span><br>
+end
+
+Sage code:
+
+```python
+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{align}
+r(X)        &= g(X) + h(X) \\
+\\
+f(ω^{2i})   &= g(ω^{2i}) + (ω^{2i})^{n/2} h(ω^{2i}) \\
+            &= g(ω^{2i}) +                h(ω^{2i}) \\
+            &= (g + h)(ω^{2i}) \\
+\end{align}
+$$
+
+So then we can now compute $DFT_ω(f)_{k=2i} = DFT_{ω^2}(r)$
+for the even powers of $f(ω^{2i})$.
+
+### Odd Values
+
+For odd values $k = 2i + 1$
+
+$$
+\begin{align}
+s(X)        &= (g(X) - h(X))·(ω^0, …, ω^{n/2 - 1}) \\
+\\
+f(X)          &= a_0 + a_1 X + a_2 X^2 + ⋯ + a_{n - 1} X^{n - 1} \\
+              &= g(X) + X^{n/2} h(X) \\
+f(ω^{2i + 1}) &= g(ω^{2i + 1}) + (ω^{2i + 1})^{n/2} h(ω^{2i + 1}) \\
+\end{align}
+$$
+But observe that for any $n$th root of unity $ω^n = 1$ and $ω^{n/2} = -1$
+$$ (ω^{2i + 1})^{n/2} = ω^{in} ω^{n/2} = ω^{n/2} = -1 $$
+$$
+\begin{align}
+⇒ f(ω^{2i + 1}) &= g(ω^{2i + 1}) - h(ω^{2i + 1}) \\
+                &= (g - h)(ω^{2i + 1})
+\end{align}
+$$
+
+Let $\mathbf{s} = (\mathbf{g} - \mathbf{h})·(ω^0, …, ω^{n/2 - 1})$ be the
+representation for $s(X)$. Then we can see that $s(ω^{2i}) = (g - h)(ω^{2i + 1})$
+as desired.
+
+So then we can now compute $DFT_ω(f)_{k=2i + 1} = DFT_{ω^2}(s)$
+for the odd powers of $f(ω^{2i + 1})$.
+
+## Example
+
+Let $n = 8$
+$$
+\begin{align}
+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{align}
+$$
+Now vectorize $g(X), h(X)$
+$$
+\begin{align}
+\mathbf{g} &= (a_0, a_1, a_2, a_3) \\
+\mathbf{h} &= (a_4, a_5, a_6, a_7) \\
+\end{align}
+$$
+Compute reduced polynomials in vector form
+$$
+\begin{align}
+\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})·(1, ω, ω^2, ω^3) \\
+           &= (a_0 - a_4, a_1 - a_5, a_2 - a_6, a_3 - a_7)·(1, ω, ω^2, ω^3) \\
+           &= (a_0 - a_4, ω (a_1 - a_5), ω^2 (a_2 - a_6), ω^3 (a_3 - a_7)) \\
+\end{align}
+$$
+Convert them to polynomials from the vectors. We also expand them out below
+for completeness.
+$$
+\begin{align}
+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) + ω (a_1 - a_5) X + ω^2 (a_2 - a_6) X^2 + ω^3 (a_3 - a_7) X^3 \\
+\end{align}
+$$
+Compute
+$$ \textrm{DFT}_{ω^2}(4, r(X)), \textrm{DFT}_{ω^2}(4, s(X)) $$
+The values returned will be
+$$
+(r(1), s(1), r(ω^2), s(ω^2), r(ω^4), s(ω^4), r(ω^6), s(ω^6))
+=
+(f(1), f(ω), f(ω^2), f(ω^3), f(ω^4), f(ω^5), f(ω^6), f(ω^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 $ω^i$.
+
+We expect that $s(X)$ on the domain $(1, ω^2, ω^4, ω^6)$ produces the values
+$(f(1), f(ω^2), f(ω^4), f(ω^6))$, while $r(X)$ on the same domain produces
+$(f(ω), f(ω^3), f(ω^5), f(ω^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 ⇒ ω^{4k} = 1$,
+$$
+\begin{align}
+r(X)       &= (a_0 + a_4) + (a_1 + a_5) X + (a_2 + a_6) X^2 + (a_3 + a_7) X^3 \\
+r(ω^{2i})  &= (a_0 + a_4) + (a_1 + a_5) ω^{2i} + (a_2 + a_6) ω^{4i} + (a_3 + a_7) ω^{6i} \\
+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}) \\
+           &= (a_0 + a_1 ω^k + a_2 ω^{2k} + a_3 ω^{3k}) +        (a_4     + a_5 ω^k   + a_6 ω^{2k} + a_7 ω^{3k}) \\
+           &= (a_0 + a_4) + (a_1 + a_5) ω^k + (a_2 + a_6) ω^{2k} + (a_3 + a_7) ω^{3k} \\
+           &= f(ω^{2i}) \\
+           &= (a_0 + a_4) + (a_1 + a_5) ω^{2i} + (a_2 + a_6) ω^{4i} + (a_3 + a_7) ω^{6i} \\
+           &= r(ω^{2i})
+\end{align}
+$$
+
+#### Odd Values
+
+For $k = 2i + 1$ odd, we have a similar relation where $4k = 8i + 4$, so
+$ω^{4k} = ω^4$. But observe that $ω^4 = -1$.
+$$
+\begin{align}
+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 \\
+s(ω^{2i})  &= (a_0 - a_4) + (a_1 - a_5) ω^{2i + 1} + (a_2 - a_6) ω^{4i + 2} + (a_3 - a_7) ω^{6i + 3} \\
+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}) \\
+           &= (a_0 + a_1 ω^k + a_2 ω^{2k} + a_3 ω^{3k}) -        (a_4     + a_5 ω^k   + a_6 ω^{2k} + a_7 ω^{3k}) \\
+           &= f(ω^{2i + 1}) \\
+           &= (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}) \\
+           &= 
+           (a_0  - a_4)
+           + (a_1 - a_5) ω^{2i + 1}
+           + (a_2 - a_6) ω^{4i + 2}
+           + (a_3 - a_7) ω^{6i + 3} \\
+           &= s(ω^{2i})
+\end{align}
+$$
+