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

242_roots_unity_9
9th roots of unity. Source: Kyle Miller’s article.

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 nth 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 nth root of unity. If r is a primitive nth 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 pth 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 pth 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 pth 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 pth root of unity in F^*.

About the code

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 pth 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<zz>(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 pth 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 pth 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 pth 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?

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s