Add experimental CRT-based SNF fast path

AGENTS.md

@@ -34,25 +34,27 @@ Use `modularsnf.ring` primitives: `gcdex`, `div`, `quo`, `stab`, `ann`.
 
 ## Testing
 
-* Structural assertions (`S = U A V`, diagonal shape, divisibility chain,
-  unimodularity) over hardcoded array equality.
-* Validate against SymPy integer-domain projection where feasible.
-* Run `just check` (or individually: `just lint`, `just typecheck`, `just test`) before finalizing.
+* Algorithm correctness lives in the Rust crate (`cargo test`): structural
+  assertions (`S = U A V`, diagonal shape, divisibility chain, unimodularity)
+  over random inputs, plus a Storjohann-vs-CRT cross-check.
+* The Python suite (`tests/`) covers only the PyO3 boundary: input validation,
+  type/shape marshalling, and the public-API contract on small cases.
+* Run `cargo test` and `just check` (`just lint`, `just typecheck`, `just test`)
+  before finalizing.
 
 ## File Structure
 
-* `modularsnf/ring.py` — ring primitives (Gcdex, Stab, Div, Ann).
+The SNF algorithms live in the Rust workspace; the Python package is a thin
+wrapper over the `modularsnf._rust` extension.
+
+* `crates/modularsnf/` — Rust lib crate: ring primitives, echelon/band
+  reduction, diagonalization (`snf.rs`/`band.rs`/`echelon.rs`/`diagonal.rs`),
+  and the CRT fast path (`crt.rs`).
+* `crates/modularsnf-py/` — PyO3 bindings exposed as `modularsnf._rust`.
+* `modularsnf/ring.py` — `RingZModN`, forwarding ring primitives to Rust.
 * `modularsnf/matrix.py` — `RingMatrix` data structure.
-* `modularsnf/echelon.py` — triangularization (Lemma 3.1).
-* `modularsnf/band.py` — band reduction (Lemmas 7.3, 7.4, Prop 7.1).
-* `modularsnf/diagonal.py` — diagonalization (Prop 7.7, Theorem 7.11).
-* `modularsnf/snf.py` — master pipeline (Lemma 7.14) + public
-  `smith_normal_form_mod` API.
-* `docs/algorithm.md` — mathematical foundation and algorithm description.
-* `docs/PLAN.md` — current ExecPlan.
-
-## ExecPlans
-
-For work spanning multiple modules or changing public APIs, write a plan in
-`docs/PLAN.md` before implementing. Users may say "use an ExecPlan" as
-shorthand.
+* `modularsnf/diagonal.py` — diagonalization wrappers.
+* `modularsnf/snf.py` — public `smith_normal_form_mod` API (default path).
+* `modularsnf/crt.py` — public `crt_snf` API (experimental fast path).
+* `docs/algorithm.md` — default algorithm (Storjohann band reduction).
+* `docs/crt.md` — CRT fast path design and correctness basis.

README.md

@@ -1,10 +1,10 @@
 # Smith Normal form of Integer matrices mod N (Storjohann)
 
-This is a Python module that follows the deterministic algorithms presented in Arne Storjohann's PhD Dissertation *Algorithms for Matrix Canonical Forms* (ETH No. 13922, 2000).
+A Rust implementation (with a thin Python binding) of the deterministic algorithms presented in Arne Storjohann's PhD Dissertation *Algorithms for Matrix Canonical Forms* (ETH No. 13922, 2000).
 
-It implements the Lemmas and subsequent subroutines that are necessary for calculating the SNF without exponential intermediate values.
+It implements the Lemmas and subsequent subroutines that are necessary for calculating the SNF without exponential intermediate values. The algorithms live in the Rust `modularsnf` crate; the Python package is a thin wrapper over the native extension.
 
-It validates against SymPy using a known equivalence between calculating the Smith Normal form of an integer matrix, and then taking mod N, compared to solving it natively in the ring.
+Correctness is validated in the Rust test suite: random inputs are checked against the Smith form contract ($S = UAV$, divisibility chain, unimodular transforms), and the default Storjohann path is cross-checked against the independent CRT path.
 
 ## Quick Start
 
@@ -28,22 +28,30 @@ A matrix is **unimodular** over $\mathbb{Z}/N\mathbb{Z}$ when
 $\gcd(\det(M),\, N) = 1$, the modular analogue of $|\det(U)| = 1$
 over $\mathbb{Z}$.
 
-### Lower-level API
+For details on the default algorithm (band reduction, diagonalization,
+Storjohann's lemmas), see [docs/algorithm.md](docs/algorithm.md).
 
-For direct access to `RingMatrix` objects:
+## Alternative: CRT fast path (experimental)
+
+`smith_normal_form_mod` uses Storjohann's band reduction, which works for any
+modulus `N >= 2`. For the **small-modulus, large-matrix** regime there is a
+second, experimental algorithm that is several times faster: a CRT-based path
+that factors `N = prod p^e`, solves the SNF over each local ring `Z/p^e` by
+valuation-pivoted elimination, and recombines via the Chinese Remainder Theorem.
 
 ```python
-from modularsnf.ring import RingZModN
-from modularsnf.matrix import RingMatrix
-from modularsnf.snf import smith_normal_form
+from modularsnf import crt_snf
 
-ring = RingZModN(12)
-A = RingMatrix(ring, [[2, 4], [6, 8]])
-U, V, S = smith_normal_form(A)   # note: (U, V, S) order
+S, U, V = crt_snf([[2, 4, 0],
+                   [6, 8, 3],
+                   [0, 3, 9]], modulus=36)
+# Same (S, U, V) contract as smith_normal_form_mod; S = U @ A @ V (mod 36).
 ```
 
-For details on the algorithm (band reduction, diagonalization, Storjohann's
-lemmas), see [docs/algorithm.md](docs/algorithm.md).
+It requires the prime factorization of `N`; pass it explicitly as
+`factors=[(p, e), ...]` to amortize factoring across many calls with the same
+modulus. Prefer `smith_normal_form_mod` for general moduli. See
+[docs/crt.md](docs/crt.md) for the design and correctness basis.
 
 ## Development Workflow (modern Python)