Skip to main content

Hints — Module 16: The Spectral Test in Higher Dimensions

Graduated hints, three per stage, gentlest first. Read only as far as you need; after a stage is green, read WALKTHROUGH.md.

Stage 1: dual_basis, primal_basis, check_duality​

  1. This is §2's dual pair transcribed. Lemma D says the dual lattice is generated by a lower-triangular matrix with diagonal (m, 1, …, 1), and the anchor identity U·Vᵀ = m·I is exactly what check_duality must verify — so the two builders are not independent, they are a matched pair.
  2. Build each row keeping the powers of a reduced mod m as you go (pow = pow * a % m), because the tests pin the reduced entries (RANDU's third dual row is (−393225, 0, 1), not (−4295360521, 0, 1)). Dual row i (i ≥ 1) is −(a^{i} mod m) in column 0 and a 1 in column i; the primal first row is (1, a, a², … mod m) and every other primal row is m·e_i. check_duality computes all t² dot products and must reject mismatched shapes without panicking.
  3. v[0][0] = m; for i in 1..t { pow = pow*a%m; v[i][0] = -pow; v[i][i] = 1; } and u[0][j] = pow_j; u[i][i] = m for i ≥ 1. For check_duality, first guard v.len()==t and every row length == t returning false, then check dot(&u[i], &v[j]) == if i==j { m } else { 0 } for all i, j.

Stage 2: reduce_basis​

  1. This is Algorithm R: sweep every ordered pair to a fixpoint. Two facts make it correct — the matrix identity F·Eᵀ = I keeps U·Vᵀ = m·I alive through every row operation, and Σ_i |v_i|² is a positive integer that strictly decreases on every transformation, so the sweep must halt.
  2. For each ordered pair (i, j) with i ≠ j: let n = v_j·v_j, d = v_i·v_j; if 2|d| ≤ n the pair is already reduced — skip it (the strict trigger is load-bearing; transforming on the tie 2|d| = n loops forever). Else set q = round(d/n) and apply v_i ← v_i − q·v_j paired with the mirror u_j ← u_j + q·u_i — indices crossed, sign plus. Wrap all pairs in a sweep that repeats until a full pass changes nothing.
  3. The nearest-integer quotient for n > 0 is (2*d + n).div_euclid(2*n). The coupled update is one loop: for k in 0..t { v[i][k] -= q*v[j][k]; u[j][k] += q*u[i][k]; }. Structure it loop { let mut changed=false; for j { for i { … changed=true } }; if !changed { return } }.

Stage 3: shortest_vector_squared​

  1. This is Algorithm E, and its heart is the certificate: Cramer's rule (x·V = w ⇒ x_i·det V = w·g_i) plus Cauchy–Schwarz bounds each coefficient by |x_i| ≤ √(best₀)·|g_i|/|d|, so every vector no longer than the seed lives in a finite box. No orthogonality is assumed — the bound is rigorous for any nonsingular basis; reduction only shrinks it.
  2. Seed best with the smallest row norm. Compute det V exactly by Bareiss fraction-free elimination (pivot in a nonzero row flipping the sign, return 0 if a column is all zeros). Get the adjugate column g_i as the cofactors (sign·minor). Set z_i = ⌊√(best·|g_i|²/d²)⌋, then scan every x ≠ 0 with |x_i| ≤ z_i, taking the min of |x·V|². Skip x = 0, and you may force the first nonzero coefficient positive (w and −w are equal length).
  3. Bareiss step: m[i][j] = (m[i][j]*m[k][k] - m[i][k]*m[k][j]) / prev, then prev = m[k][k]. Cofactor: delete row i, column j, recurse into the determinant, multiply by (−1)^{i+j}. Scan recursively over the box keeping the partial sum Σ_{l<k} x_l·v_l incremental and updating best = best.min(norm) only at the leaf when some coefficient was nonzero.

Stage 4: nu_t_squared and mu_t​

  1. nu_t_squared is the whole Algorithm S pipeline — nothing new, just stages 1–3 composed. mu_t is Knuth's figure of merit: the volume of the t-ball of radius ν_t per lattice point, normalized by the exact half-integer Gamma values.
  2. nu_t_squared: assert 2 ≤ t ≤ 6, build the dual and primal bases, reduce_basis, then shortest_vector_squared on the reduced V. mu_t: μ_t = π^{t/2}·ν_t^t/(Γ(t/2+1)·m), with the exact table Γ(2)=1, Γ(5/2)=3√π/4, Γ(3)=2, Γ(7/2)=15√π/8, Γ(4)=6 — a wrong Gamma constant is exactly what the t = 4, 5, 6 tests catch.
  3. let nusq = nu_t_squared(a, m, t) as f64; then match t for gamma (3 => 0.75*PI.sqrt(), 5 => 1.875*PI.sqrt(), the rest rationals), and return PI.powf(t/2) * nusq.powf(t/2) / (gamma * m as f64) — note ν_t^t = (ν_t²)^{t/2}, so raise nusq to the half-power t/2.