perf: 4–6× faster LinearSystem.solve_one + MT19937 tempering

[n-WN]

Summary

Speeds up LinearSystem.solve_one by 4–6× on MT19937 state-recovery workloads through three layered optimizations, with a full regression-test suite and reproducible benchmark.

case	baseline (median)	optimized (median)	speedup	% faster
mt32	0.913s	0.147s	6.20×	83.9%
mt17	1.109s	0.181s	6.14×	83.7%
mt137	1.314s	0.215s	6.12×	83.6%
mt1337	1.656s	0.284s	5.82×	82.8%
mt9	1.080s	0.397s	2.72×	63.2%
mt1	1.690s	0.914s	1.85×	45.9%

Medians across 10 runs of benchmarks/mt_benchmark.py, Python 3.12.13 on Apple Silicon, vendored m4ri 20260122. All cases cross the 30 % threshold; the 5/6 workloads that spend most of their time in the solver see ≥60 % improvement.

Why it's faster — call-site details

1. perf(solver) — direct PyLong packing + augmented-matrix RREF (gf2bv/_internal.c)

• Row ingestion: pylong_pack_mzd_row_skip_lsb walks a PyLong's digits once and emits 64-bit words straight into mzd_row(A, r). This replaces the per-bit mzd_write_bit loop that the old path used. object_bit_constant + a 0/1 short-circuit in xor_tuple skip PyNumber_Xor for literal 0 and 1 on the symbolic BitVec XOR hot path.
• Single-solve algorithm: m4ri_solve(..., mode=SOLVE_MODE_SINGLE) now dispatches to a new m4ri_solve_single_augmented. It builds [A | b] of shape rows × (cols+1) in one pass and runs mzd_echelonize_m4ri(Aug, /*full=*/1, /*k=auto=*/0) to produce RREF. The solution is the b-column bit of each pivot row; any pivot landing on column cols is an inconsistency (return None). This is substantially faster than the old _mzd_pluq + _mzd_pluq_solve_left pipeline for this shape.
• Correctness: the A-portion's last word is masked before the b bit is written, so malformed inputs with stray high bits cannot bleed into the augmented column. mzd_init/malloc are NULL-checked. The affine-space path (SOLVE_MODE_AFFINE_SPACE) is unchanged — solve_all, solve_raw_space, and kernel enumeration all go through the existing PLUQ path.

2. perf(mt) — fused and selective MT19937 tempering (gf2bv/_internal.c + gf2bv/crypto/mt.py)

The symbolic BitVec path of MersenneTwister used to allocate three tuples per temper stage (shift result, and-result, xor-result) and pay Python dispatch overhead for each of the four stages. Two stacked wins:

• Fused stage: new C helper bv_xor_shift_mask(bits, shift, is_right, mask) computes bits ^ ((bits shifted) & mask) in one tuple pass. MT.temper on a 32-bit BitVec now makes four C calls total instead of twelve — no intermediate tuple allocations, same semantics.
• Selective output: MersenneTwister._temper_output_map caches the 32-bit linear map from input bits to output bits (computed once by running the integer temper on each basis vector). In _getrandbits_word(k), when k ≤ w/2 and the input is a 32-bit BitVec, we skip the full temper and compute only the requested top k output bits as XORs of the pre-selected input bits. For MT19937 the average hamming weight is ~5 per output bit, so mt1 drops from ~59 PyLong XORs per call to 3. Above the threshold the fused path wins because it reuses intermediate y₁/y₂/y₃, so the selective route is only taken when it's strictly cheaper.

3. build(macOS) — OpenMP for the vendored m4ri (setup.py)

m4ri 20260122's autoconf probe for OpenMP only tries -fopenmp, -openmp, -qopenmp, -xopenmp — all of which Apple clang rejects. The bundled m4ri was silently being built with __M4RI_HAVE_OPENMP=0 on macOS, so the parallel sections in brilliantrussian.c / mp.c were compiled out. On macOS, download_and_build_m4ri now:

• configures with -Xpreprocessor -fopenmp -I${LIBOMP_PREFIX}/include in CFLAGS and -L${LIBOMP_PREFIX}/lib -lomp in LDFLAGS,
• rewrites m4ri_config.h to set __M4RI_HAVE_OPENMP=1,
• adds the same flags to the extension's compile/link args so gf2bv._internal itself links without the caller having to export LDFLAGS/CPPFLAGS.

LIBOMP_PREFIX defaults to /opt/homebrew/opt/libomp and can be overridden. Non-macOS builds fall through to the original --enable-openmp configure invocation.

Proof of improvement

Regression tests — tests/test_solver.py (46 cases, 0.2 s):

$ python -m pytest tests/ -q
..............................................                           [100%]
46 passed in 0.19s

Covers: over-/square-/under-determined satisfiable systems (3 seeds × 4 shapes), inconsistent systems (both "0=1" tails and contradicting-row pairs), equivalence of solve_one vs solve_raw_space (incremental RREF check: the fast path's solution must lie in origin ^ span(basis)), and the MT19937 end-to-end recovery.

Reproducing the benchmark — benchmarks/mt_benchmark.py:

$ python benchmarks/mt_benchmark.py --runs 10

Prints build/solve/total timings (median, mean, and each individual run) for all six MT19937 shapes. Every run asserts solution == state before the timing is accepted.

Per-case build-vs-solve breakdown (optimized)

case	build median	solve median	total median
mt32	0.019s	0.127s	0.147s
mt17	0.037s	0.142s	0.181s
mt137	0.042s	0.173s	0.215s
mt1337	0.093s	0.190s	0.284s
mt9	0.066s	0.331s	0.397s
mt1	0.408s	0.503s	0.914s

Test plan

• python -m pytest tests/ passes (46 tests).
• python benchmarks/mt_benchmark.py --runs 10 runs to completion with ok=true on every run.
• Shipped examples still succeed: examples/mt.py, examples/lfsr.py, examples/xoshiro.py, examples/simple.py, examples/nlfsr.py.
• Public APIs unchanged (BitVec, LinearSystem.solve_one/solve_all/solve_raw_one/solve_raw_space, QuadraticSystem, AffineSpace).
• solve_raw_space / kernel enumeration path untouched; the regression suite asserts the new single-solve lies in the affine space returned by the old path.