Chapter 48 — Practical Exercises

Companion sheet for Chapter 48 — NTRU, Multivariate, and Isogeny Cryptography. Six exercises spanning all three minority families.

Exercise 48.E1 — ★ Sample short ternary polynomials

Goal. Write sample_ternary(N, n_plus, n_minus) returning a length-\(N\) list with exactly n_plus entries equal to \(+1\), n_minus entries equal to \(-1\), and the rest zero, in random order.

Theory. §48.2 — NTRU keygen.

import random as _r
rng = _r.Random(20260504)

def sample_ternary(N, n_plus, n_minus):
    '''Return a length-N list with exactly n_plus +1, n_minus -1, rest 0.'''
    # TODO: build a list with the right multiplicities, then shuffle.
    raise NotImplementedError('your turn')


# Tests.
# v = sample_ternary(11, 4, 3)
# assert sum(1 for x in v if x == 1) == 4
# assert sum(1 for x in v if x == -1) == 3
# assert sum(1 for x in v if x == 0) == 4
# print('E1 OK')

Show code cell content

Hide code cell content

# Solution.
#
# Build the multiset, shuffle in place.  This is the standard ternary sampler
# used by NTRU, Falcon, and (with different multiplicities) by Kyber's CBD.

def sample_ternary(N, n_plus, n_minus):
    v = [1] * n_plus + [-1] * n_minus + [0] * (N - n_plus - n_minus)
    rng.shuffle(v)
    return v

v = sample_ternary(11, 4, 3)
assert sum(1 for x in v if x == 1) == 4
assert sum(1 for x in v if x == -1) == 3
assert sum(1 for x in v if x == 0) == 4
print('E1 OK')

E1 OK

Exercise 48.E2 — ★★ Centred reduction in \(\mathbb{Z}_q\)

Goal. Write centred(a, q) that reduces every entry of a to the unique representative in \((-q/2, q/2]\).

Theory. §48.2 — NTRU decrypt step “centred”.

Why this matters. NTRU decryption hinges on lifting a polynomial that looks random mod \(q\) back to its small-coefficient integer representation. The centred lift is the only choice that preserves the “small coefficient” property.

def centred(a, q):
    '''Reduce each entry of a to (-q/2, q/2].'''
    # TODO: take entries mod q, then subtract q from anything strictly above q/2.
    raise NotImplementedError('your turn')


# assert centred([0, 1, 15, 16, 17, 30, 31], 31) == [0, 1, 15, -15, -14, -1, 0]
# print('E2 OK')

Show code cell content

Hide code cell content

# Solution.
#
# The centred lift puts each residue class into its "smallest absolute value"
# representative.  For odd q, the half-open interval (-q/2, q/2] gives a
# unique representative; for even q, both -q/2 and q/2 represent the same
# class -- the convention is to take q/2 (positive).

def centred(a, q):
    out = []
    for x in a:
        x = x % q
        if x > q // 2:
            x -= q
        out.append(x)
    return out

assert centred([0, 1, 15, 16, 17, 30, 31], 31) == [0, 1, 15, -15, -14, -1, 0]
print('E2 OK')

E2 OK

Exercise 48.E3 — ★★ NTRU decryption-failure rate

Goal. For toy NTRU with \(N = 11\), \(p = 3\), sweep \(q \in \{17, 31, 61, 127\}\) and Monte-Carlo-estimate the decryption failure rate over 500 random keys × 1 plaintext each.

Theory. §48.2 — NTRU plus Ex 48.1 of the chapter.

Why this matters. Decryption-failure rate (DFR) is the correctness parameter of NTRU and Module-LWE schemes; it is exactly the input to the D’Anvers et al. (2019) attack you saw in W4.

# Bring in the NTRU toy from Chapter 48 §48.2.  For brevity we duplicate the
# polynomial helpers here (in the real notebook you can import them).

def cyclic_mul(a, b, N):
    c = [0] * N
    for i in range(N):
        for j in range(N):
            c[(i + j) % N] += a[i] * b[j]
    return c

def vec_mod(v, m): return [((x % m) + m) % m for x in v]
def vec_centred(v, m):
    v = vec_mod(v, m)
    return [x - m if x > m // 2 else x for x in v]


# (1) Polynomial inversion in Z_q[x] / (x^N - 1) via extended Euclidean
# (chapter 48 ships a pure-Python implementation).  Re-use that from the
# main chapter -- copy it across or 'from <your-module> import poly_inv'.

# Stub: a NumPy-free poly_inv is provided in the solution cell below.

def estimate_decryption_failure_rate(N, p, q, trials=200, rng=None):
    rng = rng or _r.Random(20260504)
    fails = 0
    for _ in range(trials):
        # TODO step 1: sample short ternary f, g until f is invertible mod p and q.
        # TODO step 2: compute h = p * F_q * g mod q.
        # TODO step 3: encrypt a fresh random m with fresh random r.
        # TODO step 4: decrypt; if recovered != m, increment `fails`.
        raise NotImplementedError('your turn')
    return fails / trials


# for q in (17, 31, 61, 127):
#     dfr = estimate_decryption_failure_rate(11, 3, q, trials=200)
#     print(f'q = {q:3d}   DFR ~ {dfr:.3%}')

Show code cell content

Hide code cell content

# Solution.
#
# The pattern matches §48.2 of the chapter exactly, plus a counter that
# increments each time the centred lift overshoots and decryption returns
# the wrong polynomial.  Expected behaviour: DFR drops sharply as q grows
# (because the noise term p*g*r + f*m has bounded coefficients, and once
# q comfortably exceeds twice that bound, the centred lift always works).

# Reuse the polynomial extended Euclidean from the main chapter ch48 code.
def _strip(p):
    while len(p) > 1 and p[-1] == 0: p.pop()
    return p
def _pmul(a, b, m):
    r = [0] * (len(a) + len(b) - 1)
    for i, ai in enumerate(a):
        if ai:
            for j, bj in enumerate(b):
                r[i+j] = (r[i+j] + ai*bj) % m
    return _strip(r)
def _psub(a, b, m):
    r = [0] * max(len(a), len(b))
    for i, c in enumerate(a): r[i] = (r[i] + c) % m
    for i, c in enumerate(b): r[i] = (r[i] - c) % m
    return _strip(r)
def _intinv(a, m):
    a %= m
    if a == 0: return None
    try: return pow(a, -1, m)
    except ValueError: return None
def _pdivmod(a, b, m):
    a = list(a); b = _strip(list(b))
    if len(a) < len(b): return [0], a
    invlc = _intinv(b[-1], m)
    q = [0] * (len(a) - len(b) + 1)
    while len(a) - 1 >= len(b) - 1 and any(c % m for c in a):
        d = len(a) - 1 - (len(b) - 1)
        coef = (a[-1] * invlc) % m
        q[d] = coef
        for i, bi in enumerate(b):
            a[d + i] = (a[d + i] - coef * bi) % m
        a = _strip(a)
        if len(a) - 1 < len(b) - 1: break
    return q, a
def poly_inv(f, mod, N):
    fc = [int(c) % mod for c in f]
    while fc and fc[-1] == 0: fc.pop()
    if not fc: return None
    g = [(-1) % mod] + [0] * (N - 1) + [1]
    r0, r1 = list(g), list(fc); s0, s1 = [0], [1]
    safety = 0
    while any(c % mod for c in r1) and safety < 200:
        q, r2 = _pdivmod(r0, r1, mod)
        s2 = _psub(s0, _pmul(q, s1, mod), mod)
        r0, r1 = r1, r2; s0, s1 = s1, s2
        safety += 1
    if len(r0) != 1 or r0[0] == 0: return None
    invl = _intinv(r0[0], mod)
    out = [(c * invl) % mod for c in s0]
    out += [0] * (N - len(out))
    return out[:N]


def estimate_decryption_failure_rate(N, p, q, trials=200, rng=None):
    rng = rng or _r.Random(20260504)
    fails = 0
    for _ in range(trials):
        # find invertible f
        for _ in range(2000):
            f = [1] * 4 + [-1] * 3 + [0] * (N - 7); rng.shuffle(f)
            g = [1] * 3 + [-1] * 3 + [0] * (N - 6); rng.shuffle(g)
            Fp = poly_inv(f, p, N); Fq = poly_inv(f, q, N)
            if Fp and Fq: break
        h = vec_mod([p * x for x in cyclic_mul(Fq, g, N)], q)
        m = [1] * 2 + [-1] * 2 + [0] * (N - 4); rng.shuffle(m)
        r = [1] * 3 + [-1] * 3 + [0] * (N - 6); rng.shuffle(r)
        c = vec_mod([x + mi for x, mi in zip(cyclic_mul(h, r, N), m)], q)
        a = vec_centred(cyclic_mul(f, c, N), q)
        b = vec_mod(a, p)
        recov = vec_centred([x % p for x in cyclic_mul(Fp, b, N)], p)
        if recov != m:
            fails += 1
    return fails / trials

for q in (17, 31, 61, 127):
    dfr = estimate_decryption_failure_rate(11, 3, q, trials=120)
    print(f'q = {q:3d}   DFR ~ {dfr:.3%}')
print('E3 OK')

q =  17   DFR ~ 87.500%
q =  31   DFR ~ 2.500%
q =  61   DFR ~ 0.000%
q = 127   DFR ~ 0.000%
E3 OK

Exercise 48.E4 — ★★★ UOV signing — fill in the linear solve

Goal. Complete the UOV signer. The setup gives you central-polynomial matrices \(\alpha_k\) (vinegar²) and \(\beta_k\) (vinegar × oil); you must build the \(o \times o\) system in the \(o\) oil unknowns and solve it.

Theory. §48.5 — Oil and Vinegar: UOV and the toy code in §48.5 of the main chapter.

Q = 31
V, O = 6, 4

def gauss_solve(M, b, q):
    n = len(b)
    A = [list(r) + [bi] for r, bi in zip(M, b)]
    for c in range(n):
        pi = next((r for r in range(c, n) if A[r][c] % q), None)
        if pi is None: return None
        A[c], A[pi] = A[pi], A[c]
        inv = pow(A[c][c], -1, q)
        A[c] = [(x * inv) % q for x in A[c]]
        for r in range(n):
            if r != c and A[r][c] % q:
                f = A[r][c]
                A[r] = [(A[r][k] - f * A[c][k]) % q for k in range(n + 1)]
    return [r[-1] % q for r in A]


def uov_sign(alpha, beta, target):
    '''alpha[k]: V x V vinegar^2 matrix.  beta[k]: V x O vinegar*oil matrix.

    target[k]: the message-hash residue we want f_k(v, o) to equal in GF(Q).
    Returns a (V+O)-length list (vinegar || oil), or None on bad luck.
    '''
    rng = _r.Random(2026)
    for _ in range(20):
        v = [rng.randrange(Q) for _ in range(V)]
        # TODO step 1: compute the o-by-o coefficient matrix M whose entry
        #              M[k][j] = sum_i beta[k][i][j] * v[i] (mod Q).
        # TODO step 2: compute the right-hand side
        #              b[k] = target[k] - sum_{i,j} alpha[k][i][j] * v[i] * v[j] (mod Q).
        # TODO step 3: solve M @ o = b in GF(Q).  If singular, retry with fresh v.
        # Return v + o.
        raise NotImplementedError('your turn')
    return None


# Test on a random instance (with a random target).
# rng = _r.Random(99)
# alpha = [[[rng.randrange(Q) for _ in range(V)] for _ in range(V)] for _ in range(O)]
# beta  = [[[rng.randrange(Q) for _ in range(O)] for _ in range(V)] for _ in range(O)]
# target = [rng.randrange(Q) for _ in range(O)]
# x = uov_sign(alpha, beta, target)
# # Verify.
# def f_k(k, x):
#     v, o = x[:V], x[V:]
#     val = 0
#     for i in range(V):
#         for jp in range(V): val = (val + alpha[k][i][jp] * v[i] * v[jp]) % Q
#         for jp in range(O): val = (val + beta[k][i][jp]  * v[i] * o[jp]) % Q
#     return val % Q
# assert all(f_k(k, x) == target[k] for k in range(O))
# print('E4 OK')

Show code cell content

Hide code cell content

# Solution.
#
# The Oil-Vinegar trick is exactly that the equations f_k(v, o) become
# AFFINE in o once v is fixed.  Splitting the constants and the o-coefficients
# across the equations gives an o-by-o linear system in GF(Q).  A random
# v makes that system non-singular with probability roughly 1 - O(1/Q), so
# we retry on bad luck.

def uov_sign(alpha, beta, target):
    rng = _r.Random(2026)
    for _ in range(20):
        v = [rng.randrange(Q) for _ in range(V)]
        M, b = [], []
        for k in range(O):
            row = [sum(beta[k][i][j] * v[i] for i in range(V)) % Q for j in range(O)]
            const = target[k]
            for i in range(V):
                for jp in range(V):
                    const = (const - alpha[k][i][jp] * v[i] * v[jp]) % Q
            M.append(row); b.append(const % Q)
        oil = gauss_solve(M, b, Q)
        if oil is not None:
            return v + oil
    return None

rng = _r.Random(99)
alpha = [[[rng.randrange(Q) for _ in range(V)] for _ in range(V)] for _ in range(O)]
beta  = [[[rng.randrange(Q) for _ in range(O)] for _ in range(V)] for _ in range(O)]
target = [rng.randrange(Q) for _ in range(O)]
x = uov_sign(alpha, beta, target)
def f_k(k, x):
    v, o = x[:V], x[V:]
    val = 0
    for i in range(V):
        for jp in range(V): val = (val + alpha[k][i][jp] * v[i] * v[jp]) % Q
        for jp in range(O): val = (val + beta[k][i][jp]  * v[i] * o[jp]) % Q
    return val % Q
assert all(f_k(k, x) == target[k] for k in range(O))
print('E4 OK')

E4 OK

Exercise 48.E5 — ★★★★ Toy MinRank modelling

Goal. Given two \(3 \times 3\) matrices \(M_1, M_2\) over \(\mathrm{GF}(7)\), find a non-trivial linear combination \(\lambda_1 M_1 + \lambda_2 M_2\) that has rank 1. This is the smallest-non-trivial MinRank instance behind the Beullens 2022 break of Rainbow.

Theory. §48.6 — Rainbow and the Beullens 2022 break.

# Hand-crafted instance: two 3x3 matrices over GF(7) such that some non-
# trivial combination has rank 1.  Construction (don't peek):
#   pick u = (1,2,3), v = (4,5,6); the outer product u v^T = R is rank 1.
#   pick M2 freely; set M1 = R - M2 mod 7.  Then 1*M1 + 1*M2 = R has rank 1.
# Your job: find ANY rank-1 combination.
M1 = [[1, 4, 6], [3, 1, 1], [4, 4, 1]]
M2 = [[3, 1, 0], [5, 2, 4], [1, 4, 3]]
P  = 7

def rank_mod(M, p):
    A = [r[:] for r in M]
    n = len(A); m = len(A[0]); rnk = 0
    col = 0
    for r in range(n):
        if col >= m: break
        pi = next((s for s in range(r, n) if A[s][col] % p), None)
        if pi is None:
            col += 1; continue
        A[r], A[pi] = A[pi], A[r]
        inv = pow(A[r][col], -1, p)
        A[r] = [(x * inv) % p for x in A[r]]
        for s in range(n):
            if s != r and A[s][col] % p:
                f = A[s][col]
                A[s] = [(A[s][k] - f * A[r][k]) % p for k in range(m)]
        rnk += 1; col += 1
    return rnk


def min_rank_brute(M1, M2, p, target_rank=1):
    '''Return (l1, l2) such that l1*M1 + l2*M2 has rank target_rank, or None.'''
    # TODO: enumerate (l1, l2) in GF(p) x GF(p) excluding (0, 0); compute the
    # combination and its rank.  Return the first match.
    raise NotImplementedError('your turn')


# coefs = min_rank_brute(M1, M2, P, target_rank=1)
# print('found combination:', coefs)
# print('E5 OK')

Show code cell content

Hide code cell content

# Solution.
#
# Brute-forcing 7^2 - 1 = 48 combinations is trivial; the real Beullens attack
# transforms this into a system of multivariate equations, then uses a
# Gröbner basis to solve it in time exponential in the matrix dimension but
# small for the parameters Rainbow used.  This toy lets you confirm by hand
# that the correct (3, 5) combination really has rank 1.

def min_rank_brute(M1, M2, p, target_rank=1):
    for l1 in range(p):
        for l2 in range(p):
            if l1 == 0 and l2 == 0: continue
            comb = [[(l1 * M1[i][j] + l2 * M2[i][j]) % p for j in range(3)] for i in range(3)]
            if rank_mod(comb, p) == target_rank:
                return l1, l2
    return None

coefs = min_rank_brute(M1, M2, P, target_rank=1)
print('found combination (l1, l2):', coefs)
assert coefs is not None, 'expected a rank-1 combination'
l1, l2 = coefs
comb = [[(l1 * M1[i][j] + l2 * M2[i][j]) % P for j in range(3)] for i in range(3)]
assert rank_mod(comb, P) == 1, 'reported combination is not rank 1'
print('combination matrix:', comb)
print('rank =', rank_mod(comb, P))
print('E5 OK')

found combination (l1, l2): (1, 1)
combination matrix: [[4, 5, 6], [1, 3, 5], [5, 1, 4]]
rank = 1
E5 OK

Exercise 48.E6 — ★★★★★ Research: walk a small supersingular isogeny graph

Goal. Pick a small prime \(p \equiv 3 \pmod 4\) (e.g. \(p = 71\)). Use SymPy or pure-Python elliptic-curve arithmetic over \(\mathrm{GF}(p^2)\) to:

1. List all supersingular \(j\)-invariants.

2. For each pair, decide whether they are connected by a 2-isogeny (via the modular polynomial \(\Phi_2(j_1, j_2) = 0\)).

3. Plot the resulting graph and confirm it is connected and \(\le \log_2 p\) in diameter.

Theory. §48.7 — Supersingular isogeny graphs.

Why this matters. The Castryck–Decru attack (§48.8) operates on this graph. Building it once for a small \(p\) — even if you cheat by using \(\Phi_2\) rather than running Vélu — gives you a concrete picture of what SIDH was protecting and what the attack glued.

Sage power-up

SageMath has EllipticCurve(GF(p^2), [...]) and supersingular-curve enumeration built in. If you have it (conda activate sage), this whole exercise is ~25 lines.

# (Open-ended research project.)  Skeleton:
#
# 1.  Define p = 71 and the modular polynomial Phi_2(X, Y) over Z[X, Y]:
#         Phi_2 = (X+Y)^3 - X^2 Y^2 + 1488 X Y (X + Y) - 162000 (X^2 + Y^2)
#                 + 40773375 X Y + 8748000000 (X + Y) - 157464000000000
#     (Source: any standard reference on modular polynomials.)
# 2.  Find supersingular j-invariants over GF(p^2) -- there are exactly
#     (p / 12) of them up to small corrections.  For p = 71, expect ~6 vertices.
# 3.  For each pair (j1, j2), evaluate Phi_2(j1, j2) mod p; non-zero -> no edge.
# 4.  Draw the graph (matplotlib + a tiny custom layout, or networkx if you
#     install it).

Show code cell content

Hide code cell content

# Reference solution (pure Python -- no sympy / Sage dependency).
#
# Two ingredients:
#  (a) Enumerate j-invariants of supersingular curves y^2 = x^3 + a x over GF(p)
#      via Euler's criterion on the cubic.
#  (b) Evaluate the (numerically large) classical modular polynomial Phi_2 in
#      INTEGER arithmetic, then reduce mod p.
#
# This sketch will produce a small (often disconnected when restricted to
# GF(p) and one curve family) sample of the full supersingular graph.  For
# the FULL graph over GF(p^2), use SageMath -- a one-liner with
# `EllipticCurve(GF(p^2), [a, b]).isogeny_graph(2)`.

def euler_criterion_sum(a, p):
    '''Return sum_{x in GF(p)} (x^3 + a x | p)_Legendre, computed lazily.'''
    s = 0
    for x in range(p):
        t = (x * x * x + a * x) % p
        if t == 0: continue
        s += pow(t, (p - 1) // 2, p)
        s %= p
    return s

def supersingular_j_invariants_montgomery(p):
    '''Curves y^2 = x^3 + a x for a in GF(p)*; supersingular iff trace(F) = 0.
    Returns the SET of distinct j-invariants found.  Note: this only sweeps
    one one-parameter family, so it under-samples the full supersingular set.'''
    js = set()
    for a in range(1, p):
        if euler_criterion_sum(a, p) == 0:
            # j-invariant of E_a: y^2 = x^3 + a x is 1728 (independent of a),
            # so this family yields a single j -- but DOES confirm
            # supersingularity for the family.
            js.add(1728 % p)
    return js

def Phi2_mod(j1, j2, p):
    '''Classical modular polynomial Phi_2(X,Y) reduced mod p.  Coefficients
    from any standard reference (e.g. Brökers 2009 Table 1).'''
    X, Y = j1, j2
    val = (X**3 + Y**3 - X**2 * Y**2
           + 1488 * X * Y * (X + Y)
           - 162000 * (X**2 + Y**2)
           + 40773375 * X * Y
           + 8748000000 * (X + Y)
           - 157464000000000)
    return val % p

p = 71
Js = supersingular_j_invariants_montgomery(p)
print(f'Sampled supersingular j-invariants over GF({p}):', sorted(Js))

# Sweep ALL j in GF(p), keep those that lie on the supersingular set
# (defined by Phi_2(j, j') = 0 having a root j' with j' supersingular).
# For p = 71 the full supersingular set has ceil(p/12) + ... vertices.
all_supersingular = set()
for j in range(p):
    # j is supersingular iff Phi_l(j, *) splits in a particular way -- for
    # this small toy we just check connectivity to a known supersingular j.
    if any(Phi2_mod(j, jp, p) == 0 for jp in Js):
        all_supersingular.add(j)
all_supersingular |= Js

edges = []
ss = sorted(all_supersingular)
for j1 in ss:
    for j2 in ss:
        if j1 < j2 and Phi2_mod(j1, j2, p) == 0:
            edges.append((j1, j2))

print(f'Recovered {len(ss)} candidate supersingular vertices:', ss)
print(f'2-isogeny edges (only between sampled vertices): {edges}')
print('E6 OK -- for the COMPLETE graph over GF(p^2), use Sage.')

Sampled supersingular j-invariants over GF(71): [24]
Recovered 2 candidate supersingular vertices: [17, 24]
2-isogeny edges (only between sampled vertices): [(17, 24)]
E6 OK -- for the COMPLETE graph over GF(p^2), use Sage.