# Finding a Primitive p-th Root of Unity in a Finite Field, with C++ Code

To this day, no method of finding a generator of $Z_p^*$ is known to be more efficient than essentially trying 2, then 3, and so on. Who cares? Well, the difficulty of breaking a certain public key cryptosystem (due to El Gamal) depends on the difficulty of working with generators of $Z_p^*$.Keith Conrad

An $n$th root of unity in a finite field $F$ is an element $r \in F$ satisfying $r^n=1$, where $n$ is an integer. If $n$ is the smallest positive integer with this property, $r$ is called a primitive $n$th root of unity. If $r$ is a primitive $n$th root of unity, then all elements in the set $\mu_n = \{1, r, r^2, \cdots, r^{n-1}\}$ are also roots of unity. Actually, the set $\mu_n$ form a cyclic group of order $n$ under multiplication, with generator $r$.

Problem: Suppose you are given a finite field $F=GF(2^d)$ of degree $d$, and you are promised that there indeed exists a primitive $p$th root of unity $r\in F$ for $p$ prime. Find $r$, and in particular, produce a C++code that finds it.

In what follows, we talk about how to find such a root and provide my C++ code; the code uses the awesome NTL library.

Let $F^*=F\setminus \{0\}$ be the multiplicative group of the finite field $F$. Let $s=2^d-1=ord(F^*)$ is the order of $F^*$. For a $p$th root $r$ of $F^*$ to exist, $p$ must divide $s$. Why? Because by Lagrange’s theorem, the order of a subgroup, here the group of the $p$th roots of unity under multiplication, must divide the order of the group, which is $s=ord(F^*)$.

So, how does one find $r$? The short answer is, via a generator $g$ of $F^*$ in three steps: (1) Randomly pick an element $g \in F, g \neq 0$. (2) check that $g^{s/q} \neq 1$ for all prime factors $q$ of $s$. (3) The element $r=g^{s/p}$ will be a $p$th root of unity in $F^*$.

The following C++ code is written using Victor Shoup’s super awesome library NTL. The code is part of another project I am working on: I need to compute the rank of submatrices of a Vandermonde matrix generated by a $p$th root of unity in the finite field $GF(2^d)$, given $p$ divides $2^d - 1$. The field $F=GF(2^d)$ is represented as the set of polynomials of degree strictly less than $d$ over $GF(2)$, which is isomorphic to the quotient ring $GF(2)[X] / (f)$ where $f$ is an irreducible polynomial of degree exactly $d$. Since $f$ is irreducible over $GF(2)$, this ring is indeed a field.

The following code implicitly uses namespace NTL. The type NTL::Vec has been typedef’d to Vector.

Part 1: Finding a generator of a cyclic group

The multiplicative group of a finite field is cyclic. (Why?) Moreover, it is a direct product of many cyclic subgroups.

There is no algebraic method to directly pinpoint a generator of a multiplicative group.  The number of generators of $F^*$ equals the number of integers in $1, 2, \cdots, s$ that are coprime to $s$. This quantity is $\phi(s)$, where $\phi(.)$ is Euler’s totient function defined as $\phi(p)=p-1, \phi(p^k)=p^k-p^{k-1},$ and $\phi(p^j \cdots q^k) = \phi(p^j) \cdots \phi(q^k)$. It is known that the ratio $\phi(n)/n$ is at least $n/6\log \log n$. (Source?) Hence a random element has a decent probability of being a generator.

Now, given an element $g\in F^*$, how do we test whether a given element $g$ generates the entire multiplicative group $F^*$ and not any subgroup? For $g$ to be a generator of $F^*$, we need $o_g = ord(g)=s$. Obviously, we need to check that $g^s=1$. Second, we need to check that for all positive integers $t$ between $2$ and $s-1$, that $g^t \neq 1$. However, there is a better of

However, there is a better way of checking this provided we are also given the list of prime factors of $s$. If $g$ generates a subgroup of $F^*$, then $o_g$ would be a proper divisor of $s$ since by Lagrange’s theorem, the order of a subgroup divides the order of the group. Since we have already seen $latx g^p=1$, we can be sure that $g^q \neq 1$ for all other prime factors of $s$. We don’t need to test all divisors of $s$; checking that $g^{s/q} \neq 1$ for all primes $q$ dividing $s$ suffices since they cover the candidate divisors.

The following code uses the function factorInt() to extract unique prime factors of $s$; this function is discussed later on.


GF2X findGenerator() {
GF2X gen;
ZZ otherFactor(groupOrder / params.p);
factors.kill();
factorInt(factors, otherFactor);
unique&lt;zz&gt;(factors);
factors.append(ZZ(params.p));

bool genGood = false;
while (not genGood) {
// try a random element as a generator
gen = GF2X::zero();
while (IsZero(gen))
random(gen, params.degree);
// gen != zero

genGood = true; // temporarily
for (auto factor : factors) {
ZZ orderOverFactor(groupOrder / factor);
if (unity == PowerMod(gen, orderOverFactor, modpoly)) {
genGood = false;
break;
}
}
}
// gen^factor != unity for all factors, including p
return gen;
}



Part 2: Finding a primitive root of unity

Next, let us find a primitive $p$th root of unity using the generator $g$.  If we set $\alpha = s/p$ and $r=g^\alpha$ where $g$ is a generator of $F^*$, then $r^p=g^{\alpha p} = g^s = 1$.  Note that  $s/p$ must be a positive integer since $p$ divides $s$. It follows that $r$ is a $p$th root of unity in the finite field $F$. Moreover, since $p$ is prime, the cyclic group $\mu_p$ has no non-trivial subgroups, meaning there cannot be any integer $n < p$ such that $r^n=1$. Hence $r$ is a primitive $p$th root of unity.

(To do: how would one proceed if $p$ was not a prime and we insisted on a primitive root instead of any root?)


// Set the unity and the monomial X element
SetCoeff(unity, 0, 1);
SetCoeff(X, 1, 1);

// Multiplicative group order
ZZ one(1L);
groupOrder = LeftShift(one, params.degree) - 1;

// unknown degree, got to build irred and the modulus
BuildIrred(irred, params.degree);
build(modpoly, irred);
GF2E::init(irred);

// p must divide the (multiplicative) group order
if (groupOrder == (ZZ) params.p)
{
// prime order
gen = rootP = X;
factors.kill();
factors.append((ZZ) params.p);
}
else {
gen = findGenerator();
// gen is a generator
// compute r=gen^(ord/p) so that r^p = gen^ord = 1
// since gen is a generator, gen^x = 1 for no x &lt; ord
rootP = PowerMod(gen, groupOrder / params.p, modpoly);
}



Part 3: Finding the unique prime factors of any integer

The above code uses a factorization function to factorize $s = ord(F^*)$. Since we don’t need the exponents of individual primes for our application, only a list of unique factors suffices. Here is a code which implements Pollard’s rho algorithm, combined with Miller’s probabilistic primality test. It’s not supersonic, but good enough.

Given an integer $n$, the algorithm basically finds one (possibly composite) factor and then proceeds recursively. First, it checks if $n$ is even (or a power of 2), in which case $2$ a factor. (We could have done this for first few primes, but I didn’t.)

Next, the obvious: it (probabilistically) checks whether $n$ is a prime, and if so, it gives up. The checking is performed using the celebrated Miller-Rabin probabilistic primality test. (Digression: I had a chanced to attend a lecture/conference from Gary Miller; he is awesome.)

Next, we randomly seed Pollard’s rho algorithm and find a factor. If it fails, we retry several times. If it still fails, we think that $n$ is probably prime. I tried 10 times, but you could do as many times as you please. The function $g(x)=x^2+1 \mod n$ is used to imitate a pseudo-random sequence.

Once Pollard’s rho algorithm finds a factor $m$ of $n$, the next is obvious: we remove all powers of $m$ from $n$. Then we recurse on $m$ and $n/m$.

void factorInt(Vec&lt;ZZ&gt;&amp; factors, ZZ n) {
if (n &lt;= 1)
return;
ZZ x, y, gcd, seed;
ZZ numTrials(10L);
ZZ trial(0L);
ZZ q;
//long n2 = conv&lt;long&gt;(n);

do {

// even
if (not IsOdd(n)) {
factors.append(ZZ(2));
while (n &gt; 2 &amp;&amp; divide(q, n, 2))
n = q;
// power of two
if (n &lt;= 2)
return;
}

if (ProbPrime(n))
break;

// n is not a prime
RandomBnd(seed, numTrials);
x = y = seed + 2;
gcd = ZZ(1L);

//n2 = conv&lt;long&gt;(n);

while (gcd == 1) {
g(x, x, n);
g(y, y, n);
g(y, y, n);
gcd = GCD( abs(x-y), n);
}
long gcd2 = conv&lt;long&gt;(gcd);
if (gcd == n) {
trial++;
continue; // retry
}
else {
// remove all powers of gcd from n
while (divide(q, n, gcd))
n = q;
//long n2 = conv&lt;long&gt;(n);

// subsequent factors
Vec&lt;ZZ&gt; childFactors;
if (n &gt; 1)
factorInt(childFactors, n);
if (gcd &gt; 1)
factorInt(childFactors, gcd);

factors.append(childFactors);

// done
return;
}
} while (n &gt; 1 &amp;&amp; trial &lt; numTrials );

// probably prime
factors.append(n);
return;

}


A brief remark on Pollard’s rho algorithm. (Please let me know if I have a gap in my argument.) Imagine the group $Z/nZ$ as a cycle of size $n$. Then, the function $g(x)$ samples a (sub)cycle $Z/mZ$ of this big cycle where $m$ divides $n$. If we are unlucky, the sampled cycle is the original one. Otherwise, the sampled cycle is pretty small (of size $m < n$). This means the set $H=\{g(x), g^2(x), \cdots, g^m(x)\}$, corresponding to the seed $x$, the element $g(x)$ being the identity element with order $m$. Hence it is a subgroup of $Z/nZ$. By Lagrange’s theorem, $m=ord(H)$ must divide $n=ord(Z/nZ)$.

Note that Pollard’s rho algorithm actually detects two elements $x^\prime, y^\prime \in H$ which are congruent modulo $ord(H)$. Then the difference $d = \vert x^\prime - y^\prime \vert$ must be a multiple of $ord(H)$.  The algorithm outputs $gcd(d, n)$, which is strictly less than $n$ if indeed $H \subset G$.

As mentioned in the Wikipedia article linked above, this algorithm can be improved by replacing multiple GCD computations with a single computation (in batch, you could say).

Bonus material: Unique elements from an NTL vector

Finally, the following code keeps only the unique elements in an NTL vector. I used a brute force  $O(n^2)$ time loop. Ideally, one should use a better search method (possibly via a binary search tree or a hash table.) However, GCC was complaining when I wanted to use a hashmap with NTL::ZZ keys/values, so I gave up for the sake of some development-time.

template &lt;typename T&gt;
void unique(Vec&lt;T&gt; &amp; src) {
std::vector&lt;T&gt; arr;
toStdVector&lt;T&gt;(arr, src);
src.kill();
while (arr.size() != 0) {
T elem = arr[0];
src.append(elem);

auto it = arr.begin();
while (it != arr.end()) {
if (*it == elem)
it = arr.erase(it);
else
it++;
}
}
//
//
//  typedef std::unordered_map&lt;T, T&gt; Map;
//  Map ht;
//  for (T&amp; item : src) {
//    if (ht.find(item) == ht.end())
//      ht[item] = item;
//  }
//  src.kill();
//  for (auto pos = ht.begin(); pos != ht.end(); pos++) {
//    src.append(pos-&gt;first);
//  }
}



I’d be glad to hear what you think. How would you have done it?