Pohlig Hellman for Discrete Logarithms

Pohlig-Hellman.py

#Full guide: https://leetarxiv.substack.com/p/pohlig-hellman-discrete-logarithms
def pohlig_hellman(alpha, beta, p, factors):
    n = p - 1  # order of the exponent group
    congruences = []
    
    # For each prime factor q^e
    for q, e in factors:
        # Compute gamma = alpha^(n/q)
        gamma = pow(alpha, n // q, p)
        
        x_i = 0
        # For each exponent from 0 to e-1
        for j in range(e):
            # Compute beta_j = (beta * alpha^(-x_i))^(n/q^(j+1)) mod p
            numerator = (beta * pow(alpha, -x_i, p)) % p
            exponent = n // (q ** (j + 1))
            beta_j = pow(numerator, exponent, p)
            
            # Find k such that gamma^k ≡ beta_j mod p
            for k in range(q):
                if pow(gamma, k, p) == beta_j:
                    x_i = x_i + k * (q ** j)
                    break
        
        congruences.append((x_i, q ** e))
    
    # Solve system of congruences using Chinese Remainder Theorem
    return chinese_remainder(congruences)

def chinese_remainder(congruences):
    """Solve system of congruences x ≡ a_i mod m_i"""
    x = 0
    M = 1
    for a, m in congruences:
        M *= m
    
    for a_i, m_i in congruences:
        M_i = M // m_i
        # Find inverse of M_i mod m_i
        inv = pow(M_i, -1, m_i)
        x = (x + a_i * M_i * inv) % M
    
    return x

# Given parameters
p = 251
alpha = 71
beta = 210

# Factorize 250 = 2 * 5^3
factors = [(2, 1), (5, 3)]

# Compute the discrete logarithm
x = pohlig_hellman(alpha, beta, p, factors)

print(f"x = log_{alpha}({beta}) mod {p} = {x}")
print(f"Verification: {alpha}^{x} mod {p} = {pow(alpha, x, p)}")