Pohlig-Hellman algorithm in Python
Reduces a discrete logarithm in a group of order N = prod l^e to logarithms in subgroups of prime order l and joins them by the CRT. The cost is set by the largest prime factor of N.
- File
pohlig_hellman.py, 78 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/pohlig-hellman/. The authors are listed in prov/PW5NFF.jsonl.
from math import isqrt
def _bsgs(g: int, h: int, n: int, p: int) -> int | None:
"""Least x >= 0 with g^x = h (mod p), or None; requires ord(g) <= n."""
m = isqrt(n - 1) + 1
table: dict[int, int] = {}
e = 1 % p
for j in range(m):
if e not in table:
table[e] = j
e = e * g % p
u = pow(g, -m, p)
gamma = h
for i in range(m):
if gamma in table:
return i * m + table[gamma]
gamma = gamma * u % p
return None
def _crt(residues: list[tuple[int, int]]) -> int:
"""The x in [0, prod m) with x = a (mod m) for pairwise coprime moduli m."""
x, mod = 0, 1
for a, m in residues:
t = (a - x) * pow(mod, -1, m) % m
x, mod = x + mod * t, mod * m
return x
def pohlig_hellman(g: int, h: int, p: int, factors: dict[int, int]) -> int | None:
"""Pohlig-Hellman discrete logarithm in the group of units modulo p.
factors maps each prime l to its exponent e in an integer N = prod l^e that
is a multiple of ord(g), for example the factorization of p - 1 for a prime
p. Return the least x >= 0 with g^x = h (mod p), or None if h is not a
power of g. The group need not be cyclic: a final check g^x = h decides.
Precondition: p >= 2, gcd(g, p) = 1, every key of factors is prime and
ord(g) divides N. Raises ValueError if p < 2.
Cost: O(log(N)^2 + sum of e (l)^(1/2)) multiplications modulo p, and
O(max l^(1/2)) stored residues.
"""
if p < 2:
raise ValueError("need p >= 2")
g %= p
h %= p
r = 1
for l, e in factors.items():
r *= l ** e
for l in factors: # @line L1
while r % l == 0 and pow(g, r // l, p) == 1 % p: # @line L2
r //= l # @line L3
if pow(h, r, p) != 1 % p: # @line L4
return None
residues: list[tuple[int, int]] = [] # @line L5
for l in factors: # @line L6
f, q = 0, r # @line L7
while q % l == 0:
f, q = f + 1, q // l
if f == 0:
continue
gl, hl = pow(g, q, p), pow(h, q, p) # @line L8
gamma = pow(gl, l ** (f - 1), p) # @line L9
xl = 0 # @line L10
for k in range(f): # @line L11
delta = pow(hl * pow(gl, -xl, p) % p, l ** (f - 1 - k), p) # @line L12
d = _bsgs(gamma, delta, l, p) # @line L13
if d is None: # @line L14
return None
xl += d * l ** k # @line L15
residues.append((xl, l ** f)) # @line L16
x = _crt(residues) # @line L17
if pow(g, x, p) != h: # @line L18
return None
return x # @line L19
test_pohlig_hellman.py 117 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/pohlig-hellman/. The authors are listed in prov/PW5NFF.jsonl.
import random
from math import gcd
import pytest
from pohlig_hellman import pohlig_hellman
def factor(n: int) -> dict[int, int]:
out: dict[int, int] = {}
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 is_prime(n: int) -> bool:
return n >= 2 and factor(n) == {n: 1}
def phi(n: int) -> int:
r = n
for q in factor(n):
r = r // q * (q - 1)
return r
def power_table(g: int, p: int) -> dict[int, int]:
seen: dict[int, int] = {}
e, x = 1 % p, 0
while e not in seen:
seen[e] = x
e, x = e * g % p, x + 1
return seen
def test_worked_example() -> None:
# Galbraith, Example 13.2.3: p = 19, g = 2, h = 5, p - 1 = 2 * 3^2, answer 16.
assert pohlig_hellman(2, 5, 19, {2: 1, 3: 2}) == 16
# HAC, Example 3.64: p = 251, alpha = 71, beta = 210, answer 197.
assert pohlig_hellman(71, 210, 251, factor(250)) == 197
# g = 4 has order 9 modulo 19, and 2 is not a power of 4.
assert pohlig_hellman(4, 2, 19, {2: 1, 3: 2}) is None
assert pohlig_hellman(4, 16, 19, {2: 1, 3: 2}) == 2
def test_edge_cases() -> None:
f22 = {2: 1, 11: 1}
assert pohlig_hellman(5, 1, 23, f22) == 0
assert pohlig_hellman(5, 5, 23, f22) == 1
assert pohlig_hellman(1, 1, 23, f22) == 0
assert pohlig_hellman(1, 5, 23, f22) is None
assert pohlig_hellman(22, 22, 23, f22) == 1
assert pohlig_hellman(5, 0, 23, f22) is None
assert pohlig_hellman(1, 1, 2, {}) == 0
assert pohlig_hellman(1, 0, 2, {}) is None
# a multiple of the order is enough: 2^4 * 3 * 11 * 7 is a multiple of 22
assert pohlig_hellman(5, 8, 23, {2: 4, 3: 1, 7: 1, 11: 1}) == 6
# 2^16 + 1 is prime and the group order is a power of two
assert pohlig_hellman(3, pow(3, 40000, 65537), 65537, {2: 16}) == 40000
with pytest.raises(ValueError):
pohlig_hellman(1, 1, 1, {})
def test_every_target_small_primes() -> None:
for p in [3, 5, 7, 13, 17, 31, 37, 73, 97]:
fp = factor(p - 1)
for g in range(1, p):
tab = power_table(g, p)
for h in range(p):
assert pohlig_hellman(g, h, p, fp) == tab.get(h), (g, h, p)
def test_random_primes_against_brute_force() -> None:
rng = random.Random(1978)
primes = [q for q in range(3, 100000) if is_prime(q)]
cases = 0
for trial in range(300):
p = rng.choice(primes)
fp = factor(p - 1)
r = rng.randrange(2, p)
if trial % 3 == 0:
g = r
else:
divs = [1]
for q, e in fp.items():
divs = [d * q ** i for d in divs for i in range(e + 1)]
g = pow(r, (p - 1) // rng.choice(divs), p)
tab = power_table(g, p)
for h in [1, g, rng.randrange(0, p), rng.randrange(1, p), pow(g, rng.randrange(len(tab)), p)]:
assert pohlig_hellman(g, h, p, fp) == tab.get(h), (g, h, p)
cases += 1
assert cases >= 1500
def test_composite_moduli_non_cyclic() -> None:
rng = random.Random(15)
for _ in range(400):
n = rng.randrange(2, 2500)
g = rng.randrange(1, n)
if gcd(g, n) != 1:
continue
fn = factor(phi(n))
tab = power_table(g, n)
for h in [1, g, rng.randrange(0, n), pow(g, rng.randrange(0, 999), n)]:
assert pohlig_hellman(g, h, n, fn) == tab.get(h % n), (g, h, n)
# (Z/15Z)^* is Z/2 x Z/4: 4 has order 2, and 11 has order 2 but is not a power of 4.
assert pohlig_hellman(4, 11, 15, {2: 3}) is None
assert pohlig_hellman(4, 4, 15, {2: 3}) == 1
What the next level needs
- the tests did not run in this build