Conversation
Add a Chinese Remainder Theorem based Smith Normal Form that avoids the recursive divide-and-conquer pipeline entirely. Factor N = prod p^e (passed in, so factoring is amortized across calls), compute the SNF over each local ring Z/p^e via valuation-pivoted Gaussian elimination, then CRT-recombine the per-prime unimodular transforms into U, V over Z/NZ. The local elimination is flat (no recursion), produces the unimodular transforms natively, and yields the divisibility chain for free. Exposed via the new rust_crt_snf PyO3 entry point alongside the existing recursive rust_smith_normal_form, which is left untouched for comparison. https://claude.ai/code/session_013PpJ2A2YvUCUk5EVtU5rse
Pure-Python reference implementation of the CRT-style SNF, a correctness harness validating it against the existing implementation as oracle, and a profiling script for head-to-head comparison against the recursive path. https://claude.ai/code/session_013PpJ2A2YvUCUk5EVtU5rse
Prior 0.2.0 release-hardening plan is fully shipped (v0.4.0: int64 contract, i128 intermediates, abi3 wheels, dual-backend tests). Replace with the research-backed roadmap for optimizing the CRT SNF path (Phases 0-4).
local_snf is 97-100% of CRT runtime (Phase 0 profiling). Its row/column clears reduced every update through posmod_i128. For prime-power moduli p^e < 2^31 the difference a-b*c fits in i64, so widen only when needed. Adds sub_mul_mod + i64 fast path in mulmod, and env-gated CRT_PROFILE phase-split instrumentation. ~1.5-2.05x faster across n=20..400, all moduli; 1.4-1.6x at n=1000-2000. Verified: invariant factors and U*A*V==S (mod N) unchanged vs dev baseline and optimize-modulus across 50 configs.
…se 2a) Restructure local_snf into Phase L (minimal-valuation column elimination -> U, upper-triangular mat) + Phase R (diagonalize -> V, reverse-order back- elimination to keep fill in unprocessed rows). This is the structural groundwork for a blocked GEMM trailing update, and already removes wasted work: the old inline row-clear updated all n rows of each column though only the pivot row was nonzero. ~1.1-1.3x faster than Phase 1. Adds a randomized cargo test (no external deps) asserting local_snf yields a valid Smith form: U*A*V == diag(p^vals) (mod q), ascending valuations, and U,V unimodular (det != 0 mod p), across 500+ cases. Verified: zero invariant- factor mismatches vs the independent dev baseline on the head-to-head sweep.
Capture the valuation-level-staging insight (global min-valuation pivoting defeats naive blocking; stage by p-adic level so each level is unit-pivot blocked LU) and the L21/U12 bookkeeping + deflation plan, so 2b can be executed cleanly against the cargo oracle.
Phase L now clears the valuation-0 bulk with blocked LU (panel width 48): in-place panel factorization, TRSM (L11^-1) for the deferred panel-row block, and a GEMM trailing update (mat[pend..,B..] -= L21@U12 mod q, same for U) with an i64 accumulator + delayed reduction (i128 fallback for large q). The scalar min-valuation path finishes the higher-valuation / rank-deficient tail; a final pass normalizes blocked pivots (left as units to keep mat and u row-consistent). Speedup grows with n (cache/blocking): ~1.3-1.5x over Phase 2a at n=400, and 1.6-2.6x over Phase 1 at n=1000-2000. Cumulative ~3x over the original CRT. Verified: two randomized cargo tests (600+ cases incl. multi-panel, rank- deficiency mod p forcing the blocked->scalar handoff, p|A skipping the blocked pass, rectangular) assert valid SNF; head-to-head shows zero invariant-factor mismatches vs the dev baseline and all transform_ok across the sweep.
The naive 2b trailing GEMM read panel rows target[k+s,col] with stride m across s -> ~pb distinct cache lines per output element. Pack U12 (post-TRSM panel rows) once into a contiguous (pb x ncols) buffer, then update each trailing row with contiguous axpy passes (autovectorizes); i128 accumulator fallback for large moduli. Tried an f64 GEMM via ndarray/matrixmultiply first: helped single-prime moduli but REGRESSED multi-prime (N=360: 0.86x) due to per-panel n*n product allocation + thin-K (pb=48). The packed i64 kernel wins across both: ~1.34x (n=1000,N=2), 1.24x (n=1000,N=360), 1.10x (n=2000) over Phase 2b. Correctness unchanged (cargo tests + head-to-head vs dev: no mismatches, all transform_ok).
Tall matrices (n > m) were padded to square and had V cropped, which could leave V non-unimodular for column-rank-deficient inputs (an all-zero column gave V = [[0]]). U A V == S still held, so the old tests missed it. Route tall matrices through the transpose, where the wide path crops U (which stays unimodular) instead of V. Add a Rust correctness oracle for smith_normal_form (U A V == S, divisibility chain, unimodularity over random inputs) plus a Storjohann-vs-CRT invariant-factor cross-check, which is what surfaced the bug. Tests use rand's seeded StdRng (dev-dependency).
- Drop the CRT_PROFILE timing instrumentation and phase-by-phase development comments from crt.rs; keep the algorithmic invariants (i64 overflow bounds, deferred normalization). - Replace section-divider comments and TRSM/GEMM shorthand with prose matching the rest of the crate; rustfmt ring.rs/snf.rs. - Seed crt.rs tests with rand's StdRng instead of a hand-rolled LCG. - Replace the dev roadmap (docs/PLAN.md) and the CRT prototype/profiling scripts with a concise design note (docs/crt.md).
The pure-Python Storjohann implementation is removed; modularsnf is now a thin wrapper over the modularsnf._rust extension. Validation and marshalling live in _core.py; snf.py and crt.py just call the binding. The public API is a single entrance -- smith_normal_form_mod, crt_snf, SNFResult -- with RingZModN/RingMatrix/diagonal dropped. The PyO3 wrapper now delegates to modularsnf::smith_normal_form (which also propagates the tall-matrix fix to Python) and exposes only rust_smith_normal_form and rust_crt_snf; the dead RustRingZModN, diagonal, and merge exports are removed.
Algorithm correctness now lives in the Rust crate. The Python suite (test_binding.py) covers only the binding: input validation, type/shape marshalling, and one end-to-end SymPy oracle on the invariant factors. Drop the old correctness/internal suites and the python/rust backend switch (conftest, helpers, pyproject perf marker, justfile).
Storjohann band reduction stays the documented default; add a CRT fast-path section pointing at docs/crt.md. Update README/AGENTS to describe the Rust-crate-plus-thin-binding layout and that correctness is validated in the Rust suite.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Add an experimental, more performant CRT-based SNF algorithm: factor the modulus and solve a local SNF over each prime power via valuation-pivoted blocked LU, then CRT-recombine.
Changelog