Advanced Number Theory, in depth: group order, primitive roots, baby-step giant-step and Pohlig-Hellman discrete logs, built on a factoring pipeline

The multiplicative group modulo a prime

Fix a prime p. The nonzero residues 1, 2, ..., p - 1 form a group under multiplication modulo p. It is closed, every element has an inverse, and it has exactly p - 1 elements. Write it as (Z/pZ)*.

The order of an element g is the smallest t > 0 with g^t = 1 (mod p). Lagrange's theorem says the order of every element divides the size of the group, p - 1. The powers of g cycle with period equal to its order, so g generates a subgroup of exactly that size.

A primitive root is an element whose order is the full p - 1, so its powers hit every nonzero residue. Every prime has one, and in fact phi(p - 1) of them. Once you have a primitive root g, every residue h can be written as h = g^x for exactly one x in [0, p - 1). That x is the discrete logarithm of h to base g. Computing g^x from x is fast. Computing x from g^x is the hard problem that Diffie-Hellman and DSA-style signatures rely on.

Here is the key structural fact. Because orders divide p - 1, the structure of the group is set by the prime factorisation of p - 1. Every algorithm below uses that factorisation, so we need it first.

The prerequisite: factoring p - 1

Finding the order of an element or testing for a primitive root needs the distinct prime factors of p - 1. For numbers up to about 10¹⁸, the standard pipeline has three stages. Trial division strips small primes. A deterministic Miller-Rabin test recognises when what's left is prime. Pollard's rho with Brent's cycle detection splits composites in about n^1/4 steps, using the rho-shaped sequence explained in the cycle detection guide. The code below is the compact version. Those two pages explain why each part works.

import math, random

def is_probable_prime(n):
    if n < 2:
        return False
    small = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)
    for q in small:
        if n % q == 0:
            return n == q
    d, s = n - 1, 0
    while d % 2 == 0:
        d //= 2
        s += 1
    for a in small:              # these 12 bases are deterministic for n < 3.18e24
        x = pow(a, d, n)
        if x in (1, n - 1):
            continue
        for _ in range(s - 1):
            x = x * x % n
            if x == n - 1:
                break
        else:
            return False
    return True

def pollard_brent(n):
    if n % 2 == 0:
        return 2
    while True:
        y, c, m = random.randrange(1, n), random.randrange(1, n), 128
        g = r = q = 1
        while g == 1:
            x = y
            for _ in range(r):
                y = (y * y + c) % n
            k = 0
            while k < r and g == 1:
                ys = y
                for _ in range(min(m, r - k)):
                    y = (y * y + c) % n
                    q = q * abs(x - y) % n
                g = math.gcd(q, n)
                k += m
            r *= 2
        if g == n:                 # batched gcd overshot: replay one step at a time
            g = 1
            while g == 1:
                ys = (ys * ys + c) % n
                g = math.gcd(abs(x - ys), n)
        if g != n:
            return g

def factorize(n, out=None):
    out = {} if out is None else out
    for q in (2, 3, 5, 7, 11, 13):
        while n % q == 0:
            out[q] = out.get(q, 0) + 1
            n //= q
    if n == 1:
        return out
    if is_probable_prime(n):
        out[n] = out.get(n, 0) + 1
        return out
    d = pollard_brent(n)
    factorize(d, out)
    factorize(n // d, out)
    return out

It splits 2⁶⁴ + 1 into 274177 * 67280421310721 instantly. For a random 2048-bit p, the large prime factor of p - 1 is out of reach, and that is what makes such groups safe.

Element order and primitive roots

Given the factorisation p - 1 = q₁^e1 ... q_k^ek, you can compute an element's order without enumerating powers. Start with t = p - 1. For each prime q, keep dividing t by q while g^t/q is still 1. Whatever survives is the order. An element is a primitive root exactly when g^(p-1)/q differs from 1 for every prime q dividing p - 1. Primitive roots are common (the fraction is phi(p - 1)/(p - 1)), so trying 2, 3, 4 and so on finds one quickly.

def order(g, p, fac):
    """Multiplicative order of g mod prime p; fac is factorize(p - 1)."""
    t = p - 1
    for q, e in fac.items():
        for _ in range(e):
            if pow(g, t // q, p) == 1:
                t //= q
            else:
                break
    return t

def primitive_root(p, fac):
    for g in range(2, p):
        if all(pow(g, (p - 1) // q, p) != 1 for q in fac):
            return g

For p = 8101, p - 1 = 8100 = 2² * 3⁴ * 5². The element 2 has order 100, 3 has order 810, 4 has order 50, and 6 is the smallest primitive root. Its order is the full 8100.

Baby-step giant-step

The naive way to solve g^x = h is to try every x, which takes O(n) steps for a group of order n. Baby-step giant-step (Shanks) trades memory for time. Let m = ceil(sqrt(n)) and write the unknown as x = i*m + j with 0 <= i, j < m. Then g^x = h becomes g^j = h * (g^-m)ⁱ. Store every baby step g^j in a hash table, then walk the giant steps h, h*g^-m, h*g^-2m, ... until one lands in the table. That costs O(sqrt n) time and O(sqrt n) memory.

def bsgs(g, h, p, n):
    """Smallest x in [0, n) with g^x = h mod p; n must be >= the order of g."""
    m = math.isqrt(n - 1) + 1
    table, e = {}, 1
    for j in range(m):               # baby steps: g^j -> j (keep the smallest j)
        table.setdefault(e, j)
        e = e * g % p
    step = pow(g, -m, p)             # g^(-m); negative exponents need Python 3.8+
    gamma = h
    for i in range(m):               # giant steps: h * g^(-i*m)
        if gamma in table:
            return i * m + table[gamma]
        gamma = gamma * step % p
    return None                      # h is not a power of g

The square root is the important part. A 64-bit group needs about 2³² table entries, which is tens of gigabytes. That's feasible but painful. A 256-bit group needs 2¹²⁸, which is impossible. When memory is the limit, Pollard's rho for logarithms gets the same O(sqrt n) time in constant memory, using a pseudo-random walk and cycle detection.

Pohlig-Hellman: divide by the factorisation

Baby-step giant-step costs the square root of the group order. Pohlig-Hellman shows it really costs the square root of the largest prime factor of the order. It works in two steps.

1. Project into each prime-power part. For each q^e dividing n = p - 1, raising both sides to n/q^e moves the problem into the subgroup of order q^e. A solution there is x mod q^e.
2. Peel off one base-q digit at a time. Write x mod q^e = d₀ + d₁q + ... with each digit in [0, q). Raising to n/q^k+1 isolates digit d_k as a logarithm in a subgroup of order just q. That subproblem is a baby-step giant-step of size sqrt(q).

Then the Chinese remainder theorem combines the residues into x mod (p - 1). The total cost is about the sum over prime factors of e * sqrt(q) group operations, plus the factoring.

def crt(residues, moduli):
    x, M = 0, 1
    for r, m in zip(residues, moduli):
        t = (r - x) * pow(M, -1, m) % m      # moduli must be pairwise coprime
        x += M * t
        M *= m
    return x % M

def pohlig_hellman(g, h, p, fac):
    """x with g^x = h mod p; g a primitive root; fac = factorize(p - 1)."""
    n = p - 1
    residues, moduli = [], []
    for q, e in fac.items():
        gq = pow(g, n // q, p)                  # an element of order exactly q
        x = 0
        for k in range(e):
            hk = pow(h * pow(g, -x, p) % p, n // q ** (k + 1), p)
            d = bsgs(gq, hk, p, q)
            if d is None:
                raise ValueError("h is not in the subgroup generated by g")
            x += d * q ** k
        residues.append(x)
        moduli.append(q ** e)
    return crt(residues, moduli)

Worked example: a discrete log mod 8101

Solve 6^x = 7531 (mod 8101). The group order is 8100 = 4 * 81 * 25, so the problem splits into three subproblems.

• Modulo 4. The base-2 digits come out as 1 and then 0, so x = 1 (mod 4).
• Modulo 81. The base-3 digits come out as 2, 0, 2, 1, so x = 2 + 0*3 + 2*9 + 1*27 = 47 (mod 81).
• Modulo 25. The base-5 digits come out as 4 and then 2, so x = 4 + 2*5 = 14 (mod 25).

The CRT combines 1 mod 4, 47 mod 81 and 14 mod 25 into x = 6689. Check it: pow(6, 6689, 8101) returns 7531. Plain BSGS over the whole group would have needed a 90-entry table and up to 90 giant steps. Pohlig-Hellman never searched a set larger than five. The ratio grows without bound as p - 1 gets smoother. Our tests ran the full pipeline on 300 random primes between 10⁶ and 10⁹ with random exponents, and every logarithm came back correct.

What this means for cryptographic parameters

The worked example points straight at the defence. If p - 1 has only small prime factors, discrete logs are easy, no matter how big p is. Real systems therefore use one of two shapes:

• Safe primes, p = 2q + 1 with q prime. Then p - 1 = 2q, and the large factor q bounds Pohlig-Hellman. The finite-field Diffie-Hellman groups standardised in RFC 7919 are safe primes.
• Prime-order subgroups. Work in a subgroup of large prime order q inside (Z/pZ)*, or on an elliptic curve whose group order is a large prime times a small cofactor.

Generic algorithms such as BSGS and Pollard's rho cost sqrt(q). Over finite fields, index-calculus methods do better: they run in subexponential time by exploiting the fact that integers factor into small primes. That's why finite-field groups need 2048-bit or larger primes, while elliptic curves without that structure reach similar security at around 256 bits. RSA leans on the factoring half of this story. The extended Euclidean algorithm supplies its private exponent.

Failure modes

• Small-subgroup attacks. A Diffie-Hellman peer that accepts any public value can be handed an element of small order, which leaks the secret modulo that order. That's Pohlig-Hellman used as a weapon. Validate that received values lie in the intended prime-order subgroup.
• Base that is not a generator. If g is not a primitive root, h may have no logarithm at all, and a solution is only unique modulo the order of g. Check with order() first, or handle a None result.
• Overflow in fixed-width languages. In C or Rust, x * y % p overflows once p passes 2³². Use 128-bit intermediates or Montgomery multiplication.

Trade-offs

What to do next

1. Copy the code into one file and test it: pick random primes, random exponents, and check that pohlig_hellman inverts pow.
2. Time BSGS against Pohlig-Hellman for a 40-bit safe prime and a 40-bit prime with smooth p - 1, and explain the gap.
3. Implement Pollard's rho for logarithms and compare its memory with BSGS.
4. Read RFC 7919 and verify for yourself that one of its groups is a safe prime, using is_probable_prime on p and (p - 1)/2.
5. Audit any Diffie-Hellman code you own for public-value validation against small subgroups.