# Tutorial 4: The Numeric Theoretic Transform.
# The Numeric Theoretic Transform (NTT).
In the previous tutorial [tut3](https://hackmd.io/@qozymandias/HJFqNUHBa) we covered how to convert programs into polynomials and use this representation to prove correctness. Now we must consider how to efficiently perform operations on polynomials such as "polynomial multiplication", which is a common operation, and we must consider how to do this efficiently. To understand this, we must first reacquaint ourselves with the Fast Fourier transform and learn how to apply it to our Finite Fields.

The Numeric Theoretic Transform (NTT) is a specialisation of the Discrete Fourier Transform (DFT) for integral numbers. The Fast Fourier Transform (FFT) is an efficient algorithm for computing the DFT, and therefore we will first learn about the FFT and then understand how to tweak it for integral numbers.

In the following sections, we take some time to explain the FFT before we explain the NTT. If you already understand this, then feel free to skip ahead to the NTT section.

# Understanding the FFT.
This has been detailed on another page: [Tutorial 4 - Understanding the FFT](https://hackmd.io/@qozymandias/B1sv7DoFT)

# Understanding the NTT.
## Introduction
NTT uses the same underlying FFT algorithm, however, the key difference is the number type. FFT is designed for decimal complex numbers, whereas, NTT is specifically designed for integral numbers, which means no fractional numbers and no imaginary numbers! This seems quite impossible, but, there's always a way with mathematics! The NTT achieves by leveraging finite fields! We have previously covered finite fields, but I will cover the essentials here, for more details I would point to this summary web page, [notes on rings and fields](https://hackmd.io/@YW50b255-base64/ry4jiceN6).

But first, why is it so important? Simply put, we need to perform efficient polynomial multiplication in our proving system. The NTT forward transform converts coefficient form into evaluation form, then multiplication is simply down point-wise; and finally, to get the result in coefficient form the inverse NTT is applied. This comes naturally from the fact that NTT is a variant of the FFT. The NTT uses finite field numbers, meaning there are no rounding errors because we don't allow fractional numbers (rather every nonzero element has a perfectly whole multiplicative inverse). This is essential for our polynomial proving system. Additionally, integer operations are more efficient than floating point operations and certain operations can be further optimised. The NTT is the perfect tool for the proving system, so let's get into the details that set it apart from the standard FFT.

## Understanding the maths behind NTT.
The most important part to understand of the NTT is the pre-processing. However, there are two scenarios to consider: using an existing finite field such as an elliptic curve, or "constructing" a finite field from the inputs. The former, scenario is straightforward, however, the latter is more involved, and we will focus on this to build understanding. By "constructing" a finite field, I mean selecting a prime number that suites the inputs, more on this soon, but first, an aside about finite fields for the sake of clarity.

### Prime Finite Fields versus Finite Fields.
NTT is FFT specialised for finite fields. A finite field is a finite ring with the additional property that every non-zero element has a multiplicative inverse, this is also known as a quotient ring. I have to make all the vocabulary clear here because it can get quite confusing; when they say quotient ring it is that same as finite field, additionally, I have to make a distinction between finite fields and prime finite fields. [Notes on rings and fields](https://hackmd.io/@YW50b255-base64/ry4jiceN6) does a better explanation, so I will quickly outline the difference, and you may refer to that webpage for more details.

A finite field, also known as a Galois field, is a field that contains a finite number of elements. The number of elements in a finite field is either a prime number or a power of a prime number (i.e. $p^k$ where $p$ is a prime number and $k$ is a positive integer). For every prime number $p$ and positive integer $n$, there exists a finite field with $p^n$ elements.

A prime finite field is a specific type of finite field where the number of elements (also known as the order of the field) is a prime number. The most common examples of prime finite fields are given by the integers modulo $p$ when $p$ is a prime number. These are often denoted as $Z/p$ or $F_p$.

NTT can use a finite field provided it chooses a suitable modulus and primitive root of unity, which both fall out easily because of the properties of a finite field.
Technically the modulus doesn't have to be prime, [see here](https://en.wikipedia.org/wiki/Discrete_Fourier_transform_over_a_ring#Number-theoretic_transform), but we will only focus on prime finite fields for simplicity.
### Intuition.
Now we have some idea about prime finite fields, and we can come up with a mental model to aid our intuition. Take a set of non-negative integers, determine the largest integer, and then calculate the next integer that is greater than this integer and is a prime, let this be $p$. Now, apply $p$ as a modulo operation to every integer in the set, and we have constructed a prime finite field! Remember, applying a modulo operation is simply calculating the remain of the input when divided by the modulo -- this means every integer remain the same value, however any operations with "wrap around". The field is the set of integers from $0$ to $p-1$ (inclusive), and now we must perform addition, subtraction, and multiplication operations under modulo $p$. Division can be a bit tricky: essentially, instead of raw division, the multiplicative inverse is used. Let's unravel the operations with an example.

#### Example $GF(5)$.
For example, if $p = 5$, the prime finite field would consist of the integers $0, 1, 2, 3, 4$. If we were to add $3$ and $4$ in this field, we would get $2$ (because $3+4=7$, and $7$ modulo $5$ is $2$). Naturally, you can now see how we would derive the field from a set of input integers by simply choosing the next largest prime integer from the given integers -- for example, if our input was $6, 11$ then our chosen prime would be $13$ and the field would be $GF(13)$, which would be made up over the integers $0,1,2,3,4,5,6,7,8,9,10,11,12$.

The notation for the finite field is $GF(5)$ which means Galois Field (aka a finite field) of order $5$.

##### $GF(5)$ Addition, subtraction, multiplication.
So, addition and multiplication can be perfectly defined by applying modulo to the results, this ensures the results don't "overflow" the field. Subtraction on the other hand may lead to "overflow" in the negative direction, this is only a problem when the underlying type is unsigned, so we add the modulus onto the result of the subtraction before calculating its modulo. For example, if the underlying type is signed then simply calculating $(3 - 4) \space \text{ mod } \space 5$ is sufficient, however, if the type is unsigned, then we must calculate instead $((3 - 4) + 5) \space \text{ mod } \space 5$ which provides the same result. I'm getting ahead of myself, but this is very important to note for the implementation side of things.

##### $GF(5)$ Division with the multiplicative inverse.
It is also critical to mention the multiplicative inverse, essentially given the input $a$ we must find a $b$ such that the multiplication of $a$ and $b$ equals $1$ (modulo $p$):
$$
(a * b) \space\text{ mod }\space p = 1
$$
For example, say we want to divide $3$ by $2$, first, we need to find the multiplicative inverse of $2$ in $GF(5)$. So we need to find $b$ in $GF(5)$ such that the multiplicative inverse equation is satisfied:
$$
(2 * b) \space\text{ mod }\space 5 = 1
$$
By brute force checking this equation with the elements of the field, we find $b=3$ to satisfy the equation, since $2*3 = 6$, and $6$ modulo $5$ is $1$. So, to perform the division of $3$ by $2$, we multiply $3$ by the multiplicative inverse of $2$, which is $3$. Therefore, $3/2 = 3*3 = 9$. But since we’re in $GF(5)$, we take this result modulo $5$, so $9$ modulo $5$ is $4$. Therefore, in the finite field $GF(5)$, $3/2 = 4$.

For more complex fields, we can use [extended Euclidean algorithm](https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm "Extended Euclidean algorithm"). However, there is a trick (which we will explain later) that allows us to directly calculate the multiplicative inverse with using the power operator. The inverse of $a$ is:
$$
\frac{1}{a} = a^{-1}
$$
However, we cannot use division or calculate negative powers directly. So instead we must, employ [Fermat's Little Theorem](https://en.wikipedia.org/wiki/Fermat%27s_little_theorem) as a workaround (since the nonzero elements of $GF(p)$ form a [finite group](https://en.wikipedia.org/wiki/Finite_group "Finite group") with respect to multiplication):
$$
\begin{align}
a^{p − 1} \space\text{ mod }\space p \equiv 1 && (\text{for }\space a \ne 0) \\
\end{align}
$$
Now multiply both sides by $a^{-1}$
$$
\begin{align}
a^{p − 1} * a^{-1} \space\text{ mod }\space p \equiv 1 * a^{-1}
\end{align}
$$
Thus the inverse of $a$ is
$$
\begin{align}
a^{p − 2} \space\text{ mod }\space p && (\text{where } \space p = P^N) \\
\end{align}
$$
$P$ is a prime number (the characteristic of the field) and $N$ is a positive integer (the dimension of the field).

If we return to our original example, calculating the inverse of $2$:
$$
(2 * b) \space\text{ mod }\space 5 = 1
$$
Then:
$$
\begin{align}
\text{inverse}(2) &= 2^{5 − 2} \space\text{ mod }\space 5 \\
&= 3 \\
\end{align}
$$
Which is the same result as our brute force method!

#### $n$'th primitive root of unity.
##### Recap of complex $n$'th primitive root of unity.
Let's continue to build our intuition. We have just defined our prime field, $GF(p)$, where our elements are drawn from the set $[0, p - 1]$. Now, let's figure out how to derive the primitive root of unity. In the FFT section, the complex $n$'th roots of unity, are the $n$ complex solutions to the equation:
$$
w^n = 1
$$
Where the primitive root of unity would be:
$$
w_n = e^{2\pi{}i/n}
$$
And we could calculate the $n$'th roots of unity from the primitive root:
$$
1, w^1, w^2, w^3, \ldots, w^{n-2}, w^{n-1}
$$
This makes sense for the unit circle because every point is on the circumference of the circle, so its "real" part would always equal 1 (aka the modulus) where the angle changes as the points go around the curve (where the angle is evenly spaced by $\frac{2\pi}{n}$).

##### Primitive root of unity of finite field.
Now let's think about it within the context of the finite field, well first off we are mimicking the unit circle used by complex numbers by constructing a finite field with a prime modulus. By this I mean we have essentially taken a set of integers and "bent" them around into a circle, thus making them a cyclic group; that is, they "wrap around" like a clock thanks to the modulo! Previously, I have shown how determining a minimum prime modulus and applying modulo arithmetic for the set of integers creates a prime finite field, and thus we are essentially "rolling around" in a circle!

This means we have formed a cyclic [finite group](https://en.wikipedia.org/wiki/Finite_group "Finite group") with respect to multiplication. Because all nonzero elements can be expressed as powers of a single element, called a primitive element of the field -- this is the generator of the multiplicative group. This we have touched on before in terms of an elliptic curve field where $g$ is a generator for the group and every element is some power of $g$.

For example, the prime finite field $GF(5)$ has a generator of $2$ and $3$. Because a generator needs to be able to generate all the nonzero elements when raised to different powers:
$$
\begin{align}
2^1 \space\text{ mod }\space 5 &= 2 \\
2^2 \space\text{ mod }\space 5 &= 4 \\
2^3 \space\text{ mod }\space 5 &= 3 \\
2^4 \space\text{ mod }\space 5 &= 1 \\
\end{align}
$$
This is an important realisation because the complex $n$'th roots of unity are also a cyclic finite group! Now, let's think about what could be the generator for the complex numbers. The group in that regard would simply be the $n$ roots of unity:
$$
1, w^1, w^2, w^3, \ldots, w^{n-2}, w^{n-1}
$$
These are given by the formula:
$$
\begin{align}
w &= e^{2\pi{}k \space i/n} \\
&& \text{where } k=0, 1, \ldots, n-1 \\
\end{align}
$$
Notice we take $e^{2\pi i / n}$ to powers $0,1,2, \ldots, n-1$ and we can generate all the elements of the group! Thus, the generator $g$ of a finite field is playing the role of the $e^{2\pi i / n}$ for complex numbers.

Now it's important to note, the $n$'th primitive root of unity is closely tied to the structure of the finite field. Let's be clear, the complex $n$th roots of unity form this group, so the primitive root is the generator (aka primitive element). However, for a finite field we have just constructed we don't have such a nice property, rather, we have to brute force check our group against the definition of the primitive root.

> **Definition Primitive root of unity of finite field.**
> If $n$ is a positive integer, an $n$'th primitive root of unity is a solution to the equation:
> $$
x^n = 1
$$
and where
> $$
\begin{align}
x^m & \ne 1 \\
&& (\text{where }\space 0 < m < n) \\
\end{align}
$$
If $a$ is an $n$'th primitive root of unity in a field $F$, then $F$ contains all the $n$ roots of unity, which are:
> $$
1, a^1, a^2, \ldots, a^{n-1}
$$
The field $GF(p)$ contains an $n$'th primitive root of unity if and only if $n$ is a divisor of $p-1$; if $n$ is a divisor of $p-1$, then the number of primitive $n$'th roots of unity in $GF(p)$ is $\varphi(n)$ (Euler's totient function). The number of $n$'th roots of unity in $GF(p)$ is $gcd(n,\space p -1 )$

I don't know about you, but that just goes over my head, so let's try to break it down. Let's go back to the idea of the generator, $g$, for our field $GF(p)$, we can define the field elements as powers of $g$ (modulo $p$):
$$
(g^0)\text{ mod }\space p \space,
(g^1)\text{ mod }\space p \space,
(g^2)\text{ mod }\space p \space,
\ldots,
(g^{p-1})\text{ mod }\space p
$$
So our generator will be able to generate all the nonzero elements, that is, $p - 1$ distinct elements. But we know this just maps to all the integers in our set because we choose $g$ for this very purpose.
$$
1, 2, 3, 4, \ldots, (p-1)
$$
Let's remember that this is the actual set. By the way we constructed this finite field, we know all the elements will be not divisible by $p$ because we selected $p$ to be the next largest prime integer. Therefore, we can use [Fermat's Little Theorem](https://en.wikipedia.org/wiki/Fermat%27s_little_theorem), which states that states if $p$ is a prime number and $x$ is an integer that is not divisible by $p$, then $x^{p-1}$ will leave a remainder of $1$ when divided by $p$, in other words, we have the following equivalence:
$$
x^{p-1} \space\text{ mod }\space p \equiv 1
$$
This brings us close to the original root of unity equation. Recall, in the FFT section, the complex $n$th roots of unity, are the $n$ complex solutions to the equation:
$$
w^n = 1
$$
Where the primitive root of unity is:
$$
w_n = e^{2\pi{}i/n}
$$
Let's plug this into the above equation:
$$
\begin{align}
(e^{2\pi{}i/n})^n &= e^{2\pi i } \\
&= 1 \\
\end{align}
$$
Which makes sense when considering the polar form. However, our finite field equation is slightly different from the complex definition:
$$
w^n = 1
$$
Let's start with our new equivalence from Fermat's little theorem:
$$
x^{p-1} \space\text{ mod }\space p \equiv 1
$$
Let's replace $x$ with a value from our set, let's use the generator, aka the primitive element.
$$
g^{p-1} \space\text{ mod }\space p \equiv 1
$$
In order to fit the form of the root of unity equation we have to remember the element is taken to the power of $n$, $w^n = 1$, therefore we have to account for this in the equation. So, in order for this to work we have to divide the exponent by $n$ (how many $n$ roots of unity we desire) $n$ is determined by the number of coefficients in the original polynomial.
$$
w = g^{ ( p - 1) / n}
$$
Plugging this into $w^n = 1$
$$
\begin{align}
(g^{ ( p - 1) / n})^n &= 1 \\
g^{p-1} &= 1 \\
\end{align}
$$
Which we know is true because of Fermat's little theorem. Therefore, the $n$'th primitive root of unity can be defined by the generator:
$$
w_n = g^{ ( p - 1) / n}
$$
Where $n$ will be the length of the original polynomial. Now we have "kicked the can down the road", all the difficulty lies in determining the $g$ of our field, this is tricky, in the next section I have outlined how to do it algorithmically.

### Formal procedure
The following steps have been paraphrased from [nayuki's excellent implementation and explanation](https://www.nayuki.io/page/number-theoretic-transform-integer-dft).

Given a polynomial with $n$ coefficients:
$$
P = c_0 x^0 + c_1 x^1 + \ldots + c_{n-1} x^{n-1}
$$
1. Choose a minimum modulus $m$ that is prime, such that $1\le n \le m$ and every coefficient $\{c_0, c_1, \ldots, c_{n-1}\}$ exists in the range $[0, m)$
2. Determine the generator $g$ within the range $[1, m)$, in $GF(m)$, an element $a \in GF(m)$ is called a primitive element if it is a primitive $(m − 1)$'th root of unity in $GF(m)$
for each $i$ in $[1, m)$ {
if $i^{m-1} \space\text{ mod }\space m = 1$
and for each $j$ in $\text{unique\_prime\_factors}(m-1)$ {
$i^{(m - 1)/j} \space\text{ mod }\space m \ne 1$
}
then $i$ is generator $g$
}
3. Compute $n$'th primitive root with the generator, $g^{(m-1)/n} \space\text{ mod }\space m$
4. The rest of the NTT procedure is identical to the forward and inverse transforms previous described in the FFT section.

# Implementing the NTT.
## Pre-processing of NTT; creating the modulus type.
Essentially, NTT consists of two parts:
1. Construct a prime finite field from inputs with an underlying type unsigned whole number type (e.g. unsigned 64 bit) and determine the modulus. This means we avoid fractions completely while maintaining correctness.
2. Apply core FFT algorithm with the primitive $n$'th root of unity of the finite field domain, which can be calculated given the modulus.

So we must perform pre-processing of inputs to determine the modulus and the primitive root. The following are sample functions that show how to implement the pre-processing.
These are heavily influenced by [nayuki's excellent implementation and explanation](https://www.nayuki.io/page/number-theoretic-transform-integer-dft), please check that page out if you are confused, they have also provided an interactive tool so you can play with the numbers.

Computing the modulus:

// Returns the smallest modulus mod such that mod = i * len + 1
// for some integer i >= 1, mod > len, and mod is prime.
// Although the loop might run for a long time and create arbitrarily large numbers,
// Dirichlet's theorem guarantees that such a prime number must exist.
pub fn find_modulus(len: usize, minimum: i128) -> i128 {
    let vl = len as i128;
    let mut start = (minimum - 1 + vl - 1) / vl;
    if start < 1 {
        start = 1;
    }
    let mut n = start * vl + 1;
    loop {
        if is_prime(n) {
            return n;
        }
        n += vl;
    }
}

Computing the primitive root:

// Returns an arbitrary primitive degree-th root of unity modulo mod.
// totient must be a multiple of degree. If mod is prime, an answer must exist.
pub fn find_primitive_root(degree: i128, totient: i128, modulus: i128) -> i128 {
    let gen = find_generator(totient, modulus);
    pow_mod(gen, (totient / degree) as i128, modulus)
}

// Returns an arbitrary generator of the multiplicative group of integers modulo mod.
// totient must equal the Euler phi function of mod. If mod is prime, an answer must exist.
pub fn find_generator(totient: i128, modulus: i128) -> i128 {
    for i in 1..modulus {
        if is_primitive_root(i, totient, modulus) {
            return i;
        }
    }
}

pub fn is_primitive_root(val: i128, totient: i128, modulus: i128) -> bool {
    pow_mod(val, totient, modulus) == 1
        && unique_prime_factors(totient)
            .into_iter()
            .all(|p| pow_mod(val, totient / p, modulus) != 1)
}

// Returns a list of unique prime factors of the given integer in
// ascending order. For example, unique_prime_factors(60) = vec![2, 3, 5].
pub fn unique_prime_factors(mut n: i128) -> Vec<i128> {
    let mut result = Vec::new();
    for i in 2.. {
        if i * i > n {
            break;
        }
        if n % i == 0 {
            result.push(i);
            while n % i == 0 {
                n /= i;
            }
        }
    }
    if n > 1 {
        result.push(n);
    }
    result
}

## Core operations of Modulus type
Additionally, our new integer type must be a modulus integer type, that is for every operation: addition, subtraction, multiplication, and multiplicative inverse, must apply the modulus to the result using the modulus operator. For example, the core operations can be defined as follows:

fn add(lhs, rhs) { check_bounds(lhs + rhs) % modulo }
fn mul(lhs, rhs) { check_bounds(lhs * rhs) % modulo }
fn sub(lhs, rhs) { check_bounds(lhs - rhs) % modulo }

If any of the core operations exceed the bounds, we must either subtract or add the modulo. For example, if add or mul results exceed the modulo, then negate the modulo value; if the sub results exceed the modulo, then add the modulo. For example, the pseudocode would look like this:

fn check_bounds(input : F) -> F {
    if (input >= modulo) {
        return (input - modulo); 
    }
    if (input < 0) {
        return (input + modulo); 
    }
    return input;
}

## Core NTT implementation
This is the same as the "definitive FFT pseudocode", detailed on another page: [Tutorial 4 - Understanding the FFT](https://hackmd.io/@qozymandias/B1sv7DoFT).
However, we must instead use our Modulo integer type we have defined thus far, other than that the key difference is the primitive root has changed as we have now derived it from the finite field.

### Bringing it all together
The following example program applies the NTT to the inputs and the INTT to the result.
The program demonstrates how to glue the operations together; first step is the pre-processing in which the modulus and primitive root are found, followed by applying the core NTT transforms with the root found.

let inputs: Vec<i128> = vec![11, 42, 31, 43, -11, 12, 78, 37];
let mut max_val = 0;
for &x in inputs.iter() {
    if x > max_val {
        max_val = x;
    }
}

let min_mod = max_val * max_val * inputs.len() as i128 + 1;
let modulus = find_modulus(inputs.len(), min_mod);
let root = find_primitive_root(inputs.len() as i128, modulus - 1, modulus);

let temp0 = ntt_recursive(&inputs, root, modulus);
let res = intt_recursive(&temp0, root, modulus);

# NTT implementation for a predefined Finite Field
All of the above pre-processing can be rather tricky and inefficient for generic inputs, so instead it's much easier to use an existing [Finite Field library](https://docs.rs/ff/latest/ff/) and then map our inputs into this field.

## Pre-processing stage of NTT implementation
Below is the "pre-processing" stage, where we figure out the primitive root of the field. This is quite obscure because it doesn't follow anything I said so far; this is because it uses existing constants of the library which have been exposed for this very purpose.

The method to determine the primitive root for our inputs is rather strange, however, I will explain. The `F::root_of_unity` is given, however it's not "our" primitive root, the documentation is helpful but obscure:

/// Returns the `2^s` root of unity.
///
/// It can be calculated by exponentiating `Self::multiplicative_generator` by `t`,
/// where `t = (modulus - 1) >> Self::S`.
fn root_of_unity() -> Self;

So it's for `2^s`! Where `s` is `F::S` and this is another constant defined in the library.

/// An integer `s` satisfying the equation `2^s * t = modulus - 1` with `t` odd.
///
/// This is the number of leading zero bits in the little-endian bit representation of
/// `modulus - 1`.
const S: u32;

We can think of `2^(F::S)` as the "upper bound" of the finite field, and most importantly, we have its primitive root, `(2^(F::S))`'th primitive root of unity.

The halo2 library contains a detailed usage of these constants to derive the primitive root for our input polynomial length -- this can be found [here](https://github.com/DelphinusLab/halo2-gpu-specific/blob/a2c3fba52ab2a46e08a3f161f352cf0ecd43c08b/halo2_proofs/src/poly/domain.rs#L44) -- my version takes cues from it but is much simpler for learning purposes. Let's take a look at the comments to better understand the derivation.

        // Get omega, the 2^{k}'th root of unity (i.e. n'th root of unity)
        // The loop computes omega = extended_omega ^ {2 ^ (extended_k - k)}
        //           = (omega^{2 ^ (S - extended_k)})  ^ {2 ^ (extended_k - k)}
        //           = omega ^ {2 ^ (S - k)}.
        // Notice that omega ^ {2^k} = omega ^ {2^S} = 1.

The key information to gleam is that the $n$'th primitive root of unity is given by `omega ^ (2 ^ (S - k))`. This was very tricky for me to understand, so here is my interpretation of it: The primitive root is `omega`; `omega` is given by `root_of_unity`, but this is for `2^F::S`. We must move to `2^k` (where `k = log_2(d)` and `d` is the degree of the input polynomial -- or, more specifically, the required evaluation points we need in evaluation form). Because, we want to calculate the `2^{k}`'th root of unity.

Intuitively, we know roots can be "moved" by squaring the current root. Think back to the FFT, each "level" of recursion the length of the polynomial halved and its root is squared to ensure it's still a root of unity of the halved polynomial! Therefore, to "move" `omega` and make it the `2^k` 'th primitive root of unity, we have to square it for some number of iterations. The number of iterations would be the difference between `2^F::S` and `2^k`. Since we already have `F::S` we need to determine `k`. This can be computed like so: `k = log_2(degree)`. And this is also the number of "levels" of recursion, i.e. the number of times the polynomial is halved, which makes sense considering our FFT intuition. Now we simply keep squaring `omega` by itself as many times as the difference between `F::S` and `k`.

let k = (degree as f64).log(2.0) as u32;
for _ in 0..(F::S - k) {
    omega = omega.square();
}
let primitive_root = omega;

And we have got the primitive root! Much easier to compute that building everything from scratch as we did before. Now, if we continue squaring by `k` iterations, we which reach `F::S` and `omega` will become `1` because it's a primitive root of unity. This also means as the "bottom" level of recursion `omega` will have to be squared `k` times and will equal `1` which is an essential property for FFT! We assert this fact in the sample code.

pub fn new(degree: u32) -> Self {
    let mut omega = F::root_of_unity();
    // evaluation points we need to check.

    // level = log_2(degree)
    // g = omega = root of unity
    // Omega must stasify the equation: g^(2^final_level) = 1
    //
    // w_in = g^(2^level)
    // w = w_in^idx
    let k = (degree as f64).log(2.0) as u32;
    for _ in 0..(F::S - k) {
        omega = omega.square();
    }
    let primitive_root = omega;

    if cfg!(debug_assertions) {
        for _ in (F::S - k)..F::S {
            omega = omega.square();
        }
        assert_eq!(
            omega,
            F::from(1),
            "Omega must stasify the equation: g^(2^(log_2(degree)) = 1"
        );
    }

    NTT {
        primitive_root: primitive_root,
    }
}

## "Core FFT" stage of NTT implementation
Below is the "core FFT" stage, where we figure out the primitive root of the field.

fn forward_recursive(&self, invec: &Vec<F>, root: &F, n_calls: usize) -> Vec<F> {
    let n = invec.len();
    if n == 1 {
        return vec![invec[0].clone()];
    }

    let half_n = n / 2;

    // Separate even and odd coefficients
    let (even, odd): (Vec<_>, Vec<_>) =
        invec.iter().enumerate().partition(|(i, _)| *i % 2 == 0);
    let even: Vec<_> = even.into_iter().map(|(_, v)| v.clone()).collect();
    let odd: Vec<_> = odd.into_iter().map(|(_, v)| v.clone()).collect();

    let root_squared = root.clone() * root.clone();
    // Recursively compute NTT on even and odd parts
    let even_ntt = self.forward_recursive(&even, &root_squared, n_calls + 1);
    let odd_ntt = self.forward_recursive(&odd, &root_squared, n_calls + 1);

    // Combine the results
    let mut outvec = vec![F::from(0); n];
    let mut current_root = F::from(1);
    for i in 0..half_n {
        let bfy = current_root * odd_ntt[i];
        outvec[i] = even_ntt[i] + bfy;
        outvec[i + half_n] = even_ntt[i] - bfy;
        current_root *= root;
    }
    outvec
}

fn inverse_recursive(&self, invec: &Vec<F>, root: &F) -> Vec<F> {
    let mut outvec = self.forward_recursive(invec, &(root.invert().unwrap()), 0);
    let scaler = F::from(invec.len() as u64).invert().unwrap();

    for i in 0..outvec.len() {
        outvec[i] *= scaler.clone();
    }
    outvec
}

# Epilogue: NTT code samples
Learning by doing is always recommended, so here are my "toy" programs that demonstrate NTT implementation and usage.
- https://github.com/qozymandias/ntt-algorithms
- https://github.com/qozymandias/minty