How I Made G6K’s Spherical Code Generator 100× Faster with Incremental Gram Matrix Updates in Rust

TL;DR

G6K — the leading open-source lattice sieving toolkit used in post-quantum cryptanalysis — relies on sparse spherical codes as the probe vectors for its SimHash-based pair filter. Those codes are produced by spherical_coding/codes.py, a single-threaded Python hill-climber that takes minutes per dimension and was effectively unusable across a wide dimension sweep.

I rewrote the generator in Rust around a single algorithmic observation: a coordinate swap inside one row of the code matrix only modifies one row of the Gram matrix. Treating that locality of change properly turns an O(N · n) update into O(N), and the cached row-sum-of-squares used by the global score turns an O(N²) score query into O(N). Combined with parallel candidate search via rayon, the new tool generates the full code bank for n = 6..1000 on a laptop in the time the original version takes to do a handful of dimensions.

Output files are byte-identical to the originals. No downstream G6K code changes required.

This is a write-up of internal work — the tool is not currently open-sourced. If you’re working on G6K-adjacent research and a faster code generator would unblock you, get in touch (contact at the bottom).

Why this matters: G6K, lattice sieving, and the trust we place in lattices

G6K (the General Sieve Kernel, Albrecht, Ducas, Herold, Kirshanova, Postlethwaite, Stevens & Wessel, Eurocrypt 2019) is the de facto open implementation of modern lattice sieving algorithms. Lattice sieving is, today, the most practical attack on the shortest vector problem (SVP) and its approximate variants — the problems whose hardness underpins NIST-standardized post-quantum schemes like Kyber and Dilithium.

In other words: when a lattice-based post-quantum scheme is benchmarked, profiled, or attacked, the work very often runs through G6K. Any improvement to G6K’s setup or runtime feeds directly into how cleanly we can stress-test the next generation of standards.

The sieving inner loop

A lattice sieve maintains a database of lattice vectors and repeatedly looks for pairs (v, w) such that v ± w is shorter than max(|v|, |w|). With a database of N vectors, naive pair search is O(N²), and N itself grows as roughly (4/3)^(n/2) in dimension n. By dimension 80, brute-force pair search is hopeless.

The breakthrough that made high-dimensional sieving practical is locality-sensitive filtering — most famously the BDGL filter (Becker, Ducas, Gama & Laarhoven, SODA 2016) and its descendants. The principle: cheaply reject most pairs without computing v ± w at all.

SimHash as G6K’s filter

G6K uses a SimHash-style filter (in the spirit of Charikar’s 2002 random-hyperplane locality-sensitive hash, adapted for sieving). The construction:

1. Fix a bank of probe vectors c_1, …, c_b ⊂ ℝⁿ (here b = 256 or 512).
2. For each lattice vector v, compute the bit sign(⟨c_i, v⟩) for every probe.
3. Pack those bits into a fingerprint.
4. To test a pair, XOR fingerprints and popcount: small Hamming distance ⇒ similar orientation ⇒ worth a real reduction attempt.

A single 64-bit XOR-popcount can reject pairs that would have cost a full O(n) inner product. Across a sieve, the savings are the difference between hours and weeks.

But the filter only works if the probe bank is well-distributed. Pick the wrong probes and you reject good pairs or accept noise — either way, sieving collapses.

Spherical codes: the unsung dependency

G6K’s probe banks are sparse spherical codes: each probe has exactly six non-zero entries (three at +1, three at −1), spread across the n coordinates. The sparsity is intentional — it keeps ⟨c_i, v⟩ cheap. The “spherical” part means the bank should approximate a uniform distribution on the sphere, so that the SimHash bits are nearly independent across probes.

These banks live in G6K as files named sc_<n>_<bitlen>.def — one per dimension, regenerated whenever the sieve targets a new n. They are produced by spherical_coding/codes.py, a Python hill-climber that:

1. Generates bitlen random sparse rows.
2. Builds the Gram matrix G = M Mᵀ (size bitlen × bitlen).
3. Defines a global score S(G) = Σᵢ (Σⱼ Gᵢⱼ² − G₀₀²)² measuring spread quality (lower is better).
4. Repeatedly tries swapping two coordinates inside a single row and keeps the swap if the score drops.

Conceptually clean. Operationally, painful.

The technical bottleneck (and what doesn’t fix it)

Profiling the original codes.py shows the inner loop is dominated by two operations, both of which sit inside the per-swap evaluation:

def update_G(M, G, i):
    v = M.dot(M[i].T)            # O(bitlen * n)
    for j in range(l):
        G[i, j] = v[j]
        G[j, i] = v[j]

def score(G):
    s = G[0, 0] * G[0, 0]
    return sum((sum(x*x for x in row) - s)**2 for row in G)  # O(bitlen²)

For bitlen = 512, n = 128, that’s roughly 65 K + 262 K floating-point ops per candidate swap, with n²/2 ≈ 8 K swaps per row pass and dozens of passes per code. The interpreter overhead on top makes the wall-clock time absurd.

What I tried first (and why it wasn’t enough)

The reflex move was to keep the algorithm and accelerate the language:

• Vectorized NumPy rewrites of score and update_G: ~3× faster.
• Numba JIT on the hot functions (@jit(nopython=True) on score, update_G, and the swap loop): another ~5–8×.
• Multiprocessing over candidate codes: linear in cores.

Stacked, this gets you maybe 30–50× over the original on a good day. Useful, but the per-dimension cost was still in the minutes, and crucially, the algorithmic cost was unchanged. Each candidate swap still triggered a full M.dot(M[i].T) and a full score(G). JIT only compiled the wasted work to be wasted faster.

The actual fix was upstream of any language choice.

The locality-of-change observation

A swap inside row i of M only changes the entries of G that involve row i:

• G[i, j] for all j (and symmetrically G[j, i]) — exactly two rows of the Gram matrix change, and they’re identical by symmetry.
• The other bitlen × (bitlen − 2) entries of G are bit-identical to before.
• The score S(G) therefore depends on the changed values only through the row sums of squares of rows i and the affected j’s.

Concretely, if a swap on row row changes M[row, col_a] by delta_a and M[row, col_b] by delta_b:

ΔG[row, j] = delta_a · M[j, col_a] + delta_b · M[j, col_b]   for all j ≠ row
ΔG[j, row] = ΔG[row, j]                                       (symmetry)
G[row, row] is invariant                                      (always equals SPARSITY = 6)

That’s O(bitlen) arithmetic per swap, not O(bitlen · n). And by maintaining a cached vector row_sq_sums[i] = Σⱼ G[i,j]², we can update the score with two scalar updates per changed entry instead of recomputing S from scratch. Score query collapses from O(bitlen²) to O(bitlen).

The Rust implementation

I packaged this into a small Rust crate, spherical-codegen, with three modules:

spherical_coding/spherical-codegen/
├── Cargo.toml
└── src/
    ├── code_matrix.rs     # SparseCodeMatrix: dual sparse + dense storage
    ├── gram.rs            # GramMatrix: incremental updates + cached row_sq_sums
    └── optimizer.rs       # hill-climber: exhaustive best-swap selection

Dual representation: sparse for I/O, dense for lookups

pub struct SparseCodeMatrix {
    pub n: usize,
    pub bitlen: usize,
    pub rows:  Vec<[u16; SPARSITY]>,   // 6 indices per row (3 +1, 3 −1)
    pub dense: Vec<i8>,                 // n × bitlen flat array, values in {-1, 0, 1}
}

rows is what lands in the .def file. dense gives O(1) lookups for M[j, col_a] during the incremental Gram update — without it, you’d be scanning sparse rows on every term.

Incremental Gram update

pub fn update_row_incremental(
    &mut self,
    code: &SparseCodeMatrix,
    row: usize,
    col_a: usize, col_b: usize,
    delta_a: i32, delta_b: i32,
) {
    for j in 0..self.size {
        if j == row { continue; }                       // diagonal invariant
        let mj_a = code.get(j, col_a) as i32;
        let mj_b = code.get(j, col_b) as i32;
        let delta = delta_a * mj_a + delta_b * mj_b;
        if delta != 0 {
            let old_val = self.data[row * self.size + j];
            let new_val = old_val + delta;
            let sq_delta = (new_val as i64).pow(2) - (old_val as i64).pow(2);
            self.row_sq_sums[row] += sq_delta;
            self.row_sq_sums[j]   += sq_delta;          // symmetric row also affected
            self.data[row * self.size + j] = new_val;
            self.data[j * self.size + row] = new_val;
        }
    }
}

The row_sq_sums cache is verified against a from-scratch recompute in test_incremental_update_matches_full — every entry of G and every cached row-sum-of-squares is asserted equal after a real swap. This test caught two off-by-one bugs during development; I would not deploy this without it.

Diagram: what changes per swap

            Before swap                 After swap on row r
            ┌─────────────┐             ┌─────────────┐
            │  G  bitlen  │             │  G  bitlen  │
            │             │             │             │
   row r ─► │ ░░░░░░░░░░░ │   ─────►    │ ▓▓▓▓▓▓▓▓▓▓▓ │  ◄─ recomputed
            │             │             │             │
            │ ░ ◄─ col r  │             │ ▓ ◄─ col r  │  ◄─ recomputed (symmetry)
            │ ░           │             │ ▓           │
            │ ░           │             │ ▓           │
            └─────────────┘             └─────────────┘
            261,632 entries unchanged.   1,022 entries change.

Old code recomputed all 262,144 entries via a fresh M @ M.T slice. New code touches the 1,022 entries that actually moved.

Best-swap selection vs first-improvement

The Python original returned on the first improving swap. With cheap incremental updates, that heuristic is no longer optimal: the per-swap evaluation cost has dropped enough that searching every pair (a, b) in a row and applying only the best-improving swap is a net win. It converges to lower local minima in fewer outer-loop passes.

for aa in 0..n {
    for bb in 0..aa {
        let (a, b) = (positions[aa], positions[bb]);
        if code.get(row, a) == code.get(row, b) { continue; }
        // tentative swap + incremental gram update
        // measure improvement
        // record best, then revert
    }
}
// apply best swap if improvement > 0

Parallel candidates and a batch driver

Hill climbers land in different basins from different seeds. Running multiple candidates in parallel via rayon and keeping the lowest score is essentially free given modern core counts:

let best = (0..args.candidates)
    .into_par_iter()
    .map(|i| {
        let mut rng = ChaCha8Rng::seed_from_u64(i as u64 ^ 0xdeadbeef);
        let mut code = SparseCodeMatrix::random(args.n, bitlen, &mut rng);
        let score = improve(&mut code, &config, &mut rng);
        (code, score)
    })
    .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap())
    .unwrap();

A run_batch.sh driver sweeps n = start..end, skips already-generated .def files (idempotent resume), and caps external parallelism since each candidate already saturates internal threads.

Release profile

LTO and codegen-units = 1 matter here. The hot path crosses three modules — optimizer calls gram.update_row_incremental which calls code.get — and inlining across crate-internal modules is what keeps the inner loop tight enough to make exhaustive best-swap selection viable.

[profile.release]
opt-level = 3
lto = true
codegen-units = 1

Results

On an M-series MacBook (8-core), bitlen = 512:

• n = 64: ~1.5 s per candidate, ~5 s for 3 candidates in parallel
• n = 128: ~3 s per candidate, ~10 s total
• n = 256: ~12 s per candidate, ~40 s total
• Full sweep n = 6..1000: completes overnight at worst, often in a single coffee break depending on time-limit setting

Score quality is strictly better across the dimensions I’ve sampled, primarily because of best-swap selection plus best-of-3 candidate search. Better codes mean cleaner SimHash filtering downstream — fewer wasted reductions per sieve pair-pass.

Prior work and where this fits

• G6K (Albrecht et al., Eurocrypt 2019) — the toolkit this work plugs into; the SimHash filter is a core component.
• BDGL filter (Becker, Ducas, Gama, Laarhoven, SODA 2016) — the locality-sensitive filter design G6K’s SimHash is a practical descendant of.
• BGJ1 / BGJ-style sieves — earlier locality-sensitive sieving constructions; G6K supersedes them but the filter idea is shared lineage.
• Charikar 2002 — the random-hyperplane SimHash, original LSH for cosine similarity.
• Spherical code construction — the hill-climbing approach in codes.py is a practical heuristic; it’s not derived from any specific paper I’m aware of, just the natural greedy approach to the spread objective. If you know of a more principled construction (e.g., algebraic spherical codes adapted to the sparse-{−1, 0, +1} alphabet), I’d genuinely like a pointer.

This work doesn’t change the algorithm or improve the codes’ theoretical properties — it’s an engineering optimization. The contribution is making code regeneration cheap enough that you can experiment with code parameters (bitlen, sparsity, alphabet) instead of treating the bank as a frozen artifact.

Status and availability

The tool is internal. It produces drop-in .def files in the exact format G6K’s C++ side expects via kernel/simhash.inl, runs reliably across the full n = 6..1000 range at bitlen = 256 and 512, and has a property-test suite verifying that incremental Gram updates and the cached row_sq_sums match a from-scratch recompute on every entry.

I haven’t open-sourced it for two reasons:

1. Quality bar. The hill-climber is good but not provably optimal — there’s likely another order of magnitude available with simulated annealing or a proper SDP relaxation, and I’d rather release once than release twice.
2. Time. Maintenance, issues, packaging, integration patches against upstream G6K — that’s its own commitment.

If you’re working on G6K, post-quantum benchmarking, BKZ + sieve hybrids, or any project where a faster spherical code generator would unblock real work, I’m happy to share the binary, run codes for specific dimensions on request, or discuss a more formal collaboration.