Chapter 28: Diffie-Hellman Key Exchange and the Discrete Logarithm Problem

Chapter 28: Diffie-Hellman Key Exchange and the Discrete Logarithm Problem#

28.1 Historical Context: A Revolution in Cryptography#

In 1976, Whitfield Diffie and Martin Hellman published New Directions in Cryptography (IEEE Transactions on Information Theory, Vol. 22, No. 6, pp. 644–654), a paper that fundamentally transformed the field of cryptography. Before their work, all practical cryptographic systems were symmetric: the sender and receiver had to share a secret key in advance, typically through a secure physical channel. This created a serious logistical problem – how could two parties who had never met establish a shared secret over an insecure communication channel?

Diffie and Hellman proposed an elegant solution: a key exchange protocol that allows two parties to agree on a shared secret by exchanging information over a public channel, without any prior shared secrets. The security of their protocol rests on a mathematical problem called the discrete logarithm problem (DLP), which is believed to be computationally intractable for appropriately chosen groups.

Their insight was to exploit the asymmetry between exponentiation and logarithm in certain algebraic structures. In the multiplicative group \(\mathbb{Z}_p^*\) (integers modulo a prime \(p\)), computing \(g^x \bmod p\) given \(g\) and \(x\) is efficient (via repeated squaring), but recovering \(x\) from \(g^x \bmod p\) – the discrete logarithm – appears to require time that grows super-polynomially in the size of \(p\).

Note

The concept of public-key cryptography was independently conceived by James Ellis, Clifford Cocks, and Malcolm Williamson at GCHQ (the British signals intelligence agency) as early as 1969–1973, but their work remained classified until 1997.

In 1979, Leonard Adleman made foundational contributions to the analysis of the discrete logarithm problem, establishing key complexity-theoretic results about its hardness and connecting it to other problems in computational number theory.

28.2 The Discrete Logarithm Problem#

Definition 28.1 (Discrete Logarithm Problem in \(\mathbb{Z}_p^*\))

Let \(p\) be a prime and \(g\) a generator of the cyclic group \(\mathbb{Z}_p^* = \{1, 2, \ldots, p-1\}\) under multiplication modulo \(p\). Given \(h \in \mathbb{Z}_p^*\), the discrete logarithm problem is to find the unique integer \(x\) with \(0 \le x \le p-2\) such that

\[ g^x \equiv h \pmod{p}.\]

We write \(x = \log_g h\) and call \(x\) the discrete logarithm of \(h\) with respect to base \(g\).

The group \(\mathbb{Z}_p^*\) has order \(p - 1\). When \(g\) is a primitive root modulo \(p\) (i.e., a generator of \(\mathbb{Z}_p^*\)), every element \(h \in \mathbb{Z}_p^*\) has a unique discrete logarithm in \(\{0, 1, \ldots, p-2\}\).

The DLP is easy in some groups (e.g., the additive group \(\mathbb{Z}_n\) where it reduces to division) but believed to be hard in \(\mathbb{Z}_p^*\) for large primes \(p\). The best known algorithms for solving the general DLP in \(\mathbb{Z}_p^*\) are:

Algorithm	Time Complexity
Trial multiplication	\(O(p)\)
Baby-step giant-step	\(O(\sqrt{p})\) time and space
Pollard’s rho	\(O(\sqrt{p})\) time, \(O(1)\) space
Pohlig-Hellman	\(O(\sum \sqrt{p_i})\) for order \(\prod p_i^{e_i}\)
Index calculus	\(L_p[1/2, c]\) (sub-exponential)
Number field sieve	\(L_p[1/3, c]\) (sub-exponential)

28.3 The Diffie-Hellman Key Exchange Protocol#

Definition 28.2 (Diffie-Hellman Key Exchange Protocol)

Let \(p\) be a large prime and \(g\) a generator of \(\mathbb{Z}_p^*\). Both \(p\) and \(g\) are public parameters. The protocol proceeds as follows:

Alice chooses a secret integer \(a\) uniformly at random from \(\{1, 2, \ldots, p-2\}\) and computes \(A = g^a \bmod p\). She sends \(A\) to Bob.
Bob chooses a secret integer \(b\) uniformly at random from \(\{1, 2, \ldots, p-2\}\) and computes \(B = g^b \bmod p\). He sends \(B\) to Alice.
Alice computes \(s = B^a \bmod p\).
Bob computes \(s = A^b \bmod p\).

Both parties now share the common secret \(s = g^{ab} \bmod p\).

Theorem 28.1 (Correctness of Diffie-Hellman Key Exchange)

In the Diffie-Hellman protocol, both Alice and Bob compute the same shared secret:

\[ B^a \equiv (g^b)^a \equiv g^{ba} \equiv g^{ab} \equiv (g^a)^b \equiv A^b \pmod{p}.\]

This follows directly from the commutativity of integer multiplication (\(ab = ba\)) and the laws of exponentiation in \(\mathbb{Z}_p^*\).

An eavesdropper (Eve) observes \(p\), \(g\), \(A = g^a \bmod p\), and \(B = g^b \bmod p\), but does not know \(a\) or \(b\). To compute the shared secret \(g^{ab} \bmod p\), Eve would need to solve the Computational Diffie-Hellman (CDH) problem.

28.4 Computational Assumptions#

Definition 28.3 (CDH and DDH Assumptions)

Let \(G = \langle g \rangle\) be a cyclic group of prime order \(q\).

Computational Diffie-Hellman (CDH) Assumption: Given \(g\), \(g^a\), and \(g^b\) for random \(a, b \in \mathbb{Z}_q\), it is computationally infeasible to compute \(g^{ab}\).

Decisional Diffie-Hellman (DDH) Assumption: Given \(g\), \(g^a\), \(g^b\), and an element \(h \in G\), it is computationally infeasible to determine whether \(h = g^{ab}\) or \(h\) is a random element of \(G\).

These assumptions are related by the following hierarchy of reductions:

\[ \text{DLP is hard} \;\Longrightarrow\; \text{CDH is hard} \;\Longrightarrow\; \text{DDH is hard}\]

The converse implications are not known to hold in general. In other words, it is conceivable (though not proven) that one could break the Diffie-Hellman key exchange without solving the discrete logarithm problem.

28.5 Core Utility: Modular Exponentiation#

Before implementing the full protocol, we need an efficient modular exponentiation routine. The square-and-multiply algorithm computes \(b^e \bmod m\) in \(O(\log e)\) multiplications.

import numpy as np

def mod_pow(base, exp, mod):
    """
    Compute base^exp mod mod using the square-and-multiply algorithm.
    
    Uses Python's arbitrary-precision integers to avoid overflow.
    Time complexity: O(log exp) multiplications modulo mod.
    
    Parameters
    ----------
    base : int
        The base of the exponentiation.
    exp : int
        The exponent (non-negative integer).
    mod : int
        The modulus (positive integer).
    
    Returns
    -------
    int
        base^exp mod mod.
    """
    base = int(base) % int(mod)
    exp = int(exp)
    mod = int(mod)
    
    if mod == 1:
        return 0
    
    result = 1
    while exp > 0:
        if exp & 1:            # if current bit is 1: multiply
            result = (result * base) % mod
        exp >>= 1              # shift exponent right
        base = (base * base) % mod  # square the base
    return result

# ---------- Verification ----------
# Compare with Python's built-in pow(base, exp, mod)
test_cases = [(2, 10, 1000), (3, 255, 997), (7, 1000, 104729), (2, 0, 5)]
for b, e, m in test_cases:
    ours = mod_pow(b, e, m)
    builtin = pow(b, e, m)
    assert ours == builtin, f"Mismatch: mod_pow({b},{e},{m})={ours}, pow={builtin}"
    print(f"mod_pow({b}, {e}, {m}) = {ours}  [verified]")

mod_pow(2, 10, 1000) = 24  [verified]
mod_pow(3, 255, 997) = 268  [verified]
mod_pow(7, 1000, 104729) = 83336  [verified]
mod_pow(2, 0, 5) = 1  [verified]

28.6 Finding Primitive Roots#

For the Diffie-Hellman protocol we need a generator (primitive root) of \(\mathbb{Z}_p^*\). An element \(g\) is a primitive root modulo \(p\) if and only if \(g^{(p-1)/q} \not\equiv 1 \pmod{p}\) for every prime factor \(q\) of \(p - 1\).

import numpy as np
import math

def trial_factor(n):
    """
    Return the prime factorization of n as a dict {prime: exponent}.
    Uses trial division -- suitable for moderate-sized n.
    """
    factors = {}
    d = 2
    while d * d <= n:
        while n % d == 0:
            factors[d] = factors.get(d, 0) + 1
            n //= d
        d += 1
    if n > 1:
        factors[n] = factors.get(n, 0) + 1
    return factors


def is_primitive_root(g, p, factorization=None):
    """
    Test whether g is a primitive root modulo p.
    
    Parameters
    ----------
    g : int
        Candidate generator.
    p : int
        Prime modulus.
    factorization : dict or None
        Prime factorization of p-1. Computed if not given.
    
    Returns
    -------
    bool
    """
    if factorization is None:
        factorization = trial_factor(p - 1)
    order = p - 1
    for q in factorization:
        if pow(g, order // q, p) == 1:
            return False
    return True


def find_primitive_root(p):
    """
    Find the smallest primitive root modulo prime p.
    """
    factorization = trial_factor(p - 1)
    for g in range(2, p):
        if is_primitive_root(g, p, factorization):
            return g
    raise ValueError(f"No primitive root found for p={p}")

# ---------- Demo ----------
for p in [7, 23, 101, 1009]:
    g = find_primitive_root(p)
    f = trial_factor(p - 1)
    print(f"p = {int(p):>5d},  smallest primitive root g = {g},  "
          f"factorization of p-1 = {dict(f)}")

p =     7,  smallest primitive root g = 3,  factorization of p-1 = {2: 1, 3: 1}
p =    23,  smallest primitive root g = 5,  factorization of p-1 = {2: 1, 11: 1}
p =   101,  smallest primitive root g = 2,  factorization of p-1 = {2: 2, 5: 2}
p =  1009,  smallest primitive root g = 11,  factorization of p-1 = {2: 4, 3: 2, 7: 1}

28.7 The `DiffieHellman` Class#

We now implement the full Diffie-Hellman key exchange as a reusable class.

import numpy as np
import math

def mod_pow(base, exp, mod):
    """Square-and-multiply modular exponentiation."""
    base, exp, mod = int(base), int(exp), int(mod)
    if mod == 1:
        return 0
    result = 1
    base %= mod
    while exp > 0:
        if exp & 1:
            result = (result * base) % mod
        exp >>= 1
        base = (base * base) % mod
    return result


def trial_factor(n):
    """Return the prime factorization of n as {prime: exponent}."""
    factors = {}
    d = 2
    while d * d <= n:
        while n % d == 0:
            factors[d] = factors.get(d, 0) + 1
            n //= d
        d += 1
    if n > 1:
        factors[n] = factors.get(n, 0) + 1
    return factors


def is_primitive_root(g, p, factorization=None):
    """Test whether g is a primitive root modulo p."""
    if factorization is None:
        factorization = trial_factor(p - 1)
    order = p - 1
    for q in factorization:
        if pow(g, order // q, p) == 1:
            return False
    return True


def find_primitive_root(p):
    """Find the smallest primitive root modulo prime p."""
    factorization = trial_factor(p - 1)
    for g in range(2, p):
        if is_primitive_root(g, p, factorization):
            return g
    raise ValueError(f"No primitive root found for p={p}")


class DiffieHellman:
    """
    Diffie-Hellman key exchange over Z_p*.
    
    Usage
    -----
    >>> dh = DiffieHellman.generate_params(bit_length=10)
    >>> alice_priv, alice_pub = dh.generate_keypair()
    >>> bob_priv, bob_pub = dh.generate_keypair()
    >>> s1 = dh.compute_shared_secret(alice_priv, bob_pub)
    >>> s2 = dh.compute_shared_secret(bob_priv, alice_pub)
    >>> assert s1 == s2
    """
    
    def __init__(self, p, g):
        """
        Initialize with public parameters.
        
        Parameters
        ----------
        p : int
            Prime modulus.
        g : int
            Generator of Z_p*.
        """
        self.p = p
        self.g = g
    
    @classmethod
    def generate_params(cls, p=None, bit_length=16):
        """
        Create a DiffieHellman instance with public parameters.
        
        If p is not provided, searches for a prime near 2^bit_length.
        
        Parameters
        ----------
        p : int or None
            A prime modulus. If None, one is generated.
        bit_length : int
            Approximate bit length of the prime (used if p is None).
        
        Returns
        -------
        DiffieHellman
        """
        if p is None:
            rng = np.random.default_rng()
            candidate = int(rng.integers(2**(bit_length - 1), 2**bit_length))
            if candidate % 2 == 0:
                candidate += 1
            while not cls._is_prime(candidate):
                candidate += 2
            p = candidate
        g = find_primitive_root(p)
        return cls(p, g)
    
    @staticmethod
    def _is_prime(n, trials=20):
        """Miller-Rabin primality test."""
        if n < 2:
            return False
        if n < 4:
            return True
        if n % 2 == 0:
            return False
        # Write n-1 = 2^r * d
        r, d = 0, n - 1
        while d % 2 == 0:
            r += 1
            d //= 2
        rng = np.random.default_rng(42)
        for _ in range(trials):
            a = int(rng.integers(2, n - 1))
            x = pow(a, d, n)
            if x == 1 or x == n - 1:
                continue
            for _ in range(r - 1):
                x = pow(x, 2, n)
                if x == n - 1:
                    break
            else:
                return False
        return True
    
    def generate_keypair(self, rng=None):
        """
        Generate a (private_key, public_key) pair.
        
        Parameters
        ----------
        rng : np.random.Generator or None
            Random number generator. If None, a new one is created.
        
        Returns
        -------
        (int, int)
            (private_key, public_key) where public_key = g^private_key mod p.
        """
        if rng is None:
            rng = np.random.default_rng()
        private_key = int(rng.integers(2, self.p - 1))
        public_key = mod_pow(self.g, private_key, self.p)
        return private_key, public_key
    
    def compute_shared_secret(self, my_private, other_public):
        """
        Compute the shared secret from own private key and other party's public key.
        
        Parameters
        ----------
        my_private : int
            Own private key.
        other_public : int
            Other party's public key.
        
        Returns
        -------
        int
            Shared secret g^(ab) mod p.
        """
        return mod_pow(other_public, my_private, self.p)
    
    def __repr__(self):
        return (f"DiffieHellman(p={self.p} [{self.p.bit_length()} bits], "
                f"g={self.g})")


# ---------- Quick test ----------
dh = DiffieHellman.generate_params(p=104729)
print(f"Public parameters: {dh}")

rng_a = np.random.default_rng(1)
rng_b = np.random.default_rng(2)

a_priv, a_pub = dh.generate_keypair(rng_a)
b_priv, b_pub = dh.generate_keypair(rng_b)

s_alice = dh.compute_shared_secret(a_priv, b_pub)
s_bob = dh.compute_shared_secret(b_priv, a_pub)

print(f"\nAlice:  private = {a_priv},  public = {a_pub}")
print(f"Bob:    private = {b_priv},  public = {b_pub}")
print(f"\nAlice's shared secret: {s_alice}")
print(f"Bob's   shared secret: {s_bob}")
print(f"Secrets match: {s_alice == s_bob}")

Public parameters: DiffieHellman(p=104729 [17 bits], g=12)

Alice:  private = 49557,  public = 35957
Bob:    private = 87717,  public = 62914

Alice's shared secret: 14120
Bob's   shared secret: 14120
Secrets match: True

28.8 DH Key Exchange Demo with Various Parameter Sizes#

We now demonstrate the protocol with primes of increasing size to illustrate that the computation remains efficient even as parameters grow.

import numpy as np
import math
import time

def mod_pow(base, exp, mod):
    """Square-and-multiply modular exponentiation."""
    base, exp, mod = int(base), int(exp), int(mod)
    if mod == 1:
        return 0
    result = 1
    base %= mod
    while exp > 0:
        if exp & 1:
            result = (result * base) % mod
        exp >>= 1
        base = (base * base) % mod
    return result

# A selection of primes of increasing size with known primitive roots
# Format: (p, g) where g is a primitive root mod p
test_params = [
    (23, 5),           # ~5 bits
    (1009, 11),        # ~10 bits
    (104729, 12),      # ~17 bits
    (15485863, 2),     # ~24 bits
    (2147483647, 7),   # ~31 bits (Mersenne prime 2^31 - 1)
]

rng = np.random.default_rng(2024)

print(f"{'Prime p':>15s} | {'Bits':>5s} | {'g':>3s} | {'Shared secret':>15s} | {'Time (ms)':>10s}")
print("-" * 65)

for p, g in test_params:
    a = int(rng.integers(2, p - 1))
    b = int(rng.integers(2, p - 1))
    
    t0 = time.perf_counter()
    A = mod_pow(g, a, p)
    B = mod_pow(g, b, p)
    s_a = mod_pow(B, a, p)
    s_b = mod_pow(A, b, p)
    t1 = time.perf_counter()
    
    assert s_a == s_b, "Shared secrets do not match!"
    bits = p.bit_length()
    print(f"{int(p):>15d} | {int(bits):>5d} | {int(g):>3d} | {int(s_a):>15d} | {float((t1-t0)*1000):>10.4f}")

        Prime p |  Bits |   g |   Shared secret |  Time (ms)
-----------------------------------------------------------------
             23 |     5 |   5 |               2 |     0.1840
           1009 |    10 |  11 |             300 |     0.0070
         104729 |    17 |  12 |           10651 |     0.0197
       15485863 |    24 |   2 |        11186323 |     0.0125
     2147483647 |    31 |   7 |      1782278014 |     0.0544

Note

The timings above show that modular exponentiation remains extremely fast even for 31-bit primes. In practice, DH uses primes of 2048 bits or more, where Python’s built-in pow(base, exp, mod) (which uses the same square-and-multiply algorithm with GMP optimizations) is the standard choice.

28.9 Baby-Step Giant-Step (BSGS) Algorithm#

The Baby-step Giant-step algorithm, attributed to Daniel Shanks, solves the DLP in \(O(\sqrt{n})\) time and space, where \(n\) is the order of the group. The idea is to write \(x = im + j\) where \(m = \lceil\sqrt{n}\rceil\), \(0 \le j < m\), and \(0 \le i < m\). Then:

\[ g^x = g^{im + j} \implies h \cdot g^{-j} = (g^m)^i.\]

We precompute the “baby steps” \(\{(h \cdot g^{-j}, j) : 0 \le j < m\}\) in a hash table, then search for a match among the “giant steps” \((g^m)^i\) for \(i = 0, 1, \ldots, m-1\).

import numpy as np
import math

def mod_pow(base, exp, mod):
    """Square-and-multiply modular exponentiation."""
    base, exp, mod = int(base), int(exp), int(mod)
    if mod == 1:
        return 0
    result = 1
    base %= mod
    while exp > 0:
        if exp & 1:
            result = (result * base) % mod
        exp >>= 1
        base = (base * base) % mod
    return result


def baby_step_giant_step(g, h, p):
    """
    Solve g^x = h (mod p) for x using Baby-step Giant-step.
    
    Parameters
    ----------
    g : int
        Generator (base).
    h : int
        Target value.
    p : int
        Prime modulus. The group order is n = p - 1.
    
    Returns
    -------
    int or None
        The discrete logarithm x, or None if not found.
    """
    n = p - 1  # group order
    m = math.isqrt(n) + 1  # ceiling of sqrt(n)
    
    # Baby steps: compute h * g^{-j} mod p for j = 0, 1, ..., m-1
    # g^{-1} mod p = g^{p-2} mod p  (by Fermat's little theorem)
    g_inv = mod_pow(g, p - 2, p)
    
    baby_steps = {}    # value -> j
    gamma = int(h) % p
    for j in range(m):
        baby_steps[gamma] = j
        gamma = (gamma * g_inv) % p
    
    # Giant steps: compute (g^m)^i mod p for i = 0, 1, ..., m-1
    g_m = mod_pow(g, m, p)
    gamma = 1
    for i in range(m):
        if gamma in baby_steps:
            x = (i * m + baby_steps[gamma]) % n
            # Verify
            if mod_pow(g, x, p) == h % p:
                return x
        gamma = (gamma * g_m) % p
    
    return None  # Should not happen for valid inputs


# ---------- Test BSGS ----------
test_cases = [
    (2, 8, 13),    # 2^x = 8 (mod 13) => x = 3
    (3, 13, 17),   # 3^x = 13 (mod 17) => x = ?
    (5, 3, 23),    # 5^x = 3 (mod 23) => x = ?
]

for g, h, p in test_cases:
    x = baby_step_giant_step(g, h, p)
    check = mod_pow(g, x, p)
    print(f"BSGS: {g}^x = {h} (mod {p})  =>  x = {x}  "
          f"[verify: {g}^{x} = {check} (mod {p}), correct={check == h}]")

BSGS: 2^x = 8 (mod 13)  =>  x = 3  [verify: 2^3 = 8 (mod 13), correct=True]
BSGS: 3^x = 13 (mod 17)  =>  x = 4  [verify: 3^4 = 13 (mod 17), correct=True]
BSGS: 5^x = 3 (mod 23)  =>  x = 16  [verify: 5^16 = 3 (mod 23), correct=True]

28.10 BSGS Attack: Timing vs Group Order#

We now measure how the time to solve the DLP with BSGS scales with the size of the prime. The expected complexity is \(O(\sqrt{p})\), so we should see the time grow roughly with \(\sqrt{p}\).

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
import math
import time

def mod_pow(base, exp, mod):
    """Square-and-multiply modular exponentiation."""
    base, exp, mod = int(base), int(exp), int(mod)
    if mod == 1:
        return 0
    result = 1
    base %= mod
    while exp > 0:
        if exp & 1:
            result = (result * base) % mod
        exp >>= 1
        base = (base * base) % mod
    return result


def baby_step_giant_step(g, h, p):
    """Solve g^x = h (mod p) using BSGS. Returns x or None."""
    n = p - 1
    m = math.isqrt(n) + 1
    g_inv = mod_pow(g, p - 2, p)
    baby_steps = {}
    gamma = int(h) % p
    for j in range(m):
        baby_steps[gamma] = j
        gamma = (gamma * g_inv) % p
    g_m = mod_pow(g, m, p)
    gamma = 1
    for i in range(m):
        if gamma in baby_steps:
            x = (i * m + baby_steps[gamma]) % n
            if mod_pow(g, x, p) == h % p:
                return x
        gamma = (gamma * g_m) % p
    return None


# Primes of increasing size (with known small primitive roots)
primes_and_gens = [
    (101, 2), (509, 2), (1009, 11), (5003, 2), (10007, 5),
    (50021, 2), (100003, 2), (500009, 2), (1000003, 2),
    (5000017, 3), (10000019, 2),
]

rng = np.random.default_rng(42)
prime_values = []
timings = []

print(f"{'p':>10s} | {'bits':>5s} | {'sqrt(p)':>10s} | {'Time (ms)':>10s} | {'x found':>10s}")
print("-" * 55)

for p, g in primes_and_gens:
    # Pick a random target
    x_true = int(rng.integers(1, p - 1))
    h = mod_pow(g, x_true, p)
    
    t0 = time.perf_counter()
    x_found = baby_step_giant_step(g, h, p)
    t1 = time.perf_counter()
    
    elapsed_ms = (t1 - t0) * 1000
    prime_values.append(p)
    timings.append(elapsed_ms)
    x_str = f"{int(x_found):>10d}" if x_found is not None else "     None"
    print(f"{int(p):>10d} | {int(p.bit_length()):>5d} | {int(math.isqrt(p)):>10d} | {float(elapsed_ms):>10.3f} | {x_str}")

# ---------- Plot ----------
fig, ax = plt.subplots(1, 1, figsize=(8, 5))
ax.loglog(prime_values, timings, 'o-', color='#2196F3', linewidth=2,
          markersize=7, label='BSGS measured time')

# Overlay sqrt(p) reference line
p_arr = np.array(prime_values, dtype=float)
sqrt_ref = np.sqrt(p_arr)
# Scale reference line to match the data range
scale = timings[len(timings)//2] / sqrt_ref[len(sqrt_ref)//2]
ax.loglog(prime_values, sqrt_ref * scale, '--', color='#F44336', linewidth=1.5,
          label=r'$O(\sqrt{p})$ reference')

ax.set_xlabel('Prime $p$ (log scale)', fontsize=12)
ax.set_ylabel('Time (ms, log scale)', fontsize=12)
ax.set_title('BSGS Attack: Timing vs Group Order', fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('bsgs_timing.png', dpi=150, bbox_inches='tight')
plt.show()

         p |  bits |    sqrt(p) |  Time (ms) |    x found
-------------------------------------------------------
       101 |     7 |         10 |      0.008 |          9
       509 |     9 |         22 |      0.011 |        393
      1009 |    10 |         31 |      0.011 |        660
      5003 |    13 |         70 |      0.014 |       2195
     10007 |    14 |        100 |      0.017 |       4333
     50021 |    16 |        223 |      0.036 |      42947
    100003 |    17 |        316 |      0.037 |       8595
    500009 |    19 |        707 |      0.088 |      98685
   1000003 |    20 |       1000 |      0.112 |     201470
   5000017 |    23 |       2236 |      0.416 |      None
  10000019 |    24 |       3162 |      0.343 |     979077

../_images/d53313fddbc0c6c15320bf7ed569fe15dd2550ee4a28d8fc51d66057b487b497.png

Note

The log-log plot above confirms that BSGS runtime scales as \(O(\sqrt{p})\). For a 2048-bit prime, \(\sqrt{p} \approx 2^{1024}\), which is astronomically large – well beyond the reach of any computer. This is why the DLP in \(\mathbb{Z}_p^*\) is considered computationally infeasible for cryptographic-size primes.

28.11 The Pohlig-Hellman Algorithm#

When the group order \(n = p - 1\) has only small prime factors (i.e., \(n\) is smooth), the DLP can be solved much more efficiently using the Pohlig-Hellman algorithm. If \(n = \prod q_i^{e_i}\), the algorithm:

Reduces the DLP to sub-problems in subgroups of order \(q_i^{e_i}\).
Solves each sub-problem (often with BSGS in \(O(\sqrt{q_i})\)).
Combines the results using the Chinese Remainder Theorem (CRT).

The total complexity is \(O\!\left(\sum_i e_i\left(\log n + \sqrt{q_i}\right)\right)\), which is vastly faster than \(O(\sqrt{n})\) when all \(q_i\) are small.

Warning

This is precisely why Diffie-Hellman parameters must use safe primes \(p = 2q + 1\) where \(q\) is also prime. In that case, \(p - 1 = 2q\) has a large prime factor \(q\), making Pohlig-Hellman no faster than BSGS in the full group.

import numpy as np
import math

def mod_pow(base, exp, mod):
    """Square-and-multiply modular exponentiation."""
    base, exp, mod = int(base), int(exp), int(mod)
    if mod == 1:
        return 0
    result = 1
    base %= mod
    while exp > 0:
        if exp & 1:
            result = (result * base) % mod
        exp >>= 1
        base = (base * base) % mod
    return result


def baby_step_giant_step_order(g, h, p, order):
    """
    Solve g^x = h (mod p) where g has the specified order.
    Uses brute force for small orders, BSGS for larger ones.
    """
    # For small orders, brute force is faster and more reliable
    if order <= 256:
        power = 1
        for x in range(order):
            if power == h % p:
                return x
            power = (power * g) % p
        return None
    
    m = math.isqrt(order) + 1
    g_inv = mod_pow(g, p - 2, p)
    baby = {}
    gamma = int(h) % p
    for j in range(m):
        baby[gamma] = j
        gamma = (gamma * g_inv) % p
    g_m = mod_pow(g, m, p)
    gamma = 1
    for i in range(m):
        if gamma in baby:
            x = (i * m + baby[gamma]) % order
            if mod_pow(g, x, p) == h % p:
                return x
        gamma = (gamma * g_m) % p
    return None


def extended_gcd(a, b):
    """Extended Euclidean algorithm. Returns (gcd, x, y) with a*x + b*y = gcd."""
    if a == 0:
        return b, 0, 1
    g, x1, y1 = extended_gcd(b % a, a)
    return g, y1 - (b // a) * x1, x1


def crt(residues, moduli):
    """
    Chinese Remainder Theorem: find x such that
    x = residues[i] (mod moduli[i]) for all i.
    Moduli must be pairwise coprime.
    """
    M = 1
    for m in moduli:
        M *= m
    x = 0
    for r, m in zip(residues, moduli):
        Mi = M // m
        _, inv, _ = extended_gcd(Mi % m, m)
        x += r * Mi * inv
    return x % M


def pohlig_hellman(g, h, p, factorization):
    """
    Solve g^x = h (mod p) using the Pohlig-Hellman algorithm.
    
    Parameters
    ----------
    g : int
        Generator of Z_p*.
    h : int
        Target value.
    p : int
        Prime modulus.
    factorization : dict
        Prime factorization of p-1 as {prime: exponent}.
    
    Returns
    -------
    int
        The discrete logarithm x.
    """
    n = p - 1  # group order
    residues = []
    moduli = []
    
    for qi, ei in factorization.items():
        qi_ei = qi ** ei
        # Project onto subgroup of order qi^ei
        # g_i = g^{n / qi^ei}, h_i = h^{n / qi^ei}
        exp = n // qi_ei
        g_i = mod_pow(g, exp, p)
        h_i = mod_pow(h, exp, p)
        
        # Solve DLP in subgroup of order qi^ei using the "lifting" approach
        # For simplicity, we use BSGS on the projected subgroup
        x_i = baby_step_giant_step_order(g_i, h_i, p, qi_ei)
        
        if x_i is None:
            raise ValueError(f"BSGS failed for subgroup of order {qi_ei}")
        
        residues.append(x_i)
        moduli.append(qi_ei)
    
    # Combine with CRT
    x = crt(residues, moduli)
    return x


# ---------- Test Pohlig-Hellman ----------
# Use a prime where p-1 is smooth
# p = 1009, p-1 = 1008 = 2^4 * 3^2 * 7
p = 1009
g = 11  # primitive root mod 1009
factorization = {2: 4, 3: 2, 7: 1}

# Verify: 2^4 * 3^2 * 7 = 16 * 9 * 7 = 1008 = p - 1
assert 2**4 * 3**2 * 7 == p - 1

rng = np.random.default_rng(123)
for trial in range(5):
    x_true = int(rng.integers(1, p - 1))
    h = mod_pow(g, x_true, p)
    x_found = pohlig_hellman(g, h, p, factorization)
    # x_found might differ from x_true by a multiple of ord(g),
    # but g^{x_found} should equal h
    check = mod_pow(g, x_found, p)
    print(f"Trial {trial+1}: x_true={int(x_true):>4d}, x_found={int(x_found):>4d}, "
          f"g^x_found mod p = {check}, matches h={h}: {check == h}")

Trial 1: x_true=  16, x_found=  16, g^x_found mod p = 993, matches h=993: True
Trial 2: x_true= 688, x_found= 688, g^x_found mod p = 955, matches h=955: True
Trial 3: x_true= 598, x_found= 598, g^x_found mod p = 604, matches h=604: True
Trial 4: x_true=  55, x_found=  55, g^x_found mod p = 958, matches h=958: True
Trial 5: x_true= 916, x_found= 916, g^x_found mod p = 663, matches h=663: True

28.12 Smooth-Order Vulnerability: Pohlig-Hellman in Action#

We now demonstrate the dramatic speed difference between BSGS (which treats the group as opaque) and Pohlig-Hellman (which exploits smooth group order). This demonstrates why the structure of \(p - 1\) is critical for security.

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
import math
import time

def mod_pow(base, exp, mod):
    """Square-and-multiply modular exponentiation."""
    base, exp, mod = int(base), int(exp), int(mod)
    if mod == 1:
        return 0
    result = 1
    base %= mod
    while exp > 0:
        if exp & 1:
            result = (result * base) % mod
        exp >>= 1
        base = (base * base) % mod
    return result


def trial_factor(n):
    """Return the prime factorization of n as {prime: exponent}."""
    factors = {}
    d = 2
    while d * d <= n:
        while n % d == 0:
            factors[d] = factors.get(d, 0) + 1
            n //= d
        d += 1
    if n > 1:
        factors[n] = factors.get(n, 0) + 1
    return factors


def find_primitive_root(p):
    """Find the smallest primitive root modulo prime p."""
    factorization = trial_factor(p - 1)
    for g in range(2, p):
        ok = True
        for q in factorization:
            if pow(g, (p - 1) // q, p) == 1:
                ok = False
                break
        if ok:
            return g
    raise ValueError(f"No primitive root found for p={p}")


def baby_step_giant_step(g, h, p):
    """Solve g^x = h (mod p) using BSGS."""
    n = p - 1
    m = math.isqrt(n) + 1
    g_inv = mod_pow(g, p - 2, p)
    baby = {}
    gamma = int(h) % p
    for j in range(m):
        baby[gamma] = j
        gamma = (gamma * g_inv) % p
    g_m = mod_pow(g, m, p)
    gamma = 1
    for i in range(m):
        if gamma in baby:
            x = (i * m + baby[gamma]) % n
            if mod_pow(g, x, p) == h % p:
                return x
        gamma = (gamma * g_m) % p
    return None


def baby_step_giant_step_order(g, h, p, order):
    """Solve g^x = h (mod p) where g has the specified order."""
    m = math.isqrt(order) + 1
    g_inv = mod_pow(g, p - 2, p)
    baby = {}
    gamma = int(h) % p
    for j in range(m):
        baby[gamma] = j
        gamma = (gamma * g_inv) % p
    g_m = mod_pow(g, m, p)
    gamma = 1
    for i in range(m):
        if gamma in baby:
            x = (i * m + baby[gamma]) % order
            if mod_pow(g, x, p) == h % p:
                return x
        gamma = (gamma * g_m) % p
    return None


def extended_gcd(a, b):
    """Extended Euclidean algorithm."""
    if a == 0:
        return b, 0, 1
    g, x1, y1 = extended_gcd(b % a, a)
    return g, y1 - (b // a) * x1, x1


def crt(residues, moduli):
    """Chinese Remainder Theorem."""
    M = 1
    for m in moduli:
        M *= m
    x = 0
    for r, m in zip(residues, moduli):
        Mi = M // m
        _, inv, _ = extended_gcd(Mi % m, m)
        x += r * Mi * inv
    return x % M


def pohlig_hellman(g, h, p, factorization):
    """Solve g^x = h (mod p) using Pohlig-Hellman."""
    n = p - 1
    residues = []
    moduli = []
    for qi, ei in factorization.items():
        qi_ei = qi ** ei
        exp = n // qi_ei
        g_i = mod_pow(g, exp, p)
        h_i = mod_pow(h, exp, p)
        x_i = baby_step_giant_step_order(g_i, h_i, p, qi_ei)
        if x_i is None:
            raise ValueError(f"BSGS failed for subgroup of order {qi_ei}")
        residues.append(x_i)
        moduli.append(qi_ei)
    return crt(residues, moduli)


# ---------- Compare BSGS vs Pohlig-Hellman on smooth-order primes ----------
# Primes where p-1 is smooth (all prime factors are small)
smooth_primes = [
    # p,  largest_factor_of_(p-1)
    (433, 3),           # 432 = 2^4 * 3^3
    (2017, 7),          # 2016 = 2^5 * 3^2 * 7
    (30241, 7),         # 30240 = 2^5 * 3^3 * 5 * 7
    (65521, 13),        # 65520 = 2^4 * 3^2 * 5 * 7 * 13
    (700001, 7),        # 700000 = 2^5 * 5^5 * 7
    (5005001, 13),      # 5005000 = 2^3 * 5^4 * 7 * 11 * 13
]

rng = np.random.default_rng(99)

bsgs_times = []
ph_times = []
prime_list = []

print(f"{'p':>10s} | {'p-1 factors':>30s} | {'BSGS (ms)':>10s} | {'P-H (ms)':>10s} | {'Speedup':>8s}")
print("-" * 80)

for p, _ in smooth_primes:
    g = find_primitive_root(p)
    factorization = trial_factor(p - 1)
    x_true = int(rng.integers(1, p - 1))
    h = mod_pow(g, x_true, p)
    
    # Time BSGS
    t0 = time.perf_counter()
    x_bsgs = baby_step_giant_step(g, h, p)
    t_bsgs = (time.perf_counter() - t0) * 1000
    
    # Time Pohlig-Hellman
    t0 = time.perf_counter()
    x_ph = pohlig_hellman(g, h, p, factorization)
    t_ph = (time.perf_counter() - t0) * 1000
    
    bsgs_times.append(t_bsgs)
    ph_times.append(t_ph)
    prime_list.append(p)
    
    factors_str = ' * '.join(f'{q}^{e}' if e > 1 else str(q)
                             for q, e in sorted(factorization.items()))
    speedup = t_bsgs / t_ph if t_ph > 0 else float('inf')
    print(f"{int(p):>10d} | {factors_str:>30s} | {float(t_bsgs):>10.3f} | {float(t_ph):>10.3f} | {float(speedup):>7.1f}x")

# ---------- Plot ----------
fig, ax = plt.subplots(1, 1, figsize=(9, 5))
x_pos = np.arange(len(prime_list))
width = 0.35

ax.bar(x_pos - width/2, bsgs_times, width, label='BSGS', color='#2196F3', alpha=0.85)
ax.bar(x_pos + width/2, ph_times, width, label='Pohlig-Hellman', color='#F44336', alpha=0.85)

ax.set_yscale('log')
ax.set_xticks(x_pos)
ax.set_xticklabels([str(p) for p in prime_list], rotation=30, ha='right')
ax.set_xlabel('Prime $p$ (smooth order)', fontsize=12)
ax.set_ylabel('Time (ms, log scale)', fontsize=12)
ax.set_title('DLP Attack Comparison: BSGS vs Pohlig-Hellman\n(smooth-order groups)',
             fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.savefig('pohlig_hellman_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

         p |                    p-1 factors |  BSGS (ms) |   P-H (ms) |  Speedup
--------------------------------------------------------------------------------
       433 |                      2^4 * 3^3 |      0.015 |      0.022 |     0.7x
      2017 |                  2^5 * 3^2 * 7 |      0.016 |      0.040 |     0.4x
     30241 |              2^5 * 3^3 * 5 * 7 |      0.032 |      0.043 |     0.8x
     65521 |         2^4 * 3^2 * 5 * 7 * 13 |      0.043 |      0.042 |     1.0x
    700001 |                  2^5 * 5^5 * 7 |      0.099 |      0.038 |     2.6x
   5005001 |        2^3 * 5^4 * 7 * 11 * 13 |      0.306 |      0.061 |     5.0x

../_images/819208674e083a42002bab40212b1ba573341de79c89f5302c0a7637d1ec9e9f.png

Important

The chart above shows that Pohlig-Hellman is dramatically faster than generic BSGS when \(p - 1\) is smooth. The speedup grows as the largest prime factor of \(p - 1\) stays small while \(p\) itself increases. Safe primes are one common approach to ensuring the group has a large prime-order subgroup. What matters cryptographically is that the generator has prime order \(q\) where \(q\) is large enough to resist attacks like Pohlig-Hellman and Pollard’s rho.

28.13 Safe Primes vs Unsafe Primes#

A prime \(p\) is called a safe prime if \(q = (p - 1)/2\) is also prime. In that case, \(p - 1 = 2q\) and the Pohlig-Hellman algorithm offers no advantage over BSGS (the dominant sub-problem is the DLP in a subgroup of prime order \(q \approx p/2\)).

Conversely, an unsafe prime is one where \(p - 1\) has only small prime factors. We now compare the difficulty of the DLP for safe vs unsafe primes of comparable size.

Safe primes are one common approach to ensuring a large prime-order subgroup, but they are not the only option. Standards like DSA and TLS 1.3 use Schnorr groups, where \(p - 1 = kq\) for a large prime \(q\) and a cofactor \(k\). The essential requirement is that the generator’s order is a large prime, making Pohlig-Hellman ineffective.

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
import math
import time

def mod_pow(base, exp, mod):
    """Square-and-multiply modular exponentiation."""
    base, exp, mod = int(base), int(exp), int(mod)
    if mod == 1:
        return 0
    result = 1
    base %= mod
    while exp > 0:
        if exp & 1:
            result = (result * base) % mod
        exp >>= 1
        base = (base * base) % mod
    return result


def trial_factor(n):
    """Return the prime factorization of n as {prime: exponent}."""
    factors = {}
    d = 2
    while d * d <= n:
        while n % d == 0:
            factors[d] = factors.get(d, 0) + 1
            n //= d
        d += 1
    if n > 1:
        factors[n] = factors.get(n, 0) + 1
    return factors


def find_primitive_root(p):
    """Find the smallest primitive root modulo prime p."""
    factorization = trial_factor(p - 1)
    for g in range(2, p):
        ok = True
        for q in factorization:
            if pow(g, (p - 1) // q, p) == 1:
                ok = False
                break
        if ok:
            return g
    raise ValueError(f"No primitive root found for p={p}")


def baby_step_giant_step(g, h, p):
    """Solve g^x = h (mod p) using BSGS."""
    n = p - 1
    m = math.isqrt(n) + 1
    g_inv = mod_pow(g, p - 2, p)
    baby = {}
    gamma = int(h) % p
    for j in range(m):
        baby[gamma] = j
        gamma = (gamma * g_inv) % p
    g_m = mod_pow(g, m, p)
    gamma = 1
    for i in range(m):
        if gamma in baby:
            x = (i * m + baby[gamma]) % n
            if mod_pow(g, x, p) == h % p:
                return x
        gamma = (gamma * g_m) % p
    return None


def baby_step_giant_step_order(g, h, p, order):
    """Solve g^x = h (mod p) where g has the specified order."""
    m = math.isqrt(order) + 1
    g_inv = mod_pow(g, p - 2, p)
    baby = {}
    gamma = int(h) % p
    for j in range(m):
        baby[gamma] = j
        gamma = (gamma * g_inv) % p
    g_m = mod_pow(g, m, p)
    gamma = 1
    for i in range(m):
        if gamma in baby:
            x = (i * m + baby[gamma]) % order
            if mod_pow(g, x, p) == h % p:
                return x
        gamma = (gamma * g_m) % p
    return None


def extended_gcd(a, b):
    """Extended Euclidean algorithm."""
    if a == 0:
        return b, 0, 1
    g, x1, y1 = extended_gcd(b % a, a)
    return g, y1 - (b // a) * x1, x1


def crt(residues, moduli):
    """Chinese Remainder Theorem."""
    M = 1
    for m in moduli:
        M *= m
    x = 0
    for r, m in zip(residues, moduli):
        Mi = M // m
        _, inv, _ = extended_gcd(Mi % m, m)
        x += r * Mi * inv
    return x % M


def pohlig_hellman(g, h, p, factorization):
    """Solve g^x = h (mod p) using Pohlig-Hellman."""
    n = p - 1
    residues = []
    moduli = []
    for qi, ei in factorization.items():
        qi_ei = qi ** ei
        exp = n // qi_ei
        g_i = mod_pow(g, exp, p)
        h_i = mod_pow(h, exp, p)
        x_i = baby_step_giant_step_order(g_i, h_i, p, qi_ei)
        if x_i is None:
            raise ValueError(f"BSGS failed for subgroup of order {qi_ei}")
        residues.append(x_i)
        moduli.append(qi_ei)
    return crt(residues, moduli)


# ---------- Safe primes: p = 2q + 1, q prime ----------
# These are chosen so both p and (p-1)/2 are prime
safe_primes = [23, 47, 83, 167, 359, 719, 1439, 2879, 5639, 10799]

# Unsafe primes (smooth p-1): p-1 has only small prime factors
# Manually selected: each p is prime, and p-1 is B-smooth for small B
unsafe_primes = [29, 41, 97, 181, 337, 701, 1429, 2857, 5581, 10501]

rng = np.random.default_rng(7)

safe_times = []
unsafe_times = []
safe_sizes = []
unsafe_sizes = []

print("=== SAFE PRIMES (p-1 = 2q, q prime) ===")
print(f"{'p':>8s} | {'p-1 factors':>25s} | {'P-H time (ms)':>13s}")
print("-" * 55)
for p in safe_primes:
    g = find_primitive_root(p)
    fac = trial_factor(p - 1)
    x_true = int(rng.integers(1, p - 1))
    h = mod_pow(g, x_true, p)
    t0 = time.perf_counter()
    x_found = pohlig_hellman(g, h, p, fac)
    elapsed = (time.perf_counter() - t0) * 1000
    safe_times.append(elapsed)
    safe_sizes.append(p)
    f_str = ' * '.join(f'{q}^{e}' if e > 1 else str(q) for q, e in sorted(fac.items()))
    print(f"{int(p):>8d} | {f_str:>25s} | {float(elapsed):>13.4f}")

print("\n=== UNSAFE PRIMES (smooth p-1) ===")
print(f"{'p':>8s} | {'p-1 factors':>25s} | {'P-H time (ms)':>13s}")
print("-" * 55)
for p in unsafe_primes:
    g = find_primitive_root(p)
    fac = trial_factor(p - 1)
    x_true = int(rng.integers(1, p - 1))
    h = mod_pow(g, x_true, p)
    t0 = time.perf_counter()
    x_found = pohlig_hellman(g, h, p, fac)
    elapsed = (time.perf_counter() - t0) * 1000
    unsafe_times.append(elapsed)
    unsafe_sizes.append(p)
    f_str = ' * '.join(f'{q}^{e}' if e > 1 else str(q) for q, e in sorted(fac.items()))
    print(f"{int(p):>8d} | {f_str:>25s} | {float(elapsed):>13.4f}")

# ---------- Plot ----------
fig, ax = plt.subplots(1, 1, figsize=(9, 5))
ax.semilogy(safe_sizes, safe_times, 's-', color='#4CAF50', markersize=7,
            linewidth=2, label='Safe primes ($p = 2q+1$)')
ax.semilogy(unsafe_sizes, unsafe_times, 'D-', color='#F44336', markersize=7,
            linewidth=2, label='Unsafe primes (smooth $p-1$)')

ax.set_xlabel('Prime $p$', fontsize=12)
ax.set_ylabel('Pohlig-Hellman time (ms, log scale)', fontsize=12)
ax.set_title('Safe vs Unsafe Primes: Pohlig-Hellman Attack Difficulty', fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('safe_vs_unsafe_primes.png', dpi=150, bbox_inches='tight')
plt.show()

=== SAFE PRIMES (p-1 = 2q, q prime) ===
       p |               p-1 factors | P-H time (ms)
-------------------------------------------------------
|                    2 * 11 |        0.0176
|                    2 * 23 |        0.0110
|                    2 * 41 |        0.0098
|                    2 * 83 |        0.0100
|                   2 * 179 |        0.0150
|                   2 * 359 |        0.0175
|                   2 * 719 |        0.0185
|                  2 * 1439 |        0.0177
|                  2 * 2819 |        0.0175
|                  2 * 5399 |        0.0227

=== UNSAFE PRIMES (smooth p-1) ===
       p |               p-1 factors | P-H time (ms)
-------------------------------------------------------
|                   2^2 * 7 |        0.0093
|                   2^3 * 5 |        0.0087
|                   2^5 * 3 |        0.0095
|             2^2 * 3^2 * 5 |        0.0144
|               2^4 * 3 * 7 |        0.0161
|             2^2 * 5^2 * 7 |        0.0180
|          2^2 * 3 * 7 * 17 |        0.0250
|          2^3 * 3 * 7 * 17 |        0.0258
|        2^2 * 3^2 * 5 * 31 |        0.0260
|         2^2 * 3 * 5^3 * 7 |        0.0277

../_images/f5d434947852d4ee340550adee41f81056a440e281d8d5000f7834a4f63417c4.png

28.14 Visualizing the Discrete Logarithm#

To build intuition about why the DLP is hard, we visualize the discrete exponential function \(f(x) = g^x \bmod p\) for a small prime \(p\). While the function is a bijection on \(\{0, 1, \ldots, p-2\}\), its output appears pseudo-random – there is no discernible pattern that would allow one to invert it efficiently.

../_images/45bacaaa16411c6227512f959a2529935d81a6625ba764dd886de79834045d94.png

Number of distinct values: 508 out of 508 (should be 508 since g is a primitive root)

28.15 Protocol Visualization: The DH Handshake#

The following diagram illustrates the information flow in a Diffie-Hellman key exchange, showing what is public (known to Eve) and what remains secret.

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

fig, ax = plt.subplots(figsize=(12, 7))
ax.set_xlim(0, 10)
ax.set_ylim(0, 10)
ax.axis('off')

# Colors
c_alice = '#2196F3'
c_bob = '#4CAF50'
c_public = '#FF9800'
c_secret = '#F44336'
c_shared = '#9C27B0'

# Title
ax.text(5, 9.5, 'Diffie-Hellman Key Exchange Protocol',
        ha='center', fontsize=16, fontweight='bold')

# Public parameters box
rect = mpatches.FancyBboxPatch((2.5, 8.3), 5, 0.8, boxstyle='round,pad=0.1',
                                facecolor=c_public, alpha=0.2, edgecolor=c_public)
ax.add_patch(rect)
ax.text(5, 8.7, 'Public parameters: prime $p$, generator $g$',
        ha='center', fontsize=11, color='#333')

# Alice column
ax.text(2, 7.5, 'ALICE', ha='center', fontsize=14, fontweight='bold', color=c_alice)
ax.text(2, 6.8, 'Chooses secret $a$', ha='center', fontsize=10,
        color=c_secret, style='italic')
ax.text(2, 6.2, 'Computes $A = g^a$ mod $p$', ha='center', fontsize=10)

# Bob column
ax.text(8, 7.5, 'BOB', ha='center', fontsize=14, fontweight='bold', color=c_bob)
ax.text(8, 6.8, 'Chooses secret $b$', ha='center', fontsize=10,
        color=c_secret, style='italic')
ax.text(8, 6.2, 'Computes $B = g^b$ mod $p$', ha='center', fontsize=10)

# Arrows for exchange
ax.annotate('', xy=(7.2, 5.3), xytext=(2.8, 5.3),
            arrowprops=dict(arrowstyle='->', color=c_alice, lw=2))
ax.text(5, 5.6, 'Sends $A$ (public)', ha='center', fontsize=10, color=c_alice)

ax.annotate('', xy=(2.8, 4.5), xytext=(7.2, 4.5),
            arrowprops=dict(arrowstyle='->', color=c_bob, lw=2))
ax.text(5, 4.2, 'Sends $B$ (public)', ha='center', fontsize=10, color=c_bob)

# Computation
ax.text(2, 3.3, 'Computes $s = B^a$ mod $p$', ha='center', fontsize=10)
ax.text(8, 3.3, 'Computes $s = A^b$ mod $p$', ha='center', fontsize=10)

# Shared secret
rect2 = mpatches.FancyBboxPatch((2.0, 2.0), 6, 0.8, boxstyle='round,pad=0.1',
                                 facecolor=c_shared, alpha=0.15, edgecolor=c_shared)
ax.add_patch(rect2)
ax.text(5, 2.4, 'Shared secret: $s = g^{ab}$ mod $p$',
        ha='center', fontsize=12, fontweight='bold', color=c_shared)

# Eve box
rect3 = mpatches.FancyBboxPatch((3.0, 0.5), 4, 1.0, boxstyle='round,pad=0.1',
                                 facecolor='#EEE', edgecolor='#999')
ax.add_patch(rect3)
ax.text(5, 1.2, 'EVE (eavesdropper)', ha='center', fontsize=11,
        fontweight='bold', color='#666')
ax.text(5, 0.75, 'Knows: $p, g, A, B$  |  Cannot find: $a, b, s$',
        ha='center', fontsize=9, color='#666')

plt.tight_layout()
plt.savefig('dh_protocol_diagram.png', dpi=150, bbox_inches='tight')
plt.show()

../_images/db7bf38382357bf573511fb5e5420a9590938f18df292bc4526a5329d962c64d.png

28.16 Security#

Parameters in Practice

The following table summarizes the recommended minimum key sizes for Diffie-Hellman, based on NIST and other standards bodies:

Warning

The Logjam attack (2015) demonstrated that many servers used 512-bit or 1024-bit DH parameters, and that a nation-state adversary could plausibly precompute the discrete logarithm tables for a small number of commonly used 1024-bit primes, thereby passively decrypting a large fraction of Internet traffic. This motivated the migration to 2048-bit parameters.

The Internet Key Exchange standard (IKE, RFC 2409, 1998) defines specific DH groups. Two of the most commonly used are based on \(\mathbb{Z}_p^*\) (Groups 1, 2, 5, 14–18), and more recent standards (RFC 4492) add groups based on elliptic curves which achieve equivalent security with much smaller key sizes.

28.17 Exercises#

Exercise 28.1 (DLP in additive groups). Explain why the DLP is easy in the additive group \(\mathbb{Z}_n = \{0, 1, \ldots, n-1\}\) with addition modulo \(n\). Specifically, if the generator is \(g = 1\) and the target is \(h\), describe how to find \(x\) such that \(g \cdot x \equiv h \pmod{n}\) in \(O(1)\) time. Then explain why a general additive group \(\mathbb{Z}_n\) with arbitrary generator \(g\) can be solved using the Extended Euclidean Algorithm.

Click for hint

In the additive group with generator \(g\), the “discrete logarithm” is \(x\) such that \(gx \equiv h \pmod{n}\), i.e., \(x \equiv h \cdot g^{-1} \pmod{n}\). Use the extended GCD to compute \(g^{-1} \bmod n\) when \(\gcd(g, n) = 1\).

Exercise 28.2 (Implementing DH with verification). Implement a complete Diffie-Hellman key exchange for the prime \(p = 7919\) (which is a safe prime since \((7919-1)/2 = 3959\) is also prime). Verify that:

(a) Your chosen generator \(g\) is indeed a primitive root mod \(p\).
(b) Both parties compute the same shared secret.
(c) The BSGS algorithm can recover the private keys from the public keys (since \(p\) is small).

Click for hint

To verify that \(g\) is a primitive root mod 7919, check that \(g^{(p-1)/q} \not\equiv 1 \pmod{p}\) for each prime factor \(q\) of \(p-1 = 7918 = 2 \times 3959\). So you only need to check \(g^{3959} \not\equiv 1\) and \(g^{2} \not\equiv 1\) (mod 7919).

Exercise 28.3 (Pohlig-Hellman exploration). The prime \(p = 2^{16} + 1 = 65537\) is a Fermat prime. Its predecessor \(p - 1 = 65536 = 2^{16}\) is a perfect power of 2, making it maximally smooth.

(a) Find a generator of \(\mathbb{Z}_{65537}^*\).
(b) Use Pohlig-Hellman to solve the DLP \(g^x \equiv 12345 \pmod{65537}\) and time the computation.
(c) Compare with BSGS on the same instance. What is the speedup?

Click for hint

Since \(p - 1 = 2^{16}\), the Pohlig-Hellman algorithm decomposes the DLP into 16 sub-problems in subgroups of order 2 (each trivially solvable). The total work is \(O(16) = O(\log p)\), versus BSGS at \(O(\sqrt{p}) = O(256)\).

Exercise 28.4 (Man-in-the-middle attack). The basic Diffie-Hellman protocol is vulnerable to a man-in-the-middle (MITM) attack. Describe how an active adversary Mallory can intercept and modify the communication between Alice and Bob to establish separate shared secrets with each, effectively eavesdropping on all communication. Implement a simulation of this attack.

Click for hint

Mallory intercepts \(A\) from Alice, generates her own private key \(m\), sends \(M = g^m\) to Bob (pretending to be Alice). Similarly, she intercepts \(B\) from Bob and sends \(M' = g^{m'}\) to Alice. Now Mallory shares one secret with Alice (\(g^{am}\)) and another with Bob (\(g^{m'b}\)), and can decrypt/re-encrypt all traffic.

Exercise 28.5 (Generator order and subgroup confinement). When \(p - 1 = 2q\) for a safe prime \(p\), the group \(\mathbb{Z}_p^*\) has exactly three subgroups: \(\{1\}\), the subgroup of order 2 (\(\{1, p-1\}\)), and the subgroup of order \(q\). An attacker who tricks a party into using a generator of order 2 can trivially determine the shared secret. Demonstrate this “small subgroup” attack:

(a) For \(p = 7919\), show that \(p - 1 = 7918\) is a generator of order 2.
(b) If Alice uses private key \(a\) and computes \(A = (p-1)^a \bmod p\), show that \(A\) reveals the parity of \(a\).
(c) Explain why proper implementations must verify that received public keys lie in the correct subgroup.

Click for hint

Note that \((p-1) \equiv -1 \pmod{p}\), so \((p-1)^a \equiv (-1)^a \pmod{p}\). This is \(1\) if \(a\) is even and \(p-1\) if \(a\) is odd. So \(A\) leaks one bit of information about \(a\), and the “shared secret” is always \(1\) or \(p-1\).

28.18 Summary#

In this chapter we explored the Diffie-Hellman key exchange protocol and the discrete logarithm problem upon which its security rests.

Key concepts:

The DH key exchange allows two parties to establish a shared secret over a public channel. Its correctness follows from the commutativity of exponentiation: \((g^a)^b = g^{ab} = (g^b)^a\).
The security of DH relies on the Computational Diffie-Hellman (CDH) assumption: given \(g^a\) and \(g^b\), computing \(g^{ab}\) is infeasible without knowing \(a\) or \(b\).
The Baby-step Giant-step algorithm solves the DLP in \(O(\sqrt{p})\) time and space, establishing a baseline for attack complexity.
The Pohlig-Hellman algorithm exploits the factorization of the group order. If \(p - 1\) is smooth (has only small prime factors), the DLP becomes easy. This is why DH must use safe primes \(p = 2q + 1\).
In practice, DH requires primes of at least 2048 bits. The Logjam attack demonstrated the danger of using weak parameters.
The basic DH protocol does not provide authentication and is vulnerable to man-in-the-middle attacks. In practice, it is used within authenticated protocols (TLS, IKE).

References#

W. Diffie and M. E. Hellman, “New Directions in Cryptography,” IEEE Transactions on Information Theory, vol. 22, no. 6, pp. 644–654, 1976. doi:10.1109/TIT.1976.1055638
L. M. Adleman, “A Subexponential Algorithm for the Discrete Logarithm Problem with Applications to Cryptography,” in Proc. 20th IEEE Symposium on Foundations of Computer Science, pp. 55–60, 1979.
S. Pohlig and M. Hellman, “An Improved Algorithm for Computing Logarithms over GF(p) and Its Cryptographic Significance,” IEEE Transactions on Information Theory, vol. 24, no. 1, pp. 106–110, 1978.
D. Adrian et al., “Imperfect Forward Secrecy: How Diffie-Hellman Fails in Practice” (the Logjam attack), in Proc. 22nd ACM SIGSAC Conference on Computer and Communications Security, 2015.
D. Harkins and D. Carrel, “The Internet Key Exchange (IKE),” RFC 2409, November 1998.
A. J. Menezes, P. C. van Oorschot, and S. A. Vanstone, Handbook of Applied Cryptography, Chapter 3 (“Number-Theoretic Reference Problems”), CRC Press, 1997.
D. Shanks, “Class Number, a Theory of Factorization, and Genera,” in Proc. Symposia in Pure Mathematics, vol. 20, pp. 415–440, AMS, 1971.

Chapter 28: Diffie-Hellman Key Exchange and the Discrete Logarithm Problem

Contents

Chapter 28: Diffie-Hellman Key Exchange and the Discrete Logarithm Problem#

28.1 Historical Context: A Revolution in Cryptography#

28.2 The Discrete Logarithm Problem#

28.3 The Diffie-Hellman Key Exchange Protocol#

28.4 Computational Assumptions#

28.5 Core Utility: Modular Exponentiation#

28.6 Finding Primitive Roots#

28.7 The DiffieHellman Class#

28.8 DH Key Exchange Demo with Various Parameter Sizes#

28.9 Baby-Step Giant-Step (BSGS) Algorithm#

28.10 BSGS Attack: Timing vs Group Order#

28.11 The Pohlig-Hellman Algorithm#

28.12 Smooth-Order Vulnerability: Pohlig-Hellman in Action#

28.13 Safe Primes vs Unsafe Primes#

28.14 Visualizing the Discrete Logarithm#

28.15 Protocol Visualization: The DH Handshake#

28.16 Security#

28.17 Exercises#

28.18 Summary#

References#

28.7 The `DiffieHellman` Class#