feat(cpp): ✨ select the support-form row store per propagator

AGENTS.md

@@ -96,11 +96,43 @@ Key files:
 - **`Monomial<N>`** (`cpp/monoprop/core/Monomial.h`) = `Bitset<2*N>`: ONE basis operator, two bits per
   mode/qubit. Basis-agnostic — read as a Majorana product, or as a Pauli string (JW image).
   Collections: `MonomialList<N>` (no coeffs) and `MonomialMap<N>` (monomial → real coeff).
-- **Row access** (`cpp/monoprop/detail/operator/RowAccess.h`): the one backend-agnostic vocabulary
-  (`materialize_row`, `assign_row`, `row_popcount`, `for_each_row_position`) over the dense
-  `MonomialList<N>` and the packed `detail::OperatorIndex<N>`. Any template parameterized on the row
-  store must include that header — the `OperatorIndex` overloads live in `monoprop::`, so ADL cannot
-  find them from a `monoprop::detail` argument.
+- **The row-store seam** (`cpp/monoprop/detail/operator/RowAccess.h`): a dense monomial is a transient,
+  not the storage. Four accessors — `materialize_row`, `assign_row`, `row_popcount`,
+  `for_each_row_position` — and three backends answer them: `MonomialList<N>`, `detail::OperatorIndex<N>`
+  (packed position lists) and `detail::SparseRowStore<N>` (fixed-width mode lanes plus one 2-bit-per-slot
+  `codes` word per row). Reach rows through the accessors, never through a backend's own API, and add any
+  fourth backend to `cpp/tests/row_accessor_tests.cpp`, which asserts that all of them agree through every
+  accessor. Any template parameterized on the row store must include that header — the overloads live in
+  `monoprop::`, so ADL cannot find them from a `monoprop::detail` argument.
+  `algebra/CodesAlgebra.h` is the structural algebra on a sparse row, one function per dense counterpart,
+  reading the `codes` word instead of looping over storage words, plus `sparse_toggle` — the product
+  `M ⊕ G` as one merge over two ascending lane arrays. It is exact, not an approximation:
+  `cpp/tests/codes_algebra_tests.cpp` and `codes_product_tests.cpp` assert agreement with each dense
+  version over the fixtures and randomized rows. Change one side and you must change the other. A product
+  can occupy more modes than either input; past its scratch capacity `sparse_toggle` reports `overflowed`
+  and the caller must fall back to the dense product — never truncate, because a truncated mode list still
+  carries a plausible-looking `codes` word.
+- **Which backend, and where it is bound**: a propagator uses one of the two row stores, chosen once from
+  its mode count by `SparseRowStore::preferred_for_modes()` — a build-time constant
+  (`monoprop_SPARSE_ROW_MIN_MODES`, derived in `CMakeLists.txt` from whether `ARCH_FLAG` is actually
+  emitted rather than from the option that asks for it, and deliberately not a cache entry) because what
+  moves the crossover is the target ISA. `monoprop_ROW_STORE=dense|sparse` forces it process-wide; an
+  unrecognized value throws rather than falling back, since the point of setting it is to know which
+  backend ran. `MPOperator` holds one pointer per backend with exactly one non-null and binds the live one
+  via `with_store` — **once per layer, inside `build_layer`**, never per term: the scan asks the store for
+  a row per anticommuting term, so everything downstream is templated on the store
+  (`LayerBuildEngine<N, Sink, Store>`, `fused_find_and_collect`, `probe_incoming_queries`). Off that path,
+  use `MPOperator`'s forwarding accessors; there is no accessor handing out a store, because there is no
+  one type to hand out — which is why `MonomialPropagator` exposes `for_each_term()`/`num_local_terms()`
+  rather than the `indexing()` it used to. Every C++ case runs a second time under
+  `monoprop_ROW_STORE=sparse` (the `sparse-rows` ctest label) — every fixture is below the crossover, so
+  without that the sparse backend would ship untested; `cpp/tests/row_store_selection_tests.cpp` is what
+  fails if the variable stops reaching the propagator. The two backends agree on term sets and values but
+  not on term *order*, so compare them with `just diff-baseline-sparse` (tolerance), never
+  `just diff-baseline` (byte-wise). A benchmark run records the backend it resolved to in its artifact's
+  `meta` (`monoprop_row_store` as asked, `row_store_effective` as run) and `REPORT.md` shows both: under
+  the default `auto` the setting alone does not identify the backend, and the two differ in footprint and
+  in accumulation order.
 - **`Basis` / the `Algebra` policy** (`cpp/monoprop/algebra/`): the two algebras are sibling models
   (`MajoranaAlgebra`, `PauliAlgebra` in `algebra/Algebra.h`) over shared structural primitives
   (`algebra/AlgebraCommon.h`). The propagation backbone (the scan/fold in `detail/evolution/...`) is

CMakeLists.txt

@@ -77,6 +77,20 @@ endif()
 include(${PROJECT_SOURCE_DIR}/cmake/compiler_flags/Sanitizers.cmake)
 include(${PROJECT_SOURCE_DIR}/cmake/compiler_flags/CXXFlags.cmake)
 
+# The storage width at or above which a propagator picks the support-form row store over the dense one.
+# Derived from whether ARCH_FLAG is actually emitted rather than from the option that asks for it, and
+# deliberately not a cache entry: what moves the crossover is the target ISA, so a stale cached value
+# would silently pick the wrong backend after a flag change.
+#
+# Thresholds are the first full 32-mode block where sparse is clearly faster than dense beyond
+# run-to-run noise. Expect about +/-1 block variation across machines due to cache and popcount
+# throughput.
+if(ARCH_FLAG)
+  set(monoprop_SPARSE_ROW_MIN_MODES 768)
+else()
+  set(monoprop_SPARSE_ROW_MIN_MODES 256)
+endif()
+
 # report on compiler flags in use
 message(STATUS "Configuring a ${CMAKE_BUILD_TYPE} build")
 string(TOUPPER ${CMAKE_BUILD_TYPE} _cmake_build_type_upper)
@@ -92,6 +106,10 @@ message(
   "   Build-type-specific    : ${_cmake_build_type_specific_flags}"
 )
 message(STATUS "   Vectorization flag     : ${ARCH_FLAG}")
+message(
+  STATUS
+  "   Sparse rows from       : ${monoprop_SPARSE_ROW_MIN_MODES} modes"
+)
 message(
   STATUS
   "   Project defaults       : ${CMAKE_CXX${CMAKE_CXX_STANDARD}_STANDARD_COMPILE_OPTION} ${monoprop_CXX_FLAGS}"

README.md

@@ -117,8 +117,12 @@ uv sync --all-groups --all-extras -v    # installs the workspace, incl. the benc
 uv run python -m pytest -m "not mpi"   # Python tests (serial)
 just test-mpi                          # Python + C++ tests under MPI
 just test-wide                         # Python + C++ unit tests with a 64-bit TermIndex
+just test-sparse-rows                  # Python tests with the support-form row backend forced
 ```
 
+The C++ suite runs against both row backends: `ctest` registers every case a second
+time with `monoprop_ROW_STORE=sparse`, labelled `sparse-rows`.
+
 See the [testing guide](https://docs.monoprop.algorithmiq.tech/testing)
 for the with/without-MPI details and the rank matrix.
 

benches/conftest.py

@@ -214,6 +214,7 @@ def _meta(nodes: int, ranks_per_node: int) -> dict[str, Any]:
         "nodes": nodes,
         "ranks_per_node": ranks_per_node,
         "monoprop_threads": os.environ.get("monoprop_NUM_THREADS", "default"),  # noqa: SIM112
+        "monoprop_row_store": os.environ.get("monoprop_ROW_STORE") or "auto",  # noqa: SIM112
         "cpu_count_logical": psutil.cpu_count(logical=True),
         "cpu_count_physical": psutil.cpu_count(logical=False),
         "hostname": socket.gethostname(),
@@ -238,6 +239,27 @@ def _meta(nodes: int, ranks_per_node: int) -> dict[str, Any]:
     return meta
 
 
+def _record_row_store(propagator: Any) -> None:
+    """Fold one propagator's resolved row backend into this run's metadata.
+
+    ``monoprop_ROW_STORE`` says what was asked for, not what ran: unset lets the mode width pick, and
+    the crossover it picks against is a build-time constant. The two backends accumulate a term sum in
+    different orders and have different footprints, so a report has to name the one that ran. Widths
+    differ within a run, hence so can the backend: a disagreement records as ``"mixed"`` rather than
+    letting the last propagator speak for the others.
+    """
+    if _rank() != 0:
+        return
+    # Read straight off the binding, with no getattr fallback: a benchmark whose whole job is to name
+    # the backend that ran must fail loudly against an extension that cannot say, not quietly record
+    # nothing.
+    resolved = "sparse" if propagator._simulator.rows_are_sparse else "dense"
+    seen = _RESULTS["meta"].get("row_store_effective")
+    _RESULTS["meta"]["row_store_effective"] = (
+        resolved if seen in (None, resolved) else "mixed"
+    )
+
+
 def _params(config: pytest.Config) -> dict[str, Any]:
     """Return the resolved random-problem hyperparameters (defaults included)."""
     return {
@@ -353,6 +375,7 @@ def _record_model_stats(
 ) -> None:
     """Record term count, operator memory breakdown and footprint under ``key``."""
     _record("opsize", key, {"terms": _reduce_sum(comm, propagator.size())})
+    _record_row_store(propagator)
 
     # Placement is only observable while the propagator's threads are alive.
     _record_placement(comm)
@@ -553,6 +576,7 @@ def built_graph(
 
     # Under MPI the operator is partitioned, so sum the partitions.
     _record("opsize", picture, {"terms": _reduce_sum(bench_comm, mp.size())})
+    _record_row_store(mp)
 
     # Settled RSS once the build's transients are released -- the persistent
     # footprint the per-operation peak cannot see.

cpp/include/monoprop/MonomialPropagator.h

@@ -156,14 +156,30 @@ class MonomialPropagator {
     /// graph_layers(), optimizer order) or a per-gate one (length n_gates()); on a tie, per-layer wins.
     auto set_parameter_mapping(const VecZ &parameter_mapping) -> void;
 
-    /// This rank's monomial → coefficient index. Single-partition only — see require_single_partition_.
-    auto indexing() -> detail::OperatorIndex<NumModes> & {
-        require_single_partition_("indexing()");
-        return *mp_op_.store;
+    /// This rank's terms as fn(monomial, coefficient index), in the index's own slot order.
+    /// Single-partition only — see require_single_partition_. No accessor for the store itself: which
+    /// backend holds the rows is a runtime choice (see MPOperator::with_store), so there is no one type
+    /// to hand out.
+    template <typename Fn>
+    auto for_each_term(Fn &&fn) const -> void {
+        require_single_partition_("for_each_term()");
+        mp_op_.for_each_term(std::forward<Fn>(fn));
     }
-    auto indexing() const -> const detail::OperatorIndex<NumModes> & {
-        require_single_partition_("indexing()");
-        return *mp_op_.store;
+    /// This rank's term count. Single-partition only.
+    auto num_local_terms() const -> size_t {
+        require_single_partition_("num_local_terms()");
+        return mp_op_.size();
+    }
+
+    /// Whether this propagator's rows live in the support-form backend. Which one it is is decided once
+    /// at construction (see use_sparse_rows_); this reports the answer rather than re-deriving it.
+    /// Partition-transparent: every partition of a facade makes the same choice from the same width and
+    /// the same environment, so partition 0 speaks for all of them.
+    auto rows_are_sparse() const -> bool {
+        if (is_partition_facade()) {
+            return first_partition_().rows_are_sparse();
+        }
+        return mp_op_.rows_are_sparse();
     }
 
     /// Per-layer (cos_inds, local_cycles, cross_rank_sin_send, cross_rank_sin_recv) for this
@@ -272,7 +288,7 @@ class MonomialPropagator {
 
     /// Contract the graph into the operator (Heisenberg) or state (Schrodinger). `inplace` consumes the
     /// graph and updates internal state; otherwise nothing is mutated. Core term excluded either way.
-    /// Coefficients are positioned by the owning partition's indexing(), so on a facade the result is
+    /// Coefficients are positioned by the owning partition's row index, so on a facade the result is
     /// the per-partition blocks concatenated in partition order: the same multiset as an unpartitioned
     /// run, but not positionally stable across partition counts — and the count is auto-picked from the
     /// host's core count unless pinned. Use evolved_operator_terms() when positions must mean something.
@@ -312,8 +328,13 @@ class MonomialPropagator {
     detail::MatchedEpochSet matched_scratch_;
 
     // A perf hint, never a correctness constraint: overflow spills losslessly. Sized to the cutoff's
-    // structural position bound when it has one.
-    auto packed_inline_width_() const -> size_t;
+    // structural position bound when it has one. Shared by both backends -- the bound is in physical
+    // slots, which is what an OperatorIndex inline width and a SparseRowStore slot count both count.
+    auto row_width_bound_() const -> size_t;
+
+    // Which row backend to build on, decided once per propagator. See config::Settings::row_store for
+    // why an unrecognized monoprop_ROW_STORE is a throw rather than a silent fall back to auto.
+    auto use_sparse_rows_() const -> bool;
 
     // `requested` 0 ⇒ env/auto. Returns 1 for the ordinary single-partition path.
     static auto resolve_partition_count_(size_t requested, mpi::Comm comm) -> size_t;

cpp/monoprop/Bitset.h

@@ -216,21 +216,25 @@ class Bitset {
         return os;
     }
 };
+
+// The splitmix64 finalizer. Every hash in the engine ends here, and the value routes MPI ownership, so
+// this must stay bit-identical wherever it is reached from.
+[[nodiscard]] constexpr auto splitmix_finalize(uint64_t x) noexcept -> uint64_t {
+    x ^= x >> 30;
+    x *= 0xbf58476d1ce4e5b9ULL;
+    x ^= x >> 27;
+    x *= 0x94d049bb133111ebULL;
+    x ^= x >> 31;
+    return x;
+}
 } // namespace monoprop
 
 template <typename T>
 struct SplitmixHash;
 
 template <size_t NumBits>
 struct SplitmixHash<monoprop::Bitset<NumBits>> {
-    static constexpr auto mix(uint64_t x) noexcept -> uint64_t {
-        x ^= x >> 30;
-        x *= 0xbf58476d1ce4e5b9ULL;
-        x ^= x >> 27;
-        x *= 0x94d049bb133111ebULL;
-        x ^= x >> 31;
-        return x;
-    }
+    static constexpr auto mix(uint64_t x) noexcept -> uint64_t { return monoprop::splitmix_finalize(x); }
 
     auto operator()(const monoprop::Bitset<NumBits> &bs) const noexcept -> size_t {
         constexpr size_t W = monoprop::Bitset<NumBits>::num_words();

cpp/monoprop/CMakeLists.txt

@@ -22,6 +22,7 @@ target_compile_definitions(
   PUBLIC
     $<$<TARGET_EXISTS:MPI::MPI_CXX>:monoprop_ENABLE_MPI>
     $<$<BOOL:${monoprop_WIDE_TERM_INDEX}>:monoprop_WIDE_TERM_INDEX>
+    monoprop_SPARSE_ROW_MIN_MODES=${monoprop_SPARSE_ROW_MIN_MODES}
 )
 
 # flags to prepend
@@ -100,6 +101,7 @@ target_compile_definitions(
   INTERFACE
     $<$<TARGET_EXISTS:MPI::MPI_CXX>:monoprop_ENABLE_MPI>
     $<$<BOOL:${monoprop_WIDE_TERM_INDEX}>:monoprop_WIDE_TERM_INDEX>
+    monoprop_SPARSE_ROW_MIN_MODES=${monoprop_SPARSE_ROW_MIN_MODES}
 )
 
 target_compile_features(monoprop INTERFACE cxx_std_23)

cpp/monoprop/algebra/CMakeLists.txt

@@ -6,6 +6,7 @@ target_sources(
       FILES
         "Algebra.h"
         "AlgebraCommon.h"
+        "CodesAlgebra.h"
         "MajoranaAlgebra.h"
         "PauliAlgebra.h"
 )