Finding Primitive Roots of Unity
Primitive $n$th roots of unity are field elements $\alpha \in \mathbb{F}_p$ satisfying
\[\alpha^n = 1 \text{ and } \alpha^k \neq 1 \text{ for } 0 < k < n\]These are useful for speeding up polynomial multiplication using the number theoretic transform. Most descriptions of the NTT assume the existence of a primitive $n$th root of unity and we often work over fields where they exist by construction. This post describes how to find them.
Existence
$\mathbb{F}_p$ has a primitive $n$th root of unity iff $n  p1$.

If $\alpha$ is a primitive $n$th root of unity, it generates a subgroup of order $n$ (by definition). This order must divide the order of the original group $\phi(p) = p1$.

If $n  (p1)$ then for any generator $g \in \mathbb{F}_p$, $\alpha = g^{(p1)/n}$ is a primitive $n$th root of unity since $\alpha^{n} = g^{p1} = 1$ and $\alpha^{k} = g^{\frac{k(p1)}{n}} \neq 1$ for $0 < k < n$.
Moreover, there will be $\phi(n)$ of them. Since if $\alpha$ is a primitive $n$th root of unity then $\alpha^k$ with $0 < k \leq n$ and $\gcd(n, k) = 1$ is another.
Special Primes
In practice, we assume polynomial degrees are a power of $2$. This leads to primes of the form $p = k \cdot 2^m+1$ so that $2^m  (p1)$ and $\mathbb{F}_p$ has a primitive $2^m$th root of unity.
For instance, many RLWE schemes use polynomial degree $n = 1024$ and set $p = 5 \cdot 2^{13} + 1$.
This gives us a simple algorithm for finding a primitive $2^m$th root of unity:

Sample a nonzero element $x \leftarrow \mathbb{F}_p$

If $x^{(p1)/2} = 1$ go back to step 1.
$\alpha = x^{\frac{(p1)}{2^m}}$ is a primitive $2^m$th root of unity since $\alpha^{2^m} = x^{p1} = 1$ and $x^{\frac{(p1)k}{2^m}} \neq 1 $ for $0 < k < 2^m$.
More generally, if $n = \prod_{i} p_i^{e_i}$ is the prime factorization of $n$, $\alpha = \prod_{i} x_i^{(p1) / p_i^{e_i}}$ is a primitive $n$th root of unity where $x_i^{(p1)/p_i} \neq 1$.
p = 0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001
K = GF(p)
def find_proot(n):
a = K(1)
for (i, j) in factor(n):
x = K(1)
while pow(x, (p1)/i, p) == 1:
x = K.random_element()
a *= pow(x, (p1)/(i^j), p)
return a
factors = [i^j for (i,j) in factor(p1)]
n = factors[0]
a = find_proot(n)
# a is an nth root of unity mod p
assert(a != 1 and a^n == 1)