Chapter 46 — Practical Exercises

This sheet accompanies Chapter 46 — The Post-Quantum Landscape and the Knapsack Problem. Six exercises from one-liner to small research project.

How to use

Each exercise has three parts:

1. Problem statement — the goal, plus a difficulty marker (★ very easy to ★★★★★ research project) and a link to the chapter section that covers the theory.

2. Your turn — a code cell with TODO gaps. Fill them in, run the cell, and check the assertion at the bottom.

3. Solution — hidden by default. Click “Click to show” to reveal a fully-commented reference implementation that explains every step and links back to the chapter.

The gaps deliberately point at concepts the chapter teaches; doing the exercises is a more durable way to learn the theory than re-reading the text.

Working environment

All exercises run in pure Python with the standard library plus NumPy. Run pip install numpy matplotlib once if you have not already.

Exercise 46.E1 — ★ Validate a superincreasing sequence

Goal. Write is_superincreasing(w) returning True iff \(w_k > \sum_{i<k} w_i\) for every \(k \ge 2\).

Theory. §46.5 — Superincreasing knapsacks are easy.

Why this matters. Merkle–Hellman keygen requires a superincreasing secret sequence; if you accidentally generate one that is not, the greedy decryption silently returns the wrong plaintext.

def is_superincreasing(w):
    '''Return True iff w_k > sum(w_i for i < k) for every k >= 2.'''
    # TODO: walk through w accumulating a running sum and check each w_k.
    # If any element fails the condition, return False; otherwise True.
    raise NotImplementedError('your turn')


# Tests.  Uncomment and run.
# assert is_superincreasing([2, 3, 8, 17, 31, 65, 130, 258])
# assert not is_superincreasing([2, 3, 4])           # 4 == 2+3 fails strict >
# assert not is_superincreasing([5, 3])
# assert is_superincreasing([1])                     # singleton trivially yes
# print('E1 OK')

Show code cell content

Hide code cell content

# Solution.
#
# We walk through the sequence with a running prefix sum.  At step k we test
# whether the current element strictly exceeds the sum of everything before it.
#
# Reference: §46.5 of the chapter defines the property.  The greedy decoder
# in the same section relies critically on this STRICT inequality -- if we had
# w_k == sum_{i<k} w_i, two distinct subsets would produce the same target sum
# and decryption would be ambiguous.
def is_superincreasing(w):
    s = 0
    for wk in w:
        if wk <= s:
            return False
        s += wk
    return True

assert is_superincreasing([2, 3, 8, 17, 31, 65, 130, 258])
assert not is_superincreasing([2, 3, 4])
assert not is_superincreasing([5, 3])
assert is_superincreasing([1])
print('E1 OK')

E1 OK

Exercise 46.E2 — ★★ Decrypt a Merkle–Hellman ciphertext

Goal. You are given a Merkle–Hellman secret key \((w, m, r)\) and a ciphertext \(S\). Recover the plaintext bit-vector \(\mathbf{x}\).

Theory. §46.6 — The Merkle–Hellman trapdoor.

Why this matters. This is exactly what the recipient does on every incoming ciphertext. Two non-trivial steps: (1) modular un-disguise via \(r^{-1} \bmod m\), and (2) the greedy decode.

# Given key + ciphertext.
W = [2, 7, 18, 28, 63, 127, 249, 500]   # secret superincreasing sequence
M = 1045                                # secret modulus
R = 789                                 # secret multiplier
S = 1514                                # ciphertext (an ordinary integer)

def egcd(a, b):
    if b == 0: return a, 1, 0
    g, x1, y1 = egcd(b, a % b)
    return g, y1, x1 - (a // b) * y1

def modinv(a, m):
    g, x, _ = egcd(a, m)
    if g != 1: raise ValueError('not invertible')
    return x % m


def merkle_hellman_decrypt(S, w, m, r):
    '''Recover plaintext bits from ciphertext S given secret (w, m, r).'''
    # TODO step 1: compute the modular inverse r_inv of r mod m.
    # TODO step 2: compute S' = (r_inv * S) mod m.  Note: S' should now
    #              equal sum_i x_i * w_i AS AN INTEGER -- the integer sum is
    #              less than m by construction, so reducing mod m is benign.
    # TODO step 3: greedy decode S' against the SUPERINCREASING w.  Return
    #              a list of bits.
    raise NotImplementedError('your turn')


# Test: round-trip the published values.
# bits = merkle_hellman_decrypt(S, W, M, R)
# print('Recovered bits:', bits)
# assert bits == [1, 1, 0, 1, 0, 0, 1, 1]   # the answer
# print('E2 OK')

Show code cell content

Hide code cell content

# Solution.
#
# Step 1.  By construction r and m are coprime, so r has a modular inverse.
#          We use the extended Euclidean algorithm (egcd above).  This is
#          §46.6 of the chapter.
#
# Step 2.  S' = r_inv * S mod m exposes the disguise: a_i was r * w_i mod m,
#          and S = sum_i x_i * a_i (as a plain integer).  So
#                 r_inv * S = sum_i x_i * w_i mod m.
#          Crucially, sum w_i < m by construction, so the right-hand side
#          taken mod m equals the actual integer sum -- no information lost.
#
# Step 3.  Greedy decode the easy superincreasing instance: §46.5 again.
#          Walk from the largest w_k downwards and include it in the subset
#          whenever it does not exceed the remaining target.

def merkle_hellman_decrypt(S, w, m, r):
    r_inv = modinv(r, m)
    Sp = (r_inv * S) % m                              # un-disguised target
    n = len(w)
    bits = [0] * n
    rem = Sp
    for k in range(n - 1, -1, -1):
        if w[k] <= rem:
            bits[k] = 1
            rem -= w[k]
    assert rem == 0, 'decode failed -- bad key or non-Merkle-Hellman ciphertext'
    return bits

bits = merkle_hellman_decrypt(S, W, M, R)
print('Recovered bits:', bits)
assert bits == [1, 1, 0, 1, 0, 0, 1, 1]
print('E2 OK')

Recovered bits: [1, 1, 0, 1, 0, 0, 1, 1]
E2 OK

Exercise 46.E3 — ★★★ Density vs the Lagarias–Odlyzko threshold

Goal. Implement density(public_key) and predict whether Lagarias–Odlyzko (1985) lattice reduction can break the instance.

Theory. §46.7 — Why Merkle–Hellman broke and exercise §46.8 Ex 46.2.

Why this matters. Density is the single most important parameter for predicting whether LLL will recover the plaintext. Anything below \(d \approx 0.94\) is generically broken in polynomial time; above 1, LLL expects multiple-bit collisions and the attack does not directly work.

import math

def density(a):
    '''Knapsack density d = n / log2(max(a_i)).'''
    # TODO: return n divided by log2 of the maximum entry of a.
    raise NotImplementedError('your turn')


def lo_attackable(a, threshold=0.9408):
    '''Predict whether Lagarias-Odlyzko (and Coster et al. 1992) succeeds.'''
    # TODO: return True iff density(a) < threshold.
    raise NotImplementedError('your turn')


# Toy public keys at three densities.
PK_LOW  = [0xC2DE for _ in range(8)]                                # ~16-bit a_i, n=8 -> d~0.5
PK_MED  = [(1 << 12) - 1 - i * 137 for i in range(8)]               # ~12-bit a_i -> d~0.67
PK_HIGH = [3 + i for i in range(8)]                                 # tiny a_i -> d ~ 2

# print('low  d =', density(PK_LOW),  ' attackable =', lo_attackable(PK_LOW))
# print('med  d =', density(PK_MED),  ' attackable =', lo_attackable(PK_MED))
# print('high d =', density(PK_HIGH), ' attackable =', lo_attackable(PK_HIGH))

Show code cell content

Hide code cell content

# Solution.
#
# Density is just the bit budget per knapsack entry: n / log_2(max a_i).
# Lagarias-Odlyzko 1985 proved an LLL-based attack succeeds for d < 0.6463;
# Coster, Joux, LaMacchia, Odlyzko, Schnorr, Stern (1992) raised this to
# d < 0.9408.  See the chapter §46.7.
#
# Note we use ceil(log2) -- a power-of-two max would otherwise be off by one.

def density(a):
    return len(a) / math.ceil(math.log2(max(a)))

def lo_attackable(a, threshold=0.9408):
    return density(a) < threshold

print('low  d =', round(density(PK_LOW), 3),  ' attackable =', lo_attackable(PK_LOW))
print('med  d =', round(density(PK_MED), 3),  ' attackable =', lo_attackable(PK_MED))
print('high d =', round(density(PK_HIGH), 3), ' attackable =', lo_attackable(PK_HIGH))
print('E3 OK')

low  d = 0.5  attackable = True
med  d = 0.667  attackable = True
high d = 2.0  attackable = False
E3 OK

Exercise 46.E4 — ★★★ Brute force vs greedy: timing comparison

Goal. Implement an exhaustive \(O(2^n)\) subset-sum solver and time it against the greedy decoder on superincreasing instances of growing \(n\).

Theory. §46.4 — The subset-sum / 0-1 knapsack problem.

Why this matters. Convince yourself, with concrete timings, that the NP-hard generic problem and the \(O(n)\) trapdoored problem really are in totally different complexity classes.

import time, itertools

def brute_force_subset_sum(a, S):
    '''Try every subset; return the bit-vector matching S, or None.'''
    # TODO: enumerate all 2^n bit-vectors and pick the first whose dot product
    # with `a` equals S.  itertools.product([0,1], repeat=n) is one way.
    raise NotImplementedError('your turn')


def greedy_subset_sum(w, S):
    bits = [0] * len(w)
    rem = S
    for k in range(len(w) - 1, -1, -1):
        if w[k] <= rem:
            bits[k] = 1
            rem -= w[k]
    return bits if rem == 0 else None


# Generate a superincreasing sequence and time both solvers for growing n.
def superincr(n):
    w = [2]
    while len(w) < n:
        w.append(sum(w) + 1)
    return w

# for n in (8, 12, 16, 20):
#     w = superincr(n)
#     S = sum(w[:n:2])           # take the even-indexed weights
#     t0 = time.time(); b1 = brute_force_subset_sum(w, S); tb = time.time() - t0
#     t0 = time.time(); b2 = greedy_subset_sum(w, S);      tg = time.time() - t0
#     print(f'n={n:2d}  brute = {tb*1000:7.1f} ms   greedy = {tg*1000:6.3f} ms   match={b1==b2}')

Show code cell content

Hide code cell content

# Solution.
#
# Brute force enumerates all 2^n subsets.  itertools.product is the cleanest
# implementation; for serious benchmarking you would prefer a meet-in-the-
# middle algorithm (Schroeppel-Shamir 1981) which runs in O(2^{n/2}) time and
# space.
#
# Empirically you should see brute-force time roughly DOUBLE per +1 in n,
# while greedy stays essentially flat -- the textbook exponential vs linear
# divide.  This is exactly the gap that motivates the trapdoor-knapsack idea
# (chapter §46.6) and exactly what makes that idea fragile (§46.7).

def brute_force_subset_sum(a, S):
    n = len(a)
    for bits in itertools.product([0, 1], repeat=n):
        if sum(b * ai for b, ai in zip(bits, a)) == S:
            return list(bits)
    return None

for n in (8, 12, 16, 20):
    w = superincr(n)
    S = sum(w[:n:2])
    t0 = time.time(); b1 = brute_force_subset_sum(w, S); tb = time.time() - t0
    t0 = time.time(); b2 = greedy_subset_sum(w, S);      tg = time.time() - t0
    print(f'n={n:2d}  brute = {tb*1000:7.1f} ms   greedy = {tg*1000:6.3f} ms   match={b1==b2}')

print('E4 OK')

n= 8  brute =     0.1 ms   greedy =  0.003 ms   match=True
n=12  brute =     1.7 ms   greedy =  0.002 ms   match=True
n=16  brute =    29.2 ms   greedy =  0.005 ms   match=True

n=20  brute =   504.5 ms   greedy =  0.011 ms   match=True
E4 OK

Exercise 46.E5 — ★★★★ Build the Lagarias–Odlyzko lattice + a tiny LLL

Goal. Construct the Lagarias–Odlyzko (LO) lattice for a given public-key / ciphertext pair, run a NumPy implementation of LLL on it, and decode the bit-vector from the resulting short vector.

Theory. §46.7 — Why Merkle–Hellman broke and the upstream cryptanalysis chapter Ch 40 (§40.3–40.7) for the LLL algorithm itself.

Why this matters. This is the W2 lab in miniature. The LO lattice is the bridge from “knapsack as a hard problem” to “LLL as a near-universal cryptanalytic hammer”.

import numpy as np

def build_LO_lattice(a, S):
    '''Lagarias-Odlyzko lattice: returns an (n+1) x (n+1) integer basis B.

    Rows are the basis vectors; the last column carries the sum constraint.
    A short vector with last coord 0 has front entries +/- 1 encoding bits.
    '''
    n = len(a)
    # TODO: build B with the following pattern (rows are vectors):
    #   for i in 0..n-1:  row i = 2 * e_i  +  2 * a_i * e_n
    #   for i = n      :  row n = 1_n      +  2 * S    * e_n
    # Return as a 2D NumPy array with dtype=object (so coefficients can be big ints).
    raise NotImplementedError('your turn')


def lll(B, delta=0.75):
    '''Plain LLL on integer-row basis B (NumPy 2D array).  Returns reduced B.'''
    B = np.array(B, dtype=object)
    n = B.shape[0]
    def gso():
        Bf = B.astype(float)
        Q  = np.zeros_like(Bf)
        mu = np.zeros((n, n))
        for i in range(n):
            Q[i] = Bf[i].copy()
            for j in range(i):
                denom = float(Q[j] @ Q[j])
                mu[i, j] = float(Bf[i] @ Q[j]) / denom if denom else 0.0
                Q[i] -= mu[i, j] * Q[j]
        return Q, mu
    def size_reduce(k, mu):
        for j in range(k - 1, -1, -1):
            q = round(mu[k, j])
            if q:
                B[k] = B[k] - q * B[j]
                mu[k, :j+1] -= q * mu[j, :j+1]
                mu[k, j]   -= q
    Q, mu = gso()
    k = 1
    iters = 0
    while k < n and iters < 1000 * n:
        size_reduce(k, mu)
        if (Q[k] @ Q[k]) >= (delta - mu[k, k-1]**2) * (Q[k-1] @ Q[k-1]):
            k += 1
        else:
            B[[k, k-1]] = B[[k-1, k]]
            Q, mu = gso()
            k = max(k - 1, 1)
        iters += 1
    return B


def bits_from_short_vec(v, n):
    out = []
    for x in v[:n]:
        if   int(x) ==  1: out.append(0)
        elif int(x) == -1: out.append(1)
        else: return None
    return out


def attack_merkle_hellman(a, S):
    B = build_LO_lattice(a, S)
    Bred = lll(B)
    for row in Bred:
        if int(row[-1]) != 0: continue
        b = bits_from_short_vec(row, len(a))
        if b is not None: return b
        bneg = bits_from_short_vec(-row, len(a))
        if bneg is not None: return bneg
    return None


# Build a low-density instance and break it.
# from random import Random
# rnd = Random(42)
# n = 10
# w = [rnd.randrange(2, 10)]
# for _ in range(n - 1):
#     w.append(sum(w) + rnd.randrange(1, 10))
# m = (1 << 32) + rnd.randrange(1, 1024); r = 17  # tiny example
# while True:
#     try: r_inv = pow(r, -1, m); break
#     except ValueError: r += 2
# pk = [(r * wi) % m for wi in w]
# bits = [rnd.randrange(2) for _ in range(n)]
# S = sum(b * ai for b, ai in zip(bits, pk))
# rec = attack_merkle_hellman(pk, S)
# print('density =', n / max(pk).bit_length())
# print('truth     =', bits)
# print('recovered =', rec)
# assert rec == bits
# print('E5 OK')

Show code cell content

Hide code cell content

# Solution.
#
# The lattice we build is the standard Lagarias-Odlyzko construction (1985):
#
#     b_i (i<n) = (0,...,2,...,0, 2*a_i)        <- 2 in position i, 2*a_i last
#     b_n       = (1, 1, ..., 1, 2*S)
#
# A vector v in the lattice is a Z-linear combination sum_i c_i b_i.  Setting
# the last coordinate to zero forces 2*sum_i c_i a_i + 2*c_n*S = 0 with c_n =
# +/- 1; solving gives c_i = -+ x_i where x is the plaintext bit-vector, and
# the front entries become 2c_i + c_n = +/- 1, with sign matching whether the
# bit is 0 (+1) or 1 (-1).
#
# The chapter §46.7 sketches this, and the upstream cryptanalysis course
# Chapter 40 has the full LLL treatment.

def build_LO_lattice(a, S):
    n = len(a)
    B = np.zeros((n + 1, n + 1), dtype=object)
    for i in range(n):
        B[i, i] = 2
        B[i, n] = 2 * a[i]
    for j in range(n):
        B[n, j] = 1
    B[n, n] = 2 * S
    return B

# Same lll / bits_from_short_vec / attack helpers as the "your turn" cell.

from random import Random
rnd = Random(42)
n = 10
w = [rnd.randrange(2, 10)]
for _ in range(n - 1):
    w.append(sum(w) + rnd.randrange(1, 10))
m = (1 << 32) + rnd.randrange(1, 1024); r = 17
while True:
    try: r_inv = pow(r, -1, m); break
    except ValueError: r += 2
pk = [(r * wi) % m for wi in w]
bits = [rnd.randrange(2) for _ in range(n)]
S = sum(b * ai for b, ai in zip(bits, pk))
rec = attack_merkle_hellman(pk, S)
print('density   =', round(n / max(pk).bit_length(), 3))
print('truth     =', bits)
print('recovered =', rec)
assert rec == bits
print('E5 OK')

density   = 0.667
truth     = [0, 0, 0, 0, 0, 0, 1, 0, 1, 1]
recovered = [0, 0, 0, 0, 0, 0, 1, 0, 1, 1]
E5 OK

Exercise 46.E6 — ★★★★★ Research: density vs LLL success rate

Goal. Empirically map the success rate of the Lagarias–Odlyzko attack as a function of knapsack density. Generate \(T = 50\) random Merkle–Hellman instances at each density value \(d \in \{0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0\}\), run your attack_merkle_hellman from E5, and plot the recovery rate.

Theory. Lagarias–Odlyzko 1985, refined by Coster–Joux–LaMacchia–Odlyzko– Schnorr–Stern 1992 (chapter §46.7 + the upstream Ch 40).

Why this matters. This is exactly the kind of empirical question every PQC paper has to answer for its own scheme: “at what parameter regime does the best known attack succeed, and where is the margin?”. The Coster et al. (1992) bound says success below \(d \approx 0.94\), but real-world LLL is implementation-dependent — measure your own.

Open-ended

This is a small research project. Variations to try: vary \(n\) at fixed \(d\); replace LLL with BKZ-2; replace the LO lattice with the Coster et al. 1992 variant. Bring plots to office hours.

import numpy as np
from random import Random

def random_mh_instance(n, density_target, rnd):
    '''Return (pk, ct, bits) at a target density.

    Density d ~ n / log2(max a_i).  We pick the modulus m ~ 2^(n/d) and
    the multiplier r at random coprime to m, then publish a_i = r * w_i mod m.
    '''
    # TODO step 1: generate a superincreasing secret sequence w of length n.
    # TODO step 2: compute the modulus m from the density target.
    # TODO step 3: pick a multiplier r coprime to m.
    # TODO step 4: build pk and a random plaintext, return (pk, ciphertext, bits).
    raise NotImplementedError('your turn')


# Scan densities, run your attack from E5, plot recovery rate.
# (Reuse `attack_merkle_hellman` from E5 in the same notebook.)

Show code cell content

Hide code cell content

# Solution sketch.
#
# The dominant subtlety is computing the modulus to land at a target density.
# For an n-element key with d = n / log2(max(a_i)) and max(a_i) ~ m, we want
# m ~ 2^(n/d).  But m must also exceed sum(w); we take m = max(2^(n/d), sum(w))
# plus a small random offset.
#
# Plotting matplotlib in this notebook requires `pip install matplotlib`.

import math
import matplotlib.pyplot as plt   # noqa: only used in the solution

def random_mh_instance(n, density_target, rnd):
    # superincreasing secret w
    w = [rnd.randrange(2, 10)]
    while len(w) < n:
        w.append(sum(w) + rnd.randrange(1, 10))
    target_max_log = max(int(round(n / density_target)), int(math.log2(sum(w))) + 1)
    m = max(sum(w) + rnd.randrange(1, 100), 1 << target_max_log)
    while True:
        r = rnd.randrange(2, m)
        if math.gcd(r, m) == 1: break
    pk = [(r * wi) % m for wi in w]
    bits = [rnd.randrange(2) for _ in range(n)]
    S = sum(b * ai for b, ai in zip(bits, pk))
    return pk, S, bits

def success_rate(n, d, T, rnd):
    wins = 0
    for _ in range(T):
        pk, S, bits = random_mh_instance(n, d, rnd)
        rec = attack_merkle_hellman(pk, S)
        if rec == bits:
            wins += 1
    return wins / T

# Reduce T at higher n for runtime.
n = 10
T = 30
ds  = [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
rnd = Random(2026)
rates = [success_rate(n, d, T, rnd) for d in ds]
for d, p in zip(ds, rates):
    print(f'd = {d:.2f}   success {int(p*T)}/{T} = {p:.0%}')

# Plot.
fig, ax = plt.subplots(figsize=(6, 3.2))
ax.plot(ds, rates, marker='o', color='#4f46e5')
ax.axvline(0.9408, color='#dc2626', linestyle='--', label='Coster et al. 1992 threshold')
ax.set_xlabel('density d')
ax.set_ylabel(f'LLL recovery rate (n={n}, T={T})')
ax.set_ylim(-0.05, 1.05); ax.grid(True, alpha=0.3); ax.legend()
ax.set_title('Lagarias-Odlyzko attack success vs. knapsack density')
plt.tight_layout()
plt.show()

print('E6 OK')

d = 0.20   success 30/30 = 100%
d = 0.30   success 30/30 = 100%
d = 0.40   success 30/30 = 100%
d = 0.50   success 30/30 = 100%
d = 0.60   success 30/30 = 100%
d = 0.70   success 30/30 = 100%
d = 0.80   success 30/30 = 100%
d = 0.90   success 30/30 = 100%
d = 1.00   success 30/30 = 100%

E6 OK