Discrete logarithms modulo a prime (Pohlig–Hellman)

Let p be prime and let g be a primitive root modulo p, so that every y \in (\mathbb{Z}/p\mathbb{Z})^\times is g^x for a unique x \in [0, p - 1), the discrete logarithm of y. For a general p no fast algorithm is known. When the group order p - 1 factors into small primes, the method of Pohlig and Hellman [PohHel1978] reduces the problem to logarithms in groups of prime order, which are solved by Shanks' baby-step giant-step method [Sha1971]. See also [CraPom2005].

FLINT implements this method in the nmod_discrete_log_pohlig_hellman functions.

Reduction to prime-power order#

Write p - 1 = \prod_{i=1}^{k} q_i^{e_i} with distinct primes q_i, and let c_i = (p-1)/q_i^{e_i}. The element h_i = g^{c_i} has order exactly q_i^{e_i}, and y^{c_i} = h_i^{x}. So x_i = x \bmod q_i^{e_i} is the logarithm of y^{c_i} to the base h_i, a problem in a group of order q_i^{e_i}. The x_i determine x by the Chinese remainder theorem:

where the idempotent \epsilon_i = c_i \cdot (c_i^{-1} \bmod q_i^{e_i}) satisfies \epsilon_i \equiv 1 \pmod{q_i^{e_i}} and \epsilon_i \equiv 0 \pmod{q_j^{e_j}} for j \ne i.

Digits in a prime-power group#

Fix i and write q = q_i, e = e_i, h = h_i and x_i = \sum_{j=0}^{e-1} d_j q^j with 0 \le d_j < q. The element \gamma = h^{q^{e-1}} = g^{(p-1)/q} has order q.

Logarithms of prime order#

Step 2.2 is a discrete logarithm in the group of order q generated by \gamma. FLINT uses a baby-step giant-step search with a precomputed table of bounded size.

For q = 2 no search is needed: \gamma = -1, and d_j = 0 if w = 1, d_j = 1 if w = \gamma.

For odd q, choose a stride m and let C = \lceil q / m \rceil. The table holds the pairs (\gamma^{cm}, cm) for 0 \le c < C, sorted by the first entry. Any d \in [0, q) can be written d = cm + d' with 0 \le c < C and 0 \le d' < m, so \gamma^{d} = w if and only if w \gamma^{-d'} = \gamma^{cm}.

The search returns the correct digit. Let d \in [0, q) be the true logarithm of the input w. For d' \le d, the current value is \gamma^{d - d'} with 0 \le d - d' < q. Since 0 \le cm < q for every table entry, it matches an entry exactly when d - d' = cm. The first such d' is d \bmod m, with c = \lfloor d/m \rfloor < C.

The table costs C entries and the search at most m multiplications and m binary searches. Balancing the two gives m \approx \sqrt q. FLINT starts from m = \lceil \sqrt q \rceil and doubles m until C \le 100. So the table has at most 100 entries, and the search costs O(\sqrt q) multiplications for q up to about 10^4, but between q/100 and q/50 for larger q. This keeps the precomputation small and favours primes p for which p - 1 is smooth, which is how the functions are used in FLINT.

Precomputation#

nmod_discrete_log_pohlig_hellman_precompute_prime performs the parts that depend only on p:

1. Factor p - 1 with n_factor.
2. For each q_i, compute c_i, \epsilon_i \bmod (p - 1) and q_i^{e_i - 1}.
3. Take g as the least a \ge 1 with a^{(p-1)/q_i} \ne 1 for all i, that is, the least primitive root. If p is prime, one exists. If none is found below p, an exception is thrown.
4. For each i, compute \gamma_i = g^{(p-1)/q_i}, \gamma_i^{-1}, h_i^{-1} = g^{-c_i}, the stride m_i and the sorted table.

It returns an estimate of the number of multiplications of one run: for each i, the cost of the exponentiation y^{c_i}, plus for each of the e_i digits the cost of the exponentiation in step 2.1 of Algorithm 1, m_i (1 + \ln C_i) for the search and 2 \ln q_i for the power of \beta.

The full algorithm#

The sum in step 2 is accumulated in three words and reduced once at the end. If a search in Algorithm 2 fails, or w \notin \{1, -1\} when q = 2, an exception is thrown. For an odd prime p this happens only for y = 0, which has no logarithm.

Cost. Each digit needs an exponentiation by q^{e-1-j} < p, a search and an exponentiation of \beta, so the total is O(\sum_i e_i (\log p + m_i)) multiplications, plus one exponentiation y^{c_i} for each i. When all q_i are small, this is O(\log^2 p).

Use in FLINT#

• The sparse interpolation in the multivariate gcd code of fmpz_mpoly_factor evaluates at powers of a primitive root and recovers exponents with nmod_discrete_log_pohlig_hellman_run. It chooses primes p with fmpz_next_smooth_prime, for which p - 1 currently has no prime factor above 23, so the searches are short.
• fmpz_mod_discrete_log_pohlig_hellman_run in fmpz_mod is the same method for multiprecision primes.
• n_discrete_log_bsgs in ulong_extras is a plain baby-step giant-step search, without the reduction to prime order.
• The dlog module, used for Dirichlet characters, precomputes tables for many logarithms to the same base modulo a prime power (Discrete logarithm tables).