Primitive Roots and Discrete Logarithms, in depth: orders, the generator test, baby-step giant-step, Pohlig-Hellman, Pollard rho and safe parameters

Orders and primitive roots

Work in the multiplicative group modulo a prime p: the numbers 1 to p-1 under multiplication mod p. It has p-1 elements. The order of an element g is the smallest positive k with g^k ≡ 1. Fermat's little theorem says g^p-1 ≡ 1 for every g, and Lagrange's theorem says the order of any element divides the size of the group, so the order of g always divides p-1.

A primitive root (also called a generator) is an element whose order is exactly p-1. Its powers g⁰, g¹, ..., g^p-2 are all different and so cover every nonzero residue exactly once. Every prime has primitive roots, and there are exactly φ(p-1) of them, where φ is Euler's totient. For composite moduli the picture is narrower: a primitive root modulo n exists only when n is 1, 2, 4, p^k or 2p^k for an odd prime p. Modulo 8, for instance, every odd number squares to 1, so no element has order 4 and none generates the group.

Once you have a primitive root g, every nonzero h has a unique exponent x in [0, p-2] with g^x ≡ h. That x is the discrete logarithm of h to base g, written log_g h. Like an ordinary logarithm it turns multiplication into addition, except that exponents live modulo p-1: log(ab) = log a + log b (mod p-1).

Testing a candidate generator

Checking all p-1 powers is hopeless for large p. The order of g divides p-1, so if g is not a primitive root its order is a proper divisor d of p-1, and d divides (p-1)/q for at least one prime q dividing p-1. That gives a test that needs one modular exponentiation per distinct prime factor:

g is a primitive root mod p if and only if g^(p-1)/q ≠ 1 for every prime q dividing p-1.

The catch is that you must know the factorization of p-1. For random large primes that can be as hard as the problem you are trying to make hard, which is why cryptographic primes are generated with p-1 of known shape, most often a safe prime p = 2q + 1 with q also prime.

def prime_factors(n):
    """Distinct prime factors by trial division; fine for teaching-sized n."""
    fs, d = [], 2
    while d * d <= n:
        if n % d == 0:
            fs.append(d)
            while n % d == 0:
                n //= d
        d += 1
    if n > 1:
        fs.append(n)
    return fs

def is_primitive_root(g, p):
    n = p - 1
    return all(pow(g, n // q, p) != 1 for q in prime_factors(n))

def smallest_primitive_root(p):
    return next(g for g in range(2, p) if is_primitive_root(g, p))

assert [g for g in range(2, 23) if is_primitive_root(g, 23)] == [5, 7, 10, 11, 14, 15, 17, 19, 20, 21]
assert smallest_primitive_root(1000003) == 2

Primitive roots are not rare: a random g succeeds with probability φ(p-1)/(p-1), which is never tiny in practice, so a short search finds one. For a safe prime the test shrinks to two checks, g² ≠ 1 and g^q ≠ 1, which means g is a primitive root exactly when g is not ±1 and is a quadratic non-residue.

Worked example: modulo 23

Take p = 23, so p-1 = 22 = 2 × 11. Test g = 2: 2¹¹ = 2048 = 89 × 23 + 1, so 2¹¹ ≡ 1 and 2 has order 11. It is not a primitive root. Test g = 5: 5² = 25 ≡ 2, which is not 1, and 5¹¹ ≡ 22 ≡ -1, also not 1. Both checks pass, so 5 is a primitive root. The full list is 5, 7, 10, 11, 14, 15, 17, 19, 20 and 21: ten elements, matching φ(22) = 10.

Reading the cycle gives a log table: log₅ 2 = 2, log₅ 4 = 4, log₅ 8 = 6. Check the addition rule: 2 × 4 = 8 and 2 + 4 = 6. For a 2048-bit p such a table is impossible, which is why the algorithms below exist.

Baby-step giant-step

Baby-step giant-step, due to Shanks, trades memory for time. Let n be the order of g and m = ⌈√n⌉. Any exponent x below n can be written x = i·m + j with 0 ≤ i, j < m. Then g^x = h becomes g^j = h · (g^-m)ⁱ. Store the m baby steps g^j in a hash table, then walk the giant steps h·g^-im until one lands in the table. Both phases take m steps, so time and memory are O(√n).

from math import isqrt

def bsgs(g, h, p, n=None):
    """Return x in [0, n) with pow(g, x, p) == h, or None. n = order of g (default p - 1)."""
    n = n or p - 1
    m = isqrt(n - 1) + 1                 # ceil(sqrt(n)) without floating point
    table = {}
    e = 1
    for j in range(m):                   # baby steps: g^j -> j
        table.setdefault(e, j)
        e = e * g % p
    factor = pow(g, -m, p)               # g^(-m), Python 3.8+ modular inverse
    y = h % p
    for i in range(m):                   # giant steps: h * g^(-i*m)
        if y in table:
            return i * m + table[y]
        y = y * factor % p
    return None                          # h is not in the subgroup generated by g

assert bsgs(5, 8, 23) == 6

Trace it for 5^x ≡ 8 (mod 23). n = 22, so m = 5. Baby steps for j = 0..4 are 1, 5, 2, 10, 4. 5⁵ ≡ 20, and its inverse is 15 because 20 × 15 = 300 = 13 × 23 + 1. Giant step i = 0 tests 8: not in the table. Step i = 1 tests 8 × 15 = 120 ≡ 5: in the table with j = 1. So x = 1 × 5 + 1 = 6, matching the cycle.

Compute m with an integer square root; a float sqrt on a large n can round down and leave the last exponents unreachable.

Pohlig-Hellman: divide by the factors of the order

Baby-step giant-step costs the square root of the group order. Pohlig-Hellman shows that what really counts is the square root of the largest prime factor of the order. If n = p-1 factors as a product of prime powers q^e, solve the log modulo each q^e inside a subgroup of order q, then glue the answers with the Chinese remainder theorem.

The projection is one exponentiation. Raising both sides of g^x = h to the power n/q gives (g^n/q)^x = h^n/q, where g^n/q has order q, so this equation only sees x mod q. Higher powers q^e are handled one base-q digit at a time.

def factorize(n):
    out, d = {}, 2
    while d * d <= n:
        while n % d == 0:
            out[d] = out.get(d, 0) + 1
            n //= d
        d += 1
    if n > 1:
        out[n] = out.get(n, 0) + 1
    return out

def crt(residues, moduli):
    x, m = 0, 1
    for r, mi in zip(residues, moduli):
        t = (r - x) * pow(m, -1, mi) % mi
        x, m = x + m * t, m * mi
    return x % m

def pohlig_hellman(g, h, p):
    """g must be a primitive root mod p; cost is dominated by sqrt(largest prime factor of p-1)."""
    n = p - 1
    residues, moduli = [], []
    for q, e in factorize(n).items():
        gq = pow(g, n // q, p)                     # element of order q
        x_q = 0
        for k in range(e):                         # recover x mod q^e digit by digit
            hk = pow(h * pow(g, -x_q, p) % p, n // q ** (k + 1), p)
            x_q += bsgs(gq, hk, p, q) * q ** k
        residues.append(x_q)
        moduli.append(q ** e)
    return crt(residues, moduli)

assert pohlig_hellman(5, 8, 23) == 6
p = 1000003                                         # p - 1 = 2 * 3 * 166667
for x in (0, 1, 12345, 999999):
    assert pohlig_hellman(2, pow(2, x, p), p) == x

Trace the same instance. For q = 2: 5¹¹ ≡ -1 has order 2, and 8¹¹ ≡ 1, so x ≡ 0 (mod 2). That step is just a quadratic-residue test. For q = 11: 5² = 2 has order 11 and 8² = 64 ≡ 18; the powers of 2 run 1, 2, 4, 8, 16, 9, 18, so x ≡ 6 (mod 11). The Chinese remainder theorem combines x ≡ 0 (mod 2) and x ≡ 6 (mod 11) into x = 6.

The lesson for parameter choice is direct: if p-1 is a product of small primes, the discrete log is easy no matter how large p is. A 2048-bit prime whose p-1 has no factor above 2⁴⁰ falls in about 2²⁰ steps per factor.

Pollard rho and index calculus

Pollard's rho matches baby-step giant-step's O(√n) time with constant memory. It walks a pseudo-random sequence of elements of the form g^ah^b, tracking a and b, and uses cycle detection to find two points that collide. A collision g^ah^b = g^a'h^b' gives a linear equation for x modulo the order. Expected work is about √(πn/2) group operations.

Index calculus is the reason finite-field and elliptic-curve parameters differ so much in size. In the integers modulo p, elements can be factored over a base of small primes; collect enough relations, solve a linear system for the logs of the base, and any target's log follows. The number field sieve refines this to subexponential cost. No comparable attack is known for well-chosen elliptic curves, so 256-bit curves give roughly the security of 3072-bit prime fields.

How cryptography uses this

Real protocols rarely work in the full group generated by a primitive root. With a primitive root, the public value g^x leaks x mod 2 (it is a quadratic residue exactly when x is even), and more small factors of p-1 leak more bits. Instead, protocols use a subgroup of large prime order q: with a safe prime p = 2q + 1 that is the subgroup of quadratic residues, and implementations check that a received value y satisfies 1 < y < p-1 and y^q ≡ 1. Skipping that check permits small-subgroup attacks, where an attacker sends an element of small order and learns the victim's secret modulo that order, one Pohlig-Hellman piece at a time.

Prefer standardised groups such as those in RFC 7919, or curves like X25519, but at modern sizes: Logjam showed that a precomputation against one shared 512-bit prime breaks every server using it.

Failure modes

What to do next

1. Run the code above and print the log table for p = 23 and p = 29; confirm there are φ(p-1) primitive roots for each.
2. Time BSGS against Pohlig-Hellman on primes whose p-1 is smooth versus safe primes of the same size.
3. Implement Pollard's rho with Floyd cycle detection and check its answers against BSGS.
4. Review how exponentiation stays fast and constant-time in modular exponentiation, and how inverses are computed in the extended Euclidean algorithm.
5. See the protocol built on this problem in Diffie-Hellman key exchange, and compare its hardness assumption with factoring in RSA.