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).