From 672735154d1434a24a16a56f86f2dc2a3e76591e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Sat, 29 Aug 2026 16:43:13 +0000 Subject: [PATCH] =?UTF-8?q?perf(cpp):=20=E2=9A=A1=20bind=20the=20per-term?= =?UTF-8?q?=20kernel=20and=20the=20support-form=20query=20record?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Everything the scan asks about one anticommuting term -- the product, the overlap, the rotation sign, the structural cutoff, the owner rank and the query record -- goes through the per-gate object TermProductsFor selects, so the scan itself names no representation. SparseTermProducts answers the first four off the codes word and falls back to DenseTermProducts per term when there is no row to read or no codes form of the cutoff. The storage word count becomes the third thing bound once per layer, beside the algebra and the backend: with_kernel_width turns gen.num_words() into a template parameter W, so every per-term word loop has a compile-time trip count again. Worth ~10% at two and four storage words and nothing above, so kNarrowKernelWords caps which widths get an instantiation. A store is now queried in the form it keys its rows by, so a resolve never converts. A fully paired product escapes the cutoff and no fixed-stride sparse record can hold every one, so a buffer carries a dense escape tail behind its records; an escaped record keeps its place and its stride and carries its tail index where the codes word would go. Terms and energies are unchanged: the dense baseline stays byte-identical to main's and the support-form one still agrees to rtol 1e-10. Assisted-by: ClaudeCode:claude-opus-5 --- AGENTS.md | 73 +- cpp/monoprop/Bitset.h | 72 +- cpp/monoprop/algebra/Algebra.h | 29 + cpp/monoprop/algebra/AlgebraCommon.h | 28 + .../evolution/layer_build/CMakeLists.txt | 1 + .../detail/evolution/layer_build/Common.h | 283 +++++++- .../detail/evolution/layer_build/Engine.h | 93 ++- .../detail/evolution/layer_build/Resolve.h | 21 +- .../detail/evolution/layer_build/Scan.h | 93 ++- .../evolution/layer_build/TermProduct.h | 443 ++++++++++++ cpp/monoprop/detail/operator/SparseRowStore.h | 5 +- cpp/tests/InlineWidths.h | 43 ++ cpp/tests/README.md | 6 +- cpp/tests/bitset_tests.cpp | 2 +- cpp/tests/boostAddTests.cmake | 11 +- cpp/tests/evolution_detail_tests.cpp | 4 +- cpp/tests/sparse_wire_tests.cpp | 534 +++++++++++++++ cpp/tests/term_product_tests.cpp | 636 ++++++++++++++++++ cpp/tests/word_kernel_tests.cpp | 204 ++++++ 19 files changed, 2462 insertions(+), 119 deletions(-) create mode 100644 cpp/monoprop/detail/evolution/layer_build/TermProduct.h create mode 100644 cpp/tests/InlineWidths.h create mode 100644 cpp/tests/sparse_wire_tests.cpp create mode 100644 cpp/tests/term_product_tests.cpp create mode 100644 cpp/tests/word_kernel_tests.cpp diff --git a/AGENTS.md b/AGENTS.md index bafccbc8..c7c92020 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -184,13 +184,69 @@ Key files: `CutoffEvaluator::max_mode_bound()` **plus the generator's locality**; past its 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. -- **The query record**: a query goes on the wire as a dense monomial whichever backend holds the rows, so - a support-form row is materialized before it is pushed — `query_payload_words_for(store)` - (`layer_build/Common.h`) is the one width, and `DenseQueryKeys` the one batch. A buffer is - `[nq][record 0]…[record nq-1]`, and the record count comes off that header rather than `size/stride`, - so a record form that appends anything past the last record stays readable without changing a reader. +- **The per-term kernel seam**: everything the scan asks about one anticommuting term — the product, + the overlap, the rotation sign, the structural cutoff, the owner rank and the query record — goes + through the per-gate object `TermProductsFor` selects (`layer_build/TermProduct.h`), so the + scan itself names no representation. `SparseTermProducts` answers the first four off the `codes` word + and falls back to `DenseTermProducts` per term when there is no row to read (a spilled store row, a + product past the scratch capacity) or no codes form of the cutoff (`CutoffEvaluator` recovered neither + concrete functor, e.g. under a basis change). `cpp/tests/term_product_tests.cpp` compares the two + kernels answer for answer: extend it with any new answer, or that answer ships untested. +- **The third thing bound once per layer**: the storage word count, beside the algebra and the backend. + `with_kernel_width` (`layer_build/TermProduct.h`, over `detail::with_nwords`) turns + `gen.num_words()` into a template parameter `W` at the same seam in `build_layer`, and + `fused_find_and_collect` and the `TermProductsFor` specialization + `DenseTermProductsW` are templated on it, so every per-term word loop has a compile-time trip + count and every operand's storage pointer is resolved once per gate — which is what a `Bitset` + used to give for free. Measured worth ~10% at two and four storage words and nothing at seven or eight, + so `kNarrowKernelWords` (`TermProduct.h`, 4) caps which widths get an instantiation; the cost is ~11% + of `.text`. Not a build option, unlike `monoprop_SPARSE_ROW_MIN_MODES`: the cap is a *storage word* + count, so it names the same width regime on every machine, where the sparse crossover has to follow + the target ISA. Two conditions on that number, both measured. It is the **Majorana** path: + the Pauli rotation sign already loops over the generator's non-zero words only, so `W` binds no trip + count there and the 127-qubit kicked-Ising model gains ~1%. And it scales with how much of a run is in + the per-term product at all, so a loose `lower_atol` — which rejects a term on its coefficient before + the product is computed — sees about a third of it. Three consequences for the code. + `DenseTermProductsW` specializes the cutoff only for a **length** cutoff over the **whole register**; a + support cutoff or a narrower active window keeps going through `CutoffEvaluator`, and + `uses_word_cutoff()` is how a test tells those apart. (A support arm was tried and measured: ~1% worse + everywhere, and no gain even on the Pauli models that use it, because their per-term time is not in the + cutoff.) `WordKernel` (`Bitset.h`) is the four word ops with `W` fixed that stand in for a `Bitset` + *method* — standing in for one is the membership rule, which is why the scan's fifth bound-width pass, + `fully_paired_words`, lives in `algebra/AlgebraCommon.h` beside the `cutoff_sums` it answers for and + the even-bit literal it shares with `CutoffMasks::make`. Two of the four — the fused XOR and the + AND-fold behind `parity_and` — are the *same* definitions `Bitset`'s own inline arms use + (`detail::fused_xor_words` / `and_fold_words`, declared ahead of both), because one of them decides + emitted term signs and neither may drift from the method it stands in for. `splitmix` is deliberately + a second implementation instead: that value is `monomial_hash`, so it routes MPI ownership and must + stay bit-identical, and `cpp/tests/word_kernel_tests.cpp` asserts it equal to `SplitmixHash` at every + `W` rather than by construction. `term_product_tests.cpp` compares the whole kernel against + `DenseTermProducts` over the whole inline regime, not just the capped widths — both files sweep the + regime through the one `test_utils::for_each_inline_width` in `cpp/tests/InlineWidths.h`, so the range + cannot be narrowed in one of them alone. And the kernel's precondition is that every operand is + inline, so `W` is never bound above `Bitset::kInlineWords` — which `kNarrowKernelWords` is + `static_assert`ed against in the same header, so lowering `kInlineWords` fails the compile rather than + silently specializing a spilled width. `DenseTermProductsW` is non-copyable because its three word + pointers point into its own bitsets; a copy would read and write the original's storage. + Two further bindings were tried past this seam and both measured at nothing — under 0.05% of the + instruction count on either shipping model, pinned single-threaded — because the optimizer already + hoists them out of the inlined scan loop: resolving the algebra's per-term sign inputs into a + per-gate struct the kernel holds (`A::SignContext`), and writing the query record with the word count + bound (`query_push_words`, one capacity check instead of one per word). Measure any third one the + same way before adding it; wall clock cannot see this range, and neither can an instruction count taken + with the thread pool live, which spins hard enough to inflate the total ~14x. +- **The query record**: a store is queried in the form it keys its rows by, so a resolve never converts — + `QueryKeysFor` (`layer_build/Common.h`) picks the batch, and `query_payload_words_for(store, + capacity)` the width. A buffer is `[nq][record 0]…[record nq-1][dense escape tail]`: the header, + because a tail means `size/stride` is no longer the record count; the tail, because a query is `M ⊕ G` + and a fully paired product escapes the cutoff, so no fixed-stride sparse record can hold every one. An + escaped record keeps its place and its stride, marks lane 0 with `SparseRowStore::kOverflowLane` and + carries its tail *index* where the codes word would go — an index into the tail, never an offset into + the buffer, which is what lets the fused sink widen every record without renumbering anything. Push + records through `TermProducts::push(QueryOut{records, escapes}, phase)` and finish a stream with + `append_escape_tail`; never append a record after the tail has started. A batch's *retained* keys — the ones the deferred self-miss list reads after the slots have been - refilled — are a flat word arena at the batch's own width, never a container of + refilled — are a flat word arena at the batch's own width in both forms, never a container of monomials: a `Bitset` is sized for the widest inline width whatever its own is, and one key is retained per term a layer inserts, so a `MonomialList` there carried 72 bytes where a 128-bit monomial needs 16. It is worth -1.2% user instructions and cycles on the @@ -209,8 +265,9 @@ Key files: 8%, and across sessions -4% to +10% — so a real 1% shows up there as nothing. And a multi-threaded reading inverts: the partition pool spins, so `cycles:u` *rises* on a change that lowers wall time. Pin `monoprop_NUM_THREADS=1` and count instructions. - Owner routing is `monomial_hash` everywhere, including `find_rank`, so a multi-rank run materializes - one monomial per surviving term; moving that means changing `find_rank` too. + `owner()` is still the dense `monomial_hash` on both sides, because owner routing is that hash + everywhere including `find_rank` — so a multi-rank run still materializes one monomial per surviving + term, and moving that means changing `find_rank` too. - **The anticommutation fold** (`detail::InvertedIndex`): the transpose of the row store, one column per *bit position* — not per mode. That keying is settled and measured: a mode-keyed column cannot answer a generator slot that names one Majorana of a mode, which is 66% of the Hubbard generators' slots and diff --git a/cpp/monoprop/Bitset.h b/cpp/monoprop/Bitset.h index 11697b54..a85e06c5 100644 --- a/cpp/monoprop/Bitset.h +++ b/cpp/monoprop/Bitset.h @@ -60,9 +60,10 @@ struct FusedWordCounts { size_t result_count; // popcount(a ^ b) }; -// The word passes behind Bitset's own inline arms, defined once ahead of them and reached through -// with_nwords. One of them decides emitted term signs and the other feeds a cutoff, so a second -// definition that could drift from either is not acceptable. +// The word passes shared by Bitset's own inline arms and by the per-gate WordKernel below. They live +// here, ahead of both, so each is defined once: Bitset reaches them through with_nwords and the kernel +// through its bound W, and the two must not be able to drift -- one of them decides emitted term signs +// and the other feeds a cutoff. // // `nwords` is the exact word count of every operand, passed either as a std::integral_constant (the // inline regime, where the trip count is then a compile-time one) or as a plain size_t (a spilled @@ -118,7 +119,8 @@ namespace monoprop { class Bitset { public: // The word vocabulary is public because callers reason in words: data() already hands out a - // word_type*, and callers that walk words need the same three names to say what they walk. + // word_type*, and a per-gate kernel that binds the word count needs the same three names to say + // what it binds (see detail::WordKernel). using word_type = uint64_t; static constexpr auto word_width = sizeof(word_type) * 8; static constexpr size_t kInlineWords = 8; @@ -626,6 +628,68 @@ struct SplitmixHash { } }; +namespace detail { + +// The per-term word ops with the word count supplied by the caller instead of read off the operand. +// +// The arithmetic is identical to Bitset's own methods; what differs is what the compiler knows. A +// Bitset method must load nwords_, compare it against the inline capacity and select a storage +// pointer on every call, and none of those three can be hoisted out of a loop the optimizer cannot +// see through -- which, on the per-term path, is every call. Handing a kernel a compile-time W and +// the word pointers the caller resolved once leaves a straight-line unrolled loop, which is what the +// per-width Bitset got for free. +// +// Preconditions, none of them checkable here: every pointer is a Bitset::data() of a bitset of +// exactly W words, and W <= kInlineWords so no operand is spilled. The only legal caller is one that +// bound W from a width it owns for the whole loop -- see DenseTermProductsW, which is the seam that +// binds it once per gate. +// +// Standing in for a Bitset method is also what decides membership: the scan's other bound-width word +// pass, fully_paired_words, answers a question about the algebra rather than the storage, so it lives +// beside its own oracle in AlgebraCommon.h instead. +template +struct WordKernel { + static_assert(W >= 1 && W <= Bitset::kInlineWords, + "the kernel covers the inline regime; above it the runtime loop already wins"); + + using word_type = Bitset::word_type; + + static auto clear(word_type *a) noexcept -> void { + for (size_t i = 0; i < W; ++i) { + a[i] = 0; + } + } + + // Bitset::fused_xor_into with W fixed -- the same pass, reached without the nwords_ load and + // storage-pointer select the method does per call. + static auto fused_xor_into(const word_type *a, const word_type *b, word_type *out) noexcept -> Bitset::FusedCounts { + return fused_xor_words(a, b, out, std::integral_constant{}); + } + + // Bitset::parity_and with W fixed, likewise. + [[nodiscard]] static auto parity_and(const word_type *a, const word_type *b) noexcept -> bool { + return (std::popcount(and_fold_words(a, b, std::integral_constant{})) & 1U) != 0; + } + + // SplitmixHash with W fixed. Must stay bit-identical to it: this value routes MPI ownership, so a + // divergence would move terms between ranks rather than merely run slower. Hence the W == 1 arm + // reproducing the same special case rather than folding into the loop. + [[nodiscard]] static auto splitmix(const word_type *a) noexcept -> size_t { + if constexpr (W == 1) { + return static_cast(SplitmixHash::mix(a[0])); + } + else { + uint64_t h = 0; + for (size_t i = 0; i < W; ++i) { + h ^= SplitmixHash::mix(a[i] + static_cast(i)); + } + return static_cast(h); + } + } +}; + +} // namespace detail + } // namespace monoprop namespace std { diff --git a/cpp/monoprop/algebra/Algebra.h b/cpp/monoprop/algebra/Algebra.h index aa42dff8..8948ce0a 100644 --- a/cpp/monoprop/algebra/Algebra.h +++ b/cpp/monoprop/algebra/Algebra.h @@ -24,6 +24,7 @@ #include #include "monoprop/algebra/AlgebraCommon.h" +#include "monoprop/algebra/CodesAlgebra.h" #include "monoprop/algebra/MajoranaAlgebra.h" #include "monoprop/algebra/PauliAlgebra.h" #include "monoprop/core/Monomial.h" @@ -54,6 +55,20 @@ struct MajoranaAlgebra { static auto rotation_sign(const GenContext &ctx, const Bitset &mono, const Bitset & /*new_mono*/) -> int { return mono.parity_and(ctx.interleave_mask) ? -1 : 1; } + // The same sign off word pointers the caller resolved once, with the word count bound by the + // caller rather than read off the operand. Same fold, same parity; see detail::WordKernel. + template + static auto rotation_sign_words(const GenContext &ctx, + const Bitset::word_type *mono, + const Bitset::word_type * /*new_mono*/) -> int { + return detail::WordKernel::parity_and(mono, ctx.interleave_mask.data()) ? -1 : 1; + } + // The same sign in support form. No GenContext: the interleave mask is dense by construction + // (roughly half the register), so the sparse form walks the two rows instead of carrying a mask, and + // the product row is not an argument either -- see codes_interleave_phase. + static auto codes_rotation_sign(const detail::SparseRow &mono, const detail::SparseRow &gen) -> int { + return detail::codes_interleave_phase(mono, gen); + } static auto emit_phase(int rotation_sign, size_t mono_pop, size_t gen_pop, size_t overlap) -> int { return rotation_sign * hermitian_phase(mono_pop, gen_pop, overlap); } @@ -78,6 +93,20 @@ struct PauliAlgebra { static auto rotation_sign(const GenContext &ctx, const Bitset &mono, const Bitset &new_mono) -> int { return pauli_rotation_sign(ctx.pauli_ctx, mono, new_mono); } + // W is unused here and that is the point: this sign already loops over the generator's non-zero + // words only, so there is no trip count to bind -- what the word form removes is the storage-pointer + // select that mono.word(w) repeats on every access. + template + static auto rotation_sign_words(const GenContext &ctx, + const Bitset::word_type *mono, + const Bitset::word_type *new_mono) -> int { + return pauli_rotation_sign_words(ctx.pauli_ctx, mono, new_mono); + } + // Same exponent as above off the two rows; new_mono never has to exist, since a mode the generator + // misses contributes nothing (see codes_pauli_rotation_sign). + static auto codes_rotation_sign(const detail::SparseRow &mono, const detail::SparseRow &gen) -> int { + return detail::codes_pauli_rotation_sign(mono, gen); + } // Pauli's rotation sign is already the emitted sine phase -- no Hermitian fold. static auto emit_phase(int rotation_sign, size_t /*mono_pop*/, size_t /*gen_pop*/, size_t /*overlap*/) -> int { return rotation_sign; diff --git a/cpp/monoprop/algebra/AlgebraCommon.h b/cpp/monoprop/algebra/AlgebraCommon.h index 8a399493..44d4d7f4 100644 --- a/cpp/monoprop/algebra/AlgebraCommon.h +++ b/cpp/monoprop/algebra/AlgebraCommon.h @@ -170,6 +170,34 @@ auto support_cutoff(const MonomialLike auto &mono, unsigned int cutoff) -> bool namespace detail { +// The xor_sum == 0 clause above -- "every occupied mode has both its Majoranas" -- with the storage +// word count bound at compile time, for the per-gate scan kernel that has W and the words already (see +// DenseTermProductsW). It lives here and not beside the other bound-width word ops in Bitset.h because +// "paired" is a fact about the algebra and not about the storage, and because the mask literal it shares +// with CutoffMasks::make above is easier to keep honest in one file than in two. +// +// Folded with OR and tested against zero rather than summing popcounts: the caller only ever compares +// the sum to zero, and the two agree because each per-word term is non-negative. +// +// The even-bit mask is the literal rather than an argument, which is what makes this W loads instead of +// 2W: a storage width is a whole number of words, so even_bits is this pattern in every one of +// them. Its top-word trim at a non-word-multiple width is unobservable here -- bits above the logical +// width are never set, so they pair with themselves either way. (word >> 1) & mask cannot cross a word +// because a mode's two bits are 2m and 2m+1. +// +// Whole register only: a narrower active window would need the shift cutoff_sums applies, and getting +// it wrong would silently change which terms survive. +template +[[nodiscard]] [[gnu::always_inline]] inline auto fully_paired_words(const Bitset::word_type *a) noexcept -> bool { + constexpr Bitset::word_type kEven = 0x5555555555555555ULL; + Bitset::word_type unpaired = 0; + for (size_t i = 0; i < W; ++i) { + const Bitset::word_type word = a[i]; + unpaired |= (word & kEven) ^ ((word >> 1) & kEven); + } + return unpaired == 0; +} + // Both hold their masks, so the per-term call does no width arithmetic. Real constructors rather than // aggregate initialization, deliberately: the width used to arrive free from NumModes, and both // remaining ways to get it wrong are silent. A logical width of 0 makes cutoff_sums' active window diff --git a/cpp/monoprop/detail/evolution/layer_build/CMakeLists.txt b/cpp/monoprop/detail/evolution/layer_build/CMakeLists.txt index eb4c9853..33842166 100644 --- a/cpp/monoprop/detail/evolution/layer_build/CMakeLists.txt +++ b/cpp/monoprop/detail/evolution/layer_build/CMakeLists.txt @@ -9,4 +9,5 @@ target_sources( "FusedApply.h" "Resolve.h" "Scan.h" + "TermProduct.h" ) diff --git a/cpp/monoprop/detail/evolution/layer_build/Common.h b/cpp/monoprop/detail/evolution/layer_build/Common.h index 601618b6..38f4c23f 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Common.h +++ b/cpp/monoprop/detail/evolution/layer_build/Common.h @@ -138,6 +138,12 @@ inline constexpr size_t kQueryHeaderWords = 1; [[nodiscard]] constexpr auto query_record_offset(size_t q, size_t stride) -> size_t { return kQueryHeaderWords + (q * stride); } +// Where the escape tail starts: right after the last record. Both sides derive it from the header and the +// stride, so an escape's own index into the tail is independent of either -- which is what lets the fused +// re-layout below move the records without touching the tail. +[[nodiscard]] inline auto query_tail_offset(const VecZ &buf, size_t stride) -> size_t { + return query_record_offset(query_record_count(buf), stride); +} // Fused query+value record width (R>1): the plain query record plus one trailing word holding the source's // pre-cos coeff (v_src, bit-cast from double), so query + value ride a single alltoallv instead of two. @@ -163,8 +169,25 @@ inline auto decode_value(size_t word) -> double { return std::bit_cast(word); } -// Appends one record and bumps the header, so the count is always the buffer's own rather than -// something a reader has to derive from size and stride. +// The two buffers a scan pushes a query into. A record cannot be appended once the tail has started, so +// the tail accumulates separately and is concatenated when the scan is done -- an escape's index names a +// position within the tail, so the concatenation moves nothing it refers to. +struct QueryOut { + VecZ &records; + VecZ &escapes; +}; + +// The scan's single exit: fold each stream's escapes in behind its records. +inline auto append_escape_tail(VecZ &records, VecZ &escapes) -> void { + if (escapes.empty()) { + return; + } + records.insert(records.end(), escapes.begin(), escapes.end()); + escapes.clear(); +} + +// Appends one record and bumps the header. Every push must precede the escape tail, which the scan +// guarantees by collecting escapes in a buffer of their own and concatenating once the scan is done. inline auto query_push(VecZ &buf, const Bitset &mono, int phase) -> void { assert(buf.size() >= kQueryHeaderWords && "a query buffer must be created with query_buffer()"); mpi_detail::append_monomial_words(mono, buf); @@ -182,6 +205,94 @@ inline auto query_read(const VecZ &buf, size_t q, size_t stride, Bitset &mono_ou phase_out = decode_phase(buf[base + mono_out.num_words()]); } +// A query record in support form. Deliberately the same *shape* as the dense one -- a fixed-stride +// payload followed by the phase word -- so every stride computation, alltoallv count and reader offset +// above holds with `num_words` reinterpreted as the payload word count. Only the payload differs: +// `sparse_lane_words` words of four uint16 mode lanes each plus one codes word, against the monomial's +// full word count. That is what the format is for: a 1024-mode monomial is 32 words, where a 12-slot row +// is 3 lane words plus 1. +// +// The stride needs a capacity, and it must be one every rank derives identically without communication -- +// SparseRowStore::scratch_slots_for(cutoff mode bound, widest generator), the same value the scan's +// scratch row uses, since a query carries exactly such a product. Ranks already owe each other this kind +// of agreement for the hash width (see find_rank). +// +// A row past that capacity has no sparse record: the caller must fall back to the dense one, and the +// push below asserts rather than truncating. +constexpr auto sparse_lane_words(size_t capacity) -> size_t { + return (capacity + 3) / 4; +} +constexpr auto sparse_payload_words(size_t capacity) -> size_t { + return sparse_lane_words(capacity) + 1; +} + +inline auto sparse_query_push(VecZ &buf, const SparseRow &row, size_t capacity, int phase) -> void { + assert(buf.size() >= kQueryHeaderWords && "a query buffer must be created with query_buffer()"); + const size_t lane_words = sparse_lane_words(capacity); + const size_t n = row.num_slots(); + assert(n <= capacity && "sparse_query_push row exceeds the record capacity"); + const size_t base = buf.size(); + // Zero-filled, so the lanes past the row's own are deterministic; the reader ignores them, taking the + // slot count off the codes word. + buf.resize(base + lane_words + 2, 0); + for (size_t j = 0; j < n; ++j) { + buf[base + (j / 4)] |= static_cast(row.modes[j]) << (16 * (j % 4)); + } + buf[base + lane_words] = static_cast(row.codes); + buf[base + lane_words + 1] = encode_phase(phase); + ++buf[0]; +} + +// The record left behind by a query no sparse record can hold, and it is not a corner case to be sized +// away: a query is M ⊕ G and a fully paired product escapes the cutoff, so nothing bounds its support. +// +// The record keeps its place and its stride -- which is what leaves every offset, alltoallv count and +// compaction in the engine as plain arithmetic -- and says where to find the monomial instead: lane 0 +// carries the store's own overflow marker and the codes slot carries the escape's index into the buffer's +// tail. Both conventions are SparseRowStore's for a spilled row, deliberately: one representation of "too +// wide for a codes word", not two. +// +// `escapes` accumulates the tail separately during the scan, because a record cannot be appended after +// the tail has started; the scan concatenates the two once it is done. An escape's index is its position +// in that tail, so it survives the concatenation and the fused re-layout alike. +inline auto sparse_query_push_escape(VecZ &buf, VecZ &escapes, const Bitset &mono, size_t capacity, int phase) -> void { + assert(buf.size() >= kQueryHeaderWords && "a query buffer must be created with query_buffer()"); + const size_t lane_words = sparse_lane_words(capacity); + const size_t base = buf.size(); + buf.resize(base + lane_words + 2, 0); + buf[base] = static_cast(SparseRowStore::kOverflowLane); + buf[base + lane_words] = escapes.size() / mono.num_words(); + buf[base + lane_words + 1] = encode_phase(phase); + ++buf[0]; + mpi_detail::append_monomial_words(mono, escapes); +} + +// Reading a record's shape needs the lane word count, which the reader has from the capacity; taking the +// record base as an argument keeps this usable from both the plain and the fused stride. +[[nodiscard]] inline auto sparse_record_is_escape(const VecZ &buf, size_t base) -> bool { + return buf[base] == static_cast(SparseRowStore::kOverflowLane); +} + +// lanes_out must hold `capacity` lanes. Writes only the row's own, for the same reason +// read_monomial_from_words overwrites whole words: what a previous record left beyond them cannot be read. +inline auto sparse_query_read(const VecZ &buf, + size_t q, + size_t stride, + size_t capacity, + RowMode *lanes_out, + RowCodes &codes_out, + int &phase_out) -> void { + const size_t base = query_record_offset(q, stride); + const size_t lane_words = sparse_lane_words(capacity); + assert(!sparse_record_is_escape(buf, base) && "an escaped record has no row to read"); + codes_out = static_cast(buf[base + lane_words]); + const size_t n = row_slot_count(codes_out); + for (size_t j = 0; j < n; ++j) { + lanes_out[j] = static_cast((buf[base + (j / 4)] >> (16 * (j % 4))) & 0xFFFFU); + } + phase_out = decode_phase(buf[base + lane_words + 1]); +} + // Owned storage for a batch of query keys in whichever form a store keys its rows by, plus the record // reader that fills it. Both resolve paths -- the self-resolve batch in Engine.h and the incoming probe in // Resolve.h -- want exactly this, and both used to hand-roll it: allocate once, keep the storage across @@ -204,8 +315,9 @@ class DenseQueryKeys { using key_type = Bitset; // num_bits is the monomial storage width: query_read memcpys the destination's full word count, so the - // destination is what fixes the record width. - auto configure(size_t num_bits) -> void { + // destination is what fixes the record width. capacity is the support form's row capacity, which a + // dense record has no use for -- both batches take both so the call sites need no branch. + auto configure(size_t num_bits, size_t /*capacity*/) -> void { if (extent_ != num_bits) { keys_.clear(); retained_words_.clear(); @@ -265,12 +377,159 @@ class DenseQueryKeys { size_t extent_ = 0; }; -// The payload width of one query record, in VecZ words -- the quantity every stride, alltoallv count and -// record offset in Engine.h derives from. Both backends are queried with dense monomials, so the width is -// the store's own: a support-form row still has to be materialized before it goes on the wire. -[[nodiscard]] inline auto query_payload_words_for(const auto &store) -> size_t { +// The support-form counterpart: one lane array for the whole batch plus a parallel array of keys viewing +// into it, since find_batch wants the keys contiguous and a SparseRow is only a pointer and a word. +// +// The escaped records are why the key is a SparseRowKey rather than a SparseRow: one of those queries has +// no row, so its key points at a monomial materialized out of the buffer's tail. That storage is a deque +// on purpose -- push_back must not invalidate a key handed out for an earlier slot, and a vector's would. +class SparseQueryKeys { +public: + using key_type = SparseRowKey; + + // capacity is the row capacity in slots -- the same value that fixes the record stride, since a record + // holds exactly that many lanes. num_bits sizes the escape monomials. + auto configure(size_t num_bits, size_t capacity) -> void { + if (capacity_ != capacity || num_bits_ != num_bits) { + lanes_.clear(); + keys_.clear(); + escapes_.clear(); + retained_lanes_.clear(); + retained_bases_.clear(); + retained_escapes_.clear(); + retained_.clear(); + capacity_ = capacity; + num_bits_ = num_bits; + } + } + auto ensure(size_t n) -> void { + if (keys_.size() >= n) { + return; + } + lanes_.resize(n * capacity_); + keys_.resize(n); + // Every view is rebuilt, not just the new tail: the resize above may have moved lanes_, which + // would leave the existing views pointing into freed storage. + for (size_t i = 0; i < n; ++i) { + keys_[i] = SparseRowKey{.row = SparseRow{&lanes_[i * capacity_], 0}}; + } + } + // Hygiene rather than correctness: a key always points at the entry read for it, so stale entries are + // unreachable -- they would just accumulate for the whole resolve. Dropped per batch because that is + // the granularity at which slots are refilled anyway. + auto begin_batch() -> void { escapes_.clear(); } + [[nodiscard]] auto read_record(const VecZ &buf, size_t q, size_t stride, size_t slot) -> int { + const size_t base = query_record_offset(q, stride); + const size_t lane_words = sparse_lane_words(capacity_); + if (sparse_record_is_escape(buf, base)) { + // The tail entry this record named. Its offset needs the record count and the stride, both of + // which the buffer and the caller already carry. + const size_t tail = query_tail_offset(buf, stride); + const size_t words = Bitset::words_for(num_bits_); + escapes_.emplace_back(num_bits_); + mpi_detail::read_monomial_from_words(buf, tail + (buf[base + lane_words] * words), escapes_.back()); + keys_[slot].spilled = &escapes_.back(); + return decode_phase(buf[base + lane_words + 1]); + } + int phase = 0; + keys_[slot].spilled = nullptr; + // The lane pointer stays the one ensure() set: a record's lanes are read into this slot's own run. + sparse_query_read(buf, q, stride, capacity_, &lanes_[slot * capacity_], keys_[slot].row.codes, phase); + return phase; + } + [[nodiscard]] auto data() const -> const key_type * { return keys_.data(); } + [[nodiscard]] auto operator[](size_t slot) const -> const key_type & { return keys_[slot]; } + + // See DenseQueryKeys::retain. A retained key owns its lanes here too, in a second arena, and an escaped + // one owns its monomial -- the batch's own escape storage is dropped every batch. + [[nodiscard]] auto retain(size_t slot) -> size_t { + const size_t handle = retained_.size(); + // A base is recorded for every handle, escaped or not, so retained_bases_ stays indexable by + // handle; an escaped key's is simply never read. + retained_bases_.push_back(retained_lanes_.size()); + if (keys_[slot].is_spilled()) { + retained_escapes_.push_back(*keys_[slot].spilled); + // .row left empty: is_spilled() sends every read to .spilled instead. + retained_.push_back(SparseRowKey{.row = {}, .spilled = &retained_escapes_.back()}); + return handle; + } + retained_lanes_.resize(retained_bases_.back() + capacity_); + const size_t n = keys_[slot].row.num_slots(); + std::copy_n(keys_[slot].row.modes, + n, + retained_lanes_.begin() + static_cast(retained_bases_.back())); + // The lane array grows, so a key cannot hold a pointer into it; retained() rebuilds the view. + retained_.push_back(SparseRowKey{.row = SparseRow{nullptr, keys_[slot].row.codes}}); + return handle; + } + [[nodiscard]] auto retained(size_t handle) const -> key_type { + const key_type &key = retained_[handle]; + if (key.is_spilled()) { + return key; + } + return SparseRowKey{.row = SparseRow{&retained_lanes_[retained_bases_[handle]], key.row.codes}}; + } + +private: + DefaultInitVector lanes_ = {}; + std::vector keys_ = {}; + std::deque escapes_ = {}; + // Retained keys, whose storage must outlive the batch's own (see retain). retained_bases_ is parallel + // to retained_ but indexed only for the non-escaped ones -- an escaped key's entry is unread. + DefaultInitVector retained_lanes_ = {}; + std::vector retained_bases_ = {}; + std::deque retained_escapes_ = {}; + std::vector retained_ = {}; + size_t capacity_ = 0; + size_t num_bits_ = 0; +}; + +// A query key as a dense monomial, for the handful of places that need one -- the Schrodinger fresh-insert +// scoring, which has no codes form. Returns a reference when the key already is one and a value otherwise, +// so callers bind with `const auto &` to extend the temporary, exactly as materialize_row documents. +[[nodiscard]] inline auto key_monomial(const Bitset &key, size_t /*num_bits*/) -> const Bitset & { + return key; +} +[[nodiscard]] inline auto key_monomial(const SparseRowKey &key, size_t num_bits) -> Bitset { + if (key.is_spilled()) { + return *key.spilled; + } + return sparse_row_to_bitset(key.row, num_bits); +} + +// The payload width of one query record for a store, in VecZ words -- the quantity every stride, alltoallv +// count and record offset in Engine.h derives from. An overload per store rather than one accessor on +// MPOperator, because it is a property of the wire format the store is queried through, not of the store. +// +// Dense rows put the monomial's own words on the wire. The support form will put lane words plus the codes +// word: 5 words for a 12-slot row against 33 for a 1024-mode monomial, but slightly *wider* just above the +// store's own crossover -- 5 against 4 at 96 modes, break-even near 128 modes. That band is why the record +// form is tied to the store rather than chosen per layer: a runtime record form would double the engine's +// template instantiations again, to save a word in a narrow range. +[[nodiscard]] inline auto query_payload_words_for(const OperatorIndex &store, size_t /*capacity*/) -> size_t { return Bitset::words_for(store.num_bits()); } +[[nodiscard]] inline auto query_payload_words_for(const SparseRowStore & /*store*/, size_t capacity) -> size_t { + return sparse_payload_words(capacity); +} + +// Which key batch a store's query records arrive in. Explicit specializations rather than a member +// typedef on the stores: the record codec lives here, and a store must not depend on the wire format it +// is queried through. +// +// Each store is queried in the form it keys its rows by, so a resolve never converts: the dense store +// receives monomials, the support form receives rows (and, for the queries no row can hold, the escape +// monomials its tail carries). +template +struct QueryKeysFor; +template <> +struct QueryKeysFor { + using type = DenseQueryKeys; +}; +template <> +struct QueryKeysFor { + using type = SparseQueryKeys; +}; // No monomial reconstruction: process_responses needs only the phase. inline auto query_phase(const VecZ &buf, size_t q, size_t num_words) -> int { @@ -282,6 +541,10 @@ inline auto query_value(const VecZ &buf, size_t q, size_t num_words) -> double { } // Requires v.size() == query_record_count(q): exactly one value per query record. +// +// The escape tail rides along unchanged. It can, because an escape's index names its position *within the +// tail* rather than an offset into the buffer -- so widening every record by a value word moves the tail +// without renumbering anything in it. inline auto build_fused_query_value(const VecZ &q, const std::vector &v, VecZ &out, size_t num_words) -> void { out.clear(); if (q.size() < kQueryHeaderWords) { @@ -291,7 +554,8 @@ inline auto build_fused_query_value(const VecZ &q, const std::vector &v, } const size_t W = query_words(num_words); const size_t nq = query_record_count(q); - out.reserve(kQueryHeaderWords + (nq * query_words_fused(num_words))); + const size_t tail = query_tail_offset(q, W); + out.reserve(kQueryHeaderWords + (nq * query_words_fused(num_words)) + (q.size() - tail)); out.push_back(nq); for (size_t i = 0; i < nq; ++i) { out.insert(out.end(), @@ -299,6 +563,7 @@ inline auto build_fused_query_value(const VecZ &q, const std::vector &v, q.begin() + static_cast(query_record_offset(i + 1, W))); out.push_back(encode_value(v[i])); } + out.insert(out.end(), q.begin() + static_cast(tail), q.end()); } } // namespace monoprop::detail diff --git a/cpp/monoprop/detail/evolution/layer_build/Engine.h b/cpp/monoprop/detail/evolution/layer_build/Engine.h index 051e97b0..36b38b76 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Engine.h +++ b/cpp/monoprop/detail/evolution/layer_build/Engine.h @@ -72,7 +72,11 @@ struct GraphSink { // The record stride is not a constant: it is derived from the // monomial word count now, which the sink is handed at construction. size_t num_words = 0; + // The support form's row capacity, which fixes what a record's lanes hold. Zero and unread for a dense + // record; carried by both sinks so the resolve path needs no branch to configure its key batch. + size_t record_capacity = 0; [[nodiscard]] auto stride() const -> size_t { return query_words(num_words); } + [[nodiscard]] auto capacity() const -> size_t { return record_capacity; } static auto init_response() -> Response { return std::numeric_limits::max(); } size_t R; @@ -82,8 +86,9 @@ struct GraphSink { size_t def_out_base_ = 0; std::vector in_base_; // cross-rank per-rank base into acc[s].in_entries (set in prepare) - GraphSink(size_t num_words_, size_t R_, size_t my_rank_) + GraphSink(size_t num_words_, size_t record_capacity_, size_t R_, size_t my_rank_) : num_words(num_words_), + record_capacity(record_capacity_), R(R_), my_rank(my_rank_), acc(R_) {} @@ -187,8 +192,10 @@ struct GraphSink { struct ContractSink { static constexpr bool wants_values = true; using Response = double; - size_t num_words = 0; // see GraphSink::num_words + size_t num_words = 0; // see GraphSink::num_words + size_t record_capacity = 0; // see GraphSink::record_capacity [[nodiscard]] auto stride() const -> size_t { return query_words_fused(num_words); } + [[nodiscard]] auto capacity() const -> size_t { return record_capacity; } static auto init_response() -> Response { return 0.0; } size_t R; @@ -250,7 +257,7 @@ struct ContractSink { else if (schrodinger) { // The one arm that needs a dense monomial: this scoring has no codes form, so a support-form // key materializes here. It is a fresh cross-rank Schrodinger miss, so per gate it is rare. - const auto &mono = pr.mono[g]; + const auto &mono = key_monomial(pr.mono[g], num_bits_); v_tgt = is_paired(mono) ? algebra_state_phase(basis, mono, state_mask_) : 0.0; } else { @@ -330,9 +337,11 @@ struct LayerBuildEngine { // own field rather than read off sink: Sink is a minimal concept (see RecordingSink in // evolution_detail_tests.cpp), and only the two production sinks happen to also carry num_words. size_t words_ = 0; + // The support form's row capacity (see GraphSink::record_capacity), for the key batch below. + size_t record_capacity_ = 0; // Self-resolve key batch, configured once off the operator. Separate from the incoming probe's batch: // the two are live at the same time on the resolve path. - DenseQueryKeys keys_; + typename QueryKeysFor::type keys_; LayerBuildEngine(MPOperator &local_op_, Store &store_, @@ -342,6 +351,7 @@ struct LayerBuildEngine { MatchedEpochSet &matched_scratch, size_t combined_size_, size_t payload_words, + size_t record_capacity, Sink &&sink_) : local_op(local_op_), store(store_), @@ -353,8 +363,9 @@ struct LayerBuildEngine { queries_r(R_), src_idx_r(R_), sink(std::move(sink_)), - words_(payload_words) { - keys_.configure(local_op_.num_bits()); + words_(payload_words), + record_capacity_(record_capacity) { + keys_.configure(local_op_.num_bits(), record_capacity_); matched.begin_gate(combined_size); } @@ -438,7 +449,17 @@ struct LayerBuildEngine { } ++kept; } - q.resize(query_record_offset(kept, W)); + // The escape tail follows the records down, contents untouched. Its entries keep their + // positions relative to each other, so the indices the surviving records carry stay right -- + // including the entries a dropped record orphaned, which ride along as unread padding rather + // than being renumbered. + const size_t tail = query_record_offset(nq, W); + const size_t tail_words = q.size() - tail; + const size_t kept_end = query_record_offset(kept, W); + std::copy(q.begin() + static_cast(tail), + q.end(), + q.begin() + static_cast(kept_end)); + q.resize(kept_end + tail_words); q[0] = kept; s.resize(kept); if (v != nullptr) { @@ -615,10 +636,11 @@ auto build_layer(auto &local_op, &out_cos, &fused_contract, &schrodinger](S &store) -> std::shared_ptr { - // The record width, derived once per layer and passed to everything that reads a record: the - // engine's key batch and the sink's strides. It is a function of the store alone, so every rank - // computes the same number without communication. - const size_t record_payload_words = query_payload_words_for(store); + // The record shape, derived once per layer and passed to everything that reads a record: the scan's + // kernel, the engine's key batch and the sink's strides. Both numbers are functions of the store, the + // generator and the cutoff, so every rank computes the same pair without communication. + const size_t record_capacity = sparse_record_capacity(gen, cut_eval); + const size_t record_payload_words = query_payload_words_for(store, record_capacity); FusedScanResult fused = [&] { double *const sweep_ptr = fused_scale ? fused_scale_coeffs->data() : nullptr; @@ -636,19 +658,38 @@ auto build_layer(auto &local_op, &use_fused, &sweep_ptr, &cos_build]() { - return fused_find_and_collect(local_op, - store, - gen, - cut_eval, - cut_st, - coeffs, - only_rotate_len_k, - R, - my_rank, - logical_num_modes, - /*capture_values=*/use_fused, - sweep_ptr, - cos_build); + // Third and last thing bound once per layer, beside the algebra and the + // backend: the storage word count, so the scan's per-term word loops have a + // compile-time trip count. + return with_kernel_width( + gen.num_words(), + [&local_op, + &store, + &gen, + &cut_eval, + &cut_st, + &coeffs, + &only_rotate_len_k, + &R, + &my_rank, + &logical_num_modes, + &use_fused, + &sweep_ptr, + &cos_build](std::integral_constant) { + return fused_find_and_collect(local_op, + store, + gen, + cut_eval, + cut_st, + coeffs, + only_rotate_len_k, + R, + my_rank, + logical_num_modes, + /*capture_values=*/use_fused, + sweep_ptr, + cos_build); + }); }); }(); @@ -675,6 +716,7 @@ auto build_layer(auto &local_op, matched_scratch, /*combined_size=*/store.size(), record_payload_words, + record_capacity, std::move(sink)); eng.run_exchange(/*is_leader_pass=*/true, std::move(fused.leader_queries), @@ -691,6 +733,7 @@ auto build_layer(auto &local_op, if (use_fused) { const double inv_cos = fused_scale ? 1.0 / cos_build : 1.0; // pre-cos recovery factor for hit v_tgt return run(ContractSink{.num_words = record_payload_words, + .record_capacity = record_capacity, .R = R, .my_rank = my_rank, .fc = *fused_contract, @@ -700,7 +743,7 @@ auto build_layer(auto &local_op, .schrodinger = schrodinger, .basis = basis}); } - return run(GraphSink{record_payload_words, R, my_rank}); + return run(GraphSink{record_payload_words, record_capacity, R, my_rank}); }); // Recompute metadata rides with the layer so it survives every graph transform. scaled_count is the diff --git a/cpp/monoprop/detail/evolution/layer_build/Resolve.h b/cpp/monoprop/detail/evolution/layer_build/Resolve.h index 915a349b..c0dff5dd 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Resolve.h +++ b/cpp/monoprop/detail/evolution/layer_build/Resolve.h @@ -32,9 +32,9 @@ namespace monoprop::detail { // the next index base+j in (sender,record) order, so the assignment (and multi-rank bit-exactness) cannot // drift between resolvers. Queries are source⊕G over globally-distinct sources, ⊕G injective ⇒ queries // pairwise distinct ⇒ misses distinct and absent. -// Key is the form a query record is read into -- a monomial for both stores, since a query goes on the wire -// densely. Named `Key` rather than fixed so a record form of the store's own can be added without -// re-typing every consumer. +// Key is whichever form the store being probed is keyed by (QueryKeysFor): a monomial for the dense store, +// a row-or-escape key for the support form. Named `Key` rather than fixed because the whole point of the +// two record forms is that a resolve never converts one into the other. template struct IncomingProbeT { std::vector goff; // rank_count+1 flat offsets: g = goff[s] + q @@ -52,7 +52,7 @@ struct IncomingProbeT { // The probe over a store, spelled once so callers need not name the key form. template -using IncomingProbeFor = IncomingProbeT; +using IncomingProbeFor = IncomingProbeT::type::key_type>; // Phases 1-2, read-only w.r.t. operator contents. query_stride is the per-record width: the plain query // width, or the fused one for the fused resolver. The caller runs Phase 3, then insert_incoming_misses. @@ -66,14 +66,15 @@ inline auto probe_incoming_queries(const std::vector &incoming, // seriali const MPOperator &op, Store &store, size_t rank_count, - size_t query_stride) -> IncomingProbeFor { - using Keys = DenseQueryKeys; + size_t query_stride, + size_t record_capacity) -> IncomingProbeFor { + using Keys = typename QueryKeysFor::type; IncomingProbeFor pr; pr.goff.assign(rank_count + 1, 0); for (size_t s = 0; s < rank_count; ++s) { - // Off the buffer's own header rather than size/stride, so a record form that appends anything past - // the last record stays readable here unchanged. + // Off the buffer's own header, not its size: a support-form buffer carries an escape tail after its + // records, so size/stride is not the record count. pr.goff[s + 1] = pr.goff[s] + query_record_count(incoming[s]); } pr.nq_total = pr.goff[rank_count]; @@ -111,7 +112,7 @@ inline auto probe_incoming_queries(const std::vector &incoming, // seriali // own `keys_` does not make this trade: it is a plain member of a LayerBuildEngine built fresh per // build_layer call, so it starts default-constructed every layer (see DenseQueryKeys's comment). thread_local Keys scratch; - scratch.configure(op.num_bits()); + scratch.configure(op.num_bits(), record_capacity); scratch.ensure(pr.nq_total); scratch.begin_batch(); pr.mono = std::span(scratch.data(), pr.nq_total); @@ -179,7 +180,7 @@ auto resolve_incoming(const std::vector &incoming, // serialized, one VecZ size_t combined_size, // pre-layer op size: bounds the matched set Sink &sink) -> std::vector> { using Resp = typename Sink::Response; - const auto pr = probe_incoming_queries(incoming, op, store, rank_count, sink.stride()); + const auto pr = probe_incoming_queries(incoming, op, store, rank_count, sink.stride(), sink.capacity()); std::vector> responses(rank_count); for (size_t s = 0; s < rank_count; ++s) { responses[s].assign(pr.goff[s + 1] - pr.goff[s], Sink::init_response()); diff --git a/cpp/monoprop/detail/evolution/layer_build/Scan.h b/cpp/monoprop/detail/evolution/layer_build/Scan.h index 6a59b654..1f3195fd 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Scan.h +++ b/cpp/monoprop/detail/evolution/layer_build/Scan.h @@ -30,11 +30,11 @@ #include "monoprop/core/Monomial.h" #include "monoprop/detail/evolution/CutoffContext.h" #include "monoprop/detail/evolution/layer_build/Common.h" +#include "monoprop/detail/evolution/layer_build/TermProduct.h" #include "monoprop/detail/graph_encoding/MPGraphEncodingTypes.h" #include "monoprop/detail/mpi/MPIUtils.h" #include "monoprop/detail/operator/InvertedIndex.h" #include "monoprop/detail/operator/MPOperator.h" -#include "monoprop/detail/operator/RowAccess.h" namespace monoprop::detail { @@ -163,37 +163,11 @@ inline auto rotation_dynamic_gate(std::optional only_rotate_len_k, return true; } -// phase_factor is the basis-specific sign only: Majorana interleave_phase, still to be folded with -// hermitian_phase at emit; Pauli pauli_rotation_sign, already rotation-ready. -// -// `mono` and `new_mono` are scratch owned by the caller for the whole gate, not locals: both are -// overwritten whole here, and a term must not pay for constructing them. Constructing a Bitset means -// deriving the word count from a runtime width, testing it against the inline capacity and -- above -// that capacity -- allocating; with a compile-time width all of that used to fold away to nothing, so -// the scratch is what keeps it off the per-term path. -template -[[gnu::always_inline]] inline auto emit_term_products(const auto &ham, - size_t i, - const typename A::GenContext &ctx, - MonomialLike auto &mono, - MonomialLike auto &new_mono, - size_t &overlap, - int &phase_factor) -> void { - const auto &gen = A::generator(ctx); - // for_each_row_position only sets bits, so the scratch has to start clear; reset() keeps the width, - // unlike assigning a default-constructed Bitset. - mono.reset(); - for_each_row_position(ham, i, [&mono](size_t pos) { mono.set(pos); }); - // One pass instead of two: mono ^ gen and popcount(mono & gen) are both always needed here, so - // fused_xor_into() computes them together, straight into new_mono (see Bitset::fused_xor_into). - // Its result_count (popcount of the XOR) goes unused -- the caller already has new_pop for free - // via mono_pop + gen_pop - 2*overlap -- so it is not threaded through here. - overlap = mono.fused_xor_into(gen, new_mono).overlap; - phase_factor = A::rotation_sign(ctx, mono, new_mono); -} - struct FusedScanResult { - std::vector cos_blocks; // ascending, disjoint, chunk order + std::vector cos_blocks; // ascending, disjoint, chunk order + // No escape tails here: they are folded into the query buffers before this is returned, so a + // consumer only ever sees the finished streams. They live as locals of the scan for exactly as long + // as they are pushed to. std::vector leader_queries; // size R: serialized leader queries per owner rank std::vector> leader_src; // size R: parallel to leader_queries (source op idx) std::vector follower_queries; // size R: serialized follower queries per owner rank @@ -212,10 +186,16 @@ struct FusedScanResult { // op, store, gen and cutoff_eval are deduced from their argument types; no width is named anywhere below // any more -- the monomials this builds take theirs from `gen`, which is the operator's storage width. // -// `store` is a separate argument rather than reached through `op`: build_layer has already bound the -// backend, and re-entering that dispatch here would put a branch on the per-term path, which is the one -// place it cannot go. Rows are read through the accessors, so both backends answer the same scan. -template +// `store` is a separate argument rather than reached through `op`, and its concrete type is what selects +// the per-term kernel: build_layer has already bound the backend, and re-entering that dispatch here +// would put a branch on the per-term path, which is the one place it cannot go. +// +// W is the storage word count, bound by build_layer through with_kernel_width for the same reason and at +// the same seam as the backend and the algebra; 0 means "not specialized" (see TermProductsFor). It is a +// template parameter of the scan rather than something the kernel is handed, so that the per-term code +// below stays an ordinary function body: wrapping it in a generic lambda instead measured 2-3% slower +// even on the unspecialized arm, which does no different work. +template auto fused_find_and_collect(const auto &op, const Store &store, const MonomialLike auto &gen, @@ -239,6 +219,10 @@ auto fused_find_and_collect(const auto &op, res.leader_src.assign(rank_count, std::vector{}); res.follower_queries.assign(rank_count, query_buffer()); res.follower_src.assign(rank_count, std::vector{}); + // Drained into the query buffers by append_escape_tail before this function returns, so they are + // locals; nothing outside this scan ever sees an unfinished stream. + std::vector leader_escapes(rank_count); + std::vector follower_escapes(rank_count); // Sized to R even on the early-return paths below so the fused engine's per-rank src_val_r access // is always in bounds (parallel to leader_src / follower_src). if (capture_values) { @@ -294,14 +278,16 @@ auto fused_find_and_collect(const auto &op, auto &fq = res.follower_queries; auto &fs = res.follower_src; auto &fv = res.follower_val; + auto &lesc = leader_escapes; + auto &fesc = follower_escapes; - // Per-gate, not per-term: the product scratch is overwritten whole per term, so constructing it - // per term would buy nothing and cost a width derivation, an inline-capacity test and -- above - // that capacity -- a heap allocation, every term. Both take the generator's width, which is the - // operator's storage width, so every word op below stays on Bitset's matched-width path. - const auto gen_ctx = A::make_gen_context(gen); - Bitset mono(gen.size()); - Bitset new_mono(gen.size()); + // Per-gate, not per-term: the kernel's product scratch is overwritten whole per term, so + // constructing it per term would buy nothing and cost a width derivation, an inline-capacity test + // and -- above that capacity -- a heap allocation, every term. It takes the generator's width, + // which is the operator's storage width, so every word op inside stays on Bitset's matched-width + // path. Which kernel this is follows from the store (TermProductsFor); the scan below names no + // representation. + typename TermProductsFor::type products(gen, cutoff_eval); // The dynamic gate runs before the product, so a gate-rejected term computes none. // abs_c/v_src come from the caller's coeff read, not re-read. @@ -309,28 +295,25 @@ auto fused_find_and_collect(const auto &op, if (!rotation_dynamic_gate(only_rotate_len_k, mono_pop, cut_st, abs_c)) { return; } - size_t overlap = 0; - int phase_factor = 0; - emit_term_products(store, i, gen_ctx, mono, new_mono, overlap, phase_factor); + const auto [overlap, phase_factor] = products.product(store, i); // Structural cutoff on the partner M⊕G, unless upper_atol rescues it (CutoffContext::is_above_upper). const size_t new_pop = mono_pop + gen_pop - 2 * overlap; - if (const bool struct_pass = cutoff_eval.passes_with_popcount(new_mono, new_pop); - !struct_pass && !cut_st.is_above_upper(abs_c)) { + if (const bool struct_pass = products.passes(new_pop); !struct_pass && !cut_st.is_above_upper(abs_c)) { return; } const int phase = A::emit_phase(phase_factor, mono_pop, gen_pop, overlap); // Single rank: every partner is self-owned, skip the O(W) hash; multi-rank routes by owner. - const size_t r_prime = (rank_count == 1) ? my_rank : monomial_hash(new_mono) % rank_count; + const size_t r_prime = (rank_count == 1) ? my_rank : products.owner(rank_count); const size_t source = i; if (is_follower) { - query_push(fq[r_prime], new_mono, phase); + products.push(QueryOut{fq[r_prime], fesc[r_prime]}, phase); fs[r_prime].push_back(source); if (capture_values) { fv[r_prime].push_back(v_src); } } else { - query_push(lq[r_prime], new_mono, phase); + products.push(QueryOut{lq[r_prime], lesc[r_prime]}, phase); ls[r_prime].push_back(source); if (capture_values) { lv[r_prime].push_back(v_src); @@ -361,7 +344,9 @@ auto fused_find_and_collect(const auto &op, n_foll); } if (rank_count == 1) { - const size_t record_words = query_words(new_mono.num_words()); + // The record width comes off the kernel, not off the generator: it is a property of the form + // a query is pushed in, which is the kernel's business and not the monomial's. + const size_t record_words = products.record_words(); lq[my_rank].reserve(kQueryHeaderWords + ((n_anti - n_foll) * record_words)); ls[my_rank].reserve(n_anti - n_foll); fq[my_rank].reserve(kQueryHeaderWords + (n_foll * record_words)); @@ -428,6 +413,12 @@ auto fused_find_and_collect(const auto &op, } res.cos_blocks.push_back(cos_b.finish()); } + // The one place a stream is finished. The early returns above are all before the first push, so their + // escape buffers are empty and skipping this is a no-op for them. + for (size_t r = 0; r < rank_count; ++r) { + append_escape_tail(res.leader_queries[r], leader_escapes[r]); + append_escape_tail(res.follower_queries[r], follower_escapes[r]); + } return res; } diff --git a/cpp/monoprop/detail/evolution/layer_build/TermProduct.h b/cpp/monoprop/detail/evolution/layer_build/TermProduct.h new file mode 100644 index 00000000..cc071005 --- /dev/null +++ b/cpp/monoprop/detail/evolution/layer_build/TermProduct.h @@ -0,0 +1,443 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +// The scan's per-term kernel. For one anticommuting term it answers five questions -- the product M(+)G, +// the overlap the emitted phase needs, the basis rotation sign, whether the product survives the +// structural cutoff, and (for a survivor) its owner rank and query record -- and those five are exactly +// what changes when a row stops being a bitset. So they are gathered behind one per-gate object, chosen +// off the store type by TermProductsFor, rather than spread over the scan's emit lambda. +// +// The call sequence per term is product() -> passes() -> owner() -> push(), each reading the product the +// previous left. That is stateful on purpose: the product is per-gate scratch, since a term must not pay +// to construct the storage its product goes into (see the note on the scratch monomials below). + +#include +#include +#include +#include +#include +#include +#include + +#include "monoprop/TypeAliases.h" +#include "monoprop/algebra/Algebra.h" +#include "monoprop/algebra/CodesAlgebra.h" +#include "monoprop/detail/evolution/layer_build/Common.h" +#include "monoprop/detail/operator/OperatorIndex.h" +#include "monoprop/detail/operator/SparseRowStore.h" + +namespace monoprop::detail { + +// phase_factor is the basis-specific sign only: Majorana interleave_phase, still to be folded with +// hermitian_phase at emit; Pauli pauli_rotation_sign, already rotation-ready. +// ham and the two monomials are deduced from their argument types; A stays the sole explicit template +// argument at call sites, same as before. +// +// `mono` and `new_mono` are scratch owned by the caller for the whole gate, not locals: both are +// overwritten whole here, and a term must not pay for constructing them. Constructing a Bitset means +// deriving the word count from a runtime width, testing it against the inline capacity and -- above +// that capacity -- allocating; with a compile-time width all of that used to fold away to nothing, so +// the scratch is what keeps it off the per-term path. +template +[[gnu::always_inline]] inline auto emit_term_products(const auto &ham, + size_t i, + const typename A::GenContext &ctx, + MonomialLike auto &mono, + MonomialLike auto &new_mono, + size_t &overlap, + int &phase_factor) -> void { + const auto &gen = A::generator(ctx); + // for_each_position only sets bits, so the scratch has to start clear; reset() keeps the width, + // unlike assigning a default-constructed Bitset. + mono.reset(); + for_each_row_position(ham, i, [&mono](size_t pos) { mono.set(pos); }); + // One pass instead of two: mono ^ gen and popcount(mono & gen) are both always needed here, so + // fused_xor_into() computes them together, straight into new_mono (see Bitset::fused_xor_into). + // Its result_count (popcount of the XOR) goes unused -- the caller already has new_pop for free + // via mono_pop + gen_pop - 2*overlap -- so it is not threaded through here. + overlap = mono.fused_xor_into(gen, new_mono).overlap; + phase_factor = A::rotation_sign(ctx, mono, new_mono); +} + +// What one term's product yields. The product's popcount is not here: the caller gets it for free as +// mono_pop + gen_pop - 2*overlap, where the product itself would have to count it. +struct TermProduct { + size_t overlap = 0; + int phase_factor = 0; +}; + +// Dense rows -- the representation the engine has always used, and the reference the sparse one below +// must match term for term. +template +class DenseTermProducts { +public: + // gen fixes the width of both scratch monomials, and it is the operator's storage width, so every + // word op below stays on Bitset's matched-width path. cutoff_eval is held by reference: it borrows + // the caller's CutoffFn already, so it outlives no less than this does. + DenseTermProducts(const Bitset &gen, const CutoffEvaluator &cutoff_eval) + : ctx_(A::make_gen_context(gen)), + cutoff_(&cutoff_eval), + mono_(gen.size()), + new_mono_(gen.size()) {} + + template + [[gnu::always_inline]] auto product(const Store &store, size_t i) -> TermProduct { + TermProduct out; + emit_term_products(store, i, ctx_, mono_, new_mono_, out.overlap, out.phase_factor); + return out; + } + + [[nodiscard]] auto passes(size_t new_pop) const -> bool { + return cutoff_->passes_with_popcount(new_mono_, new_pop); + } + // find_rank's expression, minus its n_ranks == 0 guard, which the scan's rank_count > 1 + // short-circuit already covers. Owner routing is monomial_hash everywhere, including that initial + // distribution, so this must not become anything else. + [[nodiscard]] auto owner(size_t rank_count) const -> size_t { return monomial_hash(new_mono_) % rank_count; } + auto push(QueryOut out, int phase) const -> void { query_push(out.records, new_mono_, phase); } + + // Record width for the per-rank query reserves, which run before the first product. + [[nodiscard]] auto record_words() const -> size_t { return query_words(new_mono_.num_words()); } + + // The product monomial, for the sparse emitter's fallback and for the differential tests. + [[nodiscard]] auto product_row() const -> const Bitset & { return new_mono_; } + +private: + typename A::GenContext ctx_; + const CutoffEvaluator *cutoff_; + Bitset mono_; + Bitset new_mono_; +}; + +// Thrown by DenseTermProductsW's constructor: with_kernel_width picks W from the store's row word +// count and the generator's own width is independent of that (a stale generator bank, a basis +// change that resized the store but not the gates, ...), so unlike the per-term word ops in +// Bitset.h -- deliberately assert-only, since Release must keep their loops bare -- this binding +// happens once per gate. A real branch there costs nothing next to the per-term work it guards, so +// it stays a check even in Release rather than silently reading past W words of gen/mono/new_mono. +class KernelWidthMismatch : public std::runtime_error { +public: + using std::runtime_error::runtime_error; +}; + +// The dense kernel with the storage word count bound at compile time, chosen once per gate by the +// scan (see with_kernel_width). Answers exactly what DenseTermProducts answers, in the same order +// and to the same values -- what differs is that every word loop inside has a known trip count and +// every operand's storage pointer is resolved once here instead of on each access. +// +// Why this is a separate class rather than a W parameter on DenseTermProducts: only the hot answers +// are worth specializing, and a fallback is still needed for W outside the inline regime, where the +// kernel's inline-operand precondition does not hold. +// +// The cutoff is specialized only for a length cutoff over the whole register. Both other cases -- a +// support cutoff, or an active window narrower than the storage width -- keep going through the +// evaluator. Not for lack of trying: a support arm folding or_sum the same way was measured and cost +// about 1% everywhere, gaining nothing even on the Pauli models that use it, because their per-term +// time is not in the cutoff. A narrow window would need a third kernel, and getting its shift wrong +// would silently change which terms survive. +template +class DenseTermProductsW { +public: + DenseTermProductsW(const Bitset &gen, const CutoffEvaluator &cutoff_eval) + : ctx_(A::make_gen_context(gen)), + cutoff_(&cutoff_eval), + mono_(gen.size()), + new_mono_(gen.size()), + gen_words_(A::generator(ctx_).data()), + mono_words_(mono_.data()), + new_words_(new_mono_.data()) { + if (gen.num_words() != W) { + throw KernelWidthMismatch( + std::format("DenseTermProductsW bound against a generator of {} words; the kernel's W must " + "be the generator's word count.", + W, + gen.num_words())); + } + if (const auto *length = cutoff_eval.length_cutoff(); length != nullptr && length->masks.whole_register()) { + length_cutoff_ = length->cutoff; + } + } + + // The word pointers below point into this object's own bitsets, so a copy would leave the copy + // reading and writing the original's storage. Nothing copies this -- it is the scan's per-gate + // local -- so the case is made unrepresentable rather than documented. + DenseTermProductsW(const DenseTermProductsW &) = delete; + auto operator=(const DenseTermProductsW &) -> DenseTermProductsW & = delete; + DenseTermProductsW(DenseTermProductsW &&) = delete; + auto operator=(DenseTermProductsW &&) -> DenseTermProductsW & = delete; + + template + [[gnu::always_inline]] auto product(const Store &store, size_t i) -> TermProduct { + WordKernel::clear(mono_words_); + // Straight to the words: Bitset::set would re-select the storage pointer for every set bit, + // and a row carries one per surviving slot. + for_each_row_position(store, i, [this](size_t pos) { + mono_words_[pos / Bitset::word_width] |= uint64_t{1} << (pos % Bitset::word_width); + }); + const auto counts = WordKernel::fused_xor_into(mono_words_, gen_words_, new_words_); + return {counts.overlap, A::template rotation_sign_words(ctx_, mono_words_, new_words_)}; + } + + [[nodiscard]] auto passes(size_t new_pop) const -> bool { + if (length_cutoff_.has_value()) { + // Same two clauses as CutoffEvaluator::passes_with_popcount for a length cutoff, in the + // same order: the popcount test proves keep without reading the monomial, and the paired + // test is the xor_sum == 0 clause that rescues a fully paired term of any length. + return new_pop <= *length_cutoff_ || fully_paired_words(new_words_); + } + return cutoff_->passes_with_popcount(new_mono_, new_pop); + } + [[nodiscard]] auto owner(size_t rank_count) const -> size_t { + return WordKernel::splitmix(new_words_) % rank_count; + } + auto push(QueryOut out, int phase) const -> void { query_push(out.records, new_mono_, phase); } + + [[nodiscard]] auto record_words() const -> size_t { return query_words(W); } + [[nodiscard]] auto product_row() const -> const Bitset & { return new_mono_; } + + // Whether passes() answers off the words or through the evaluator. Exists for the differential + // tests, which otherwise cannot tell a run that exercised the word cutoff from one that compared + // the evaluator against itself -- the same reason SparseTermProducts::fell_back() is observable. + [[nodiscard]] auto uses_word_cutoff() const -> bool { return length_cutoff_.has_value(); } + +private: + typename A::GenContext ctx_; + const CutoffEvaluator *cutoff_; + Bitset mono_; + Bitset new_mono_; + // Resolved once, in the constructor's order: each is the data() of a member above, which does not + // move for this object's lifetime because no member below is ever resized or reassigned. + const uint64_t *gen_words_; + uint64_t *mono_words_; + uint64_t *new_words_; + std::optional length_cutoff_ = std::nullopt; +}; + +// The slot capacity a support-form query record is cut to, which is also the scan's scratch product +// capacity -- a query carries exactly such a product, so one number has to serve both or a product that +// fits the scratch would not fit the record. A product occupies at most the cutoff's mode bound plus the +// generator's own modes; an absent bound means the cutoff has no codes form and no term will reach the +// toggle, but the row is sized anyway so the capacity is never zero. +// +// Every rank derives this from the same circuit and cutoff, which is what lets it fix a wire stride with no +// communication -- the same agreement find_rank already needs for the hash width. +[[nodiscard]] inline auto sparse_record_capacity(const Bitset &gen, const CutoffEvaluator &cutoff_eval) -> size_t { + return SparseRowStore::scratch_slots_for(cutoff_eval.max_mode_bound().value_or(SparseRowStore::kMaxSlots), + occupied_mode_count(gen)); +} + +// Support-form rows. The five answers split three ways: product, overlap and rotation sign come off the +// codes word and the two mode lists (sparse_toggle plus A::codes_rotation_sign, O(slots) where the dense +// form is O(storage words)); the cutoff is a popcount test that reads the row only when the bound is +// exceeded; and push() writes the row itself, so a term that survives touches a storage word only when it +// escaped the sparse form or when owner() needs the dense hash at R>1. A term the cutoff rejects touches +// none at all, which is the whole point. +// +// Three cases fall back to the dense kernel for that term, none of them rare enough to assert away: a +// spilled store row (no view exists), a product past the scratch capacity (sparse_toggle reports it +// rather than truncating), and a generator too wide for one codes word. A cutoff that is neither of the +// two concrete functors has no codes form either, and falls back for every term of the gate. +template +class SparseTermProducts { +public: + SparseTermProducts(const Bitset &gen, const CutoffEvaluator &cutoff_eval) + : fallback_(gen, cutoff_eval), + dense_(gen.size()) { + gen_lanes_.reserve(SparseRowStore::kMaxSlots); + bool gen_fits = true; + for_each_mode_slot(gen, [this, &gen_fits](size_t mode, unsigned int code) { + if (gen_lanes_.size() == SparseRowStore::kMaxSlots) { + gen_fits = false; + return; + } + gen_codes_ |= static_cast(code) << (2 * gen_lanes_.size()); + gen_lanes_.push_back(static_cast(mode)); + }); + + // Which cutoff, its bound, and the inactive-mode prefix, all fixed for the propagator's + // lifetime. active_bit_offset counts physical bits and a mode spans two, hence the halving. + if (const auto *length = cutoff_eval.length_cutoff(); length != nullptr) { + kind_ = Kind::Length; + cutoff_value_ = length->cutoff; + inactive_mode_prefix_ = length->masks.active_bit_offset / 2; + } + else if (const auto *support = cutoff_eval.support_cutoff(); support != nullptr) { + kind_ = Kind::Support; + cutoff_value_ = support->cutoff; + inactive_mode_prefix_ = support->masks.active_bit_offset / 2; + } + sparse_usable_ = gen_fits && kind_ != Kind::None; + // Shared with the record stride rather than derived here: a product that fits the scratch has to fit + // the record it is pushed into. + capacity_ = sparse_record_capacity(gen, cutoff_eval); + out_lanes_.resize(capacity_); + } + + template + [[gnu::always_inline]] auto product(const Store &store, size_t i) -> TermProduct { + dense_valid_ = false; + if (sparse_usable_ && !store.spilled(i)) { + const SparseRow mono = store.view(i); + const auto toggled = sparse_toggle(mono, generator(), std::span(out_lanes_.data(), capacity_)); + if (!toggled.overflowed) { + fallback_used_ = false; + product_ = toggled; + return {toggled.overlap, A::codes_rotation_sign(mono, generator())}; + } + } + fallback_used_ = true; + return fallback_.product(store, i); + } + + [[nodiscard]] auto passes(size_t new_pop) const -> bool { + if (fallback_used_) { + return fallback_.passes(new_pop); + } + if (kind_ == Kind::Length) { + return codes_length_passes_with_popcount(product_row(), cutoff_value_, new_pop, inactive_mode_prefix_); + } + return codes_support_passes_with_popcount(product_row(), cutoff_value_, new_pop, inactive_mode_prefix_); + } + // Still the dense hash, and it has to be: owner routing is monomial_hash everywhere, including + // find_rank's initial distribution, and the store's own probe hash is a different function for a + // different purpose. So a *multi-rank* run still materializes once per surviving term here -- moving + // that would mean changing find_rank too. A serial run never calls this (the scan short-circuits at + // rank_count == 1), so it materializes only for the terms that escape. + auto owner(size_t rank_count) -> size_t { return monomial_hash(dense_row()) % rank_count; } + + // The row when there is one, and the dense escape when there is not. Which of the two is not a tuning + // choice: a fully paired product escapes the cutoff, so nothing bounds a query's support. + auto push(QueryOut out, int phase) -> void { + if (fallback_used_) { + sparse_query_push_escape(out.records, out.escapes, dense_row(), capacity_, phase); + return; + } + sparse_query_push(out.records, product_row(), capacity_, phase); + } + [[nodiscard]] auto record_words() const -> size_t { return query_words(sparse_payload_words(capacity_)); } + + // The product in support form. Meaningless when the term fell back to the dense kernel. + [[nodiscard]] auto product_row() const -> SparseRow { return SparseRow{out_lanes_.data(), product_.codes}; } + [[nodiscard]] auto fell_back() const -> bool { return fallback_used_; } + // The record's lane capacity, which is also this kernel's scratch capacity -- see sparse_record_capacity. + [[nodiscard]] auto record_capacity() const -> size_t { return capacity_; } + +private: + enum class Kind : uint8_t { None, Length, Support }; + + [[nodiscard]] auto generator() const -> SparseRow { return SparseRow{gen_lanes_.data(), gen_codes_}; } + + // Memoized because owner() and push() both want it and only push() runs unconditionally. + auto dense_row() -> const Bitset & { + if (fallback_used_) { + return fallback_.product_row(); + } + if (!dense_valid_) { + dense_.reset(); // as in the dense kernel: the slot walk only sets bits + fill_from_sparse_row(product_row(), dense_); + dense_valid_ = true; + } + return dense_; + } + + DenseTermProducts fallback_; + std::vector gen_lanes_ = {}; + RowCodes gen_codes_ = 0; + DefaultInitVector out_lanes_ = {}; + SparseProduct product_ = {}; + Bitset dense_; + size_t capacity_ = 0; + size_t inactive_mode_prefix_ = 0; + unsigned int cutoff_value_ = 0; + Kind kind_ = Kind::None; + bool sparse_usable_ = false; + bool fallback_used_ = true; + bool dense_valid_ = false; +}; + +// Which kernel a store wants at a given bound width. Explicit specializations for the same reason +// QueryKeysFor has them: a store must not know what the scan does with its rows, and an unhandled +// store must fail to compile rather than pick a default. W == 0 is the unspecialized arm, and it is +// the only one the sparse store has -- a sparse row's per-term work is O(slots), not O(storage words), +// so is deliberately left incomplete rather than silently dense. +template +struct TermProductsFor; +template +struct TermProductsFor { + using type = DenseTermProducts; +}; +template +struct TermProductsFor { + using type = DenseTermProductsW; +}; +template +struct TermProductsFor { + using type = SparseTermProducts; +}; + +// Bind the storage word count once per gate and let the scan build a kernel that knows it. +// +// Like with_algebra and with_store, this turns one runtime gate-wide property into a compile-time one +// at a single seam. Everything after that seam is templated on it. Doing this once per gate keeps +// compile time under control; encoding Bitset width in the type was too expensive. +// +// We use tag dispatch instead of passing in a kernel object. The body needs to declare the kernel as a +// local, and passing by reference was ~3% slower on the unspecialized arm because locals optimize better. +// +// Only W in [1, kNarrowKernelWords] is specialized. Above that, the runtime loop is already as fast or +// faster, and fewer code paths are better for instruction cache. +// +// This cap is fixed here, not a build option. Unlike monoprop_SPARSE_ROW_MIN_MODES (ISA-dependent), +// this is a storage-word count and maps to the same width regime on every machine. +// +// kNarrowKernelWords = 4 means 128 storage modes, which covers all current shipped models. Measured +// benefit is about 10% at 2-4 words, fading by 7 words, for about 11% extra `.text` from the four +// specializations. Raising the cap to Bitset::kInlineWords is correct but not faster; setting it to 0 +// restores the pre-seam code path for re-measurement. +inline constexpr size_t kNarrowKernelWords = 4; +static_assert(kNarrowKernelWords <= Bitset::kInlineWords, + "a specialized kernel assumes its operands are inline; above kInlineWords they spill"); + +// Whether a store's per-term work has a storage-word trip count for W to bind at all. A trait rather +// than an is_same_v in the dispatch, so a new backend states its own answer next to the kernel it +// asks for above, instead of silently inheriting "no". +template +inline constexpr bool kBindsKernelWidth = false; +template <> +inline constexpr bool kBindsKernelWidth = true; + +// W as the dispatch will actually pass it: itself while the build specializes that width, 0 once it +// does not. Mapping the arm rather than shortening the dispatch is what makes raising the cap a +// one-line change with no new arm. +template +inline constexpr size_t kCappedKernelWidth = W <= kNarrowKernelWords ? W : 0; + +template +[[gnu::always_inline]] inline auto with_kernel_width(size_t num_words, F &&f) -> decltype(auto) { + // The width regime, not the arm count, is what with_nwords' contract asks a caller to gate on: + // above kInlineWords the words are on the heap, so there is no width to bind even in principle -- + // the kernel's precondition is that every operand is inline. + if constexpr (kBindsKernelWidth) { + if (num_words >= 1 && num_words <= Bitset::kInlineWords) { + return with_nwords(num_words, [&f](std::integral_constant) -> decltype(auto) { + return f(std::integral_constant>{}); + }); + } + } + return f(std::integral_constant{}); +} + +} // namespace monoprop::detail diff --git a/cpp/monoprop/detail/operator/SparseRowStore.h b/cpp/monoprop/detail/operator/SparseRowStore.h index 6303f9f2..50a67b83 100644 --- a/cpp/monoprop/detail/operator/SparseRowStore.h +++ b/cpp/monoprop/detail/operator/SparseRowStore.h @@ -130,7 +130,8 @@ inline auto for_each_mode_slot(const Bitset &mono, Fn &&fn) -> void { } // Occupied modes in a dense monomial, via the same slot walk as sparse_row_hash/dense_row_equals below -- -// what a spilled row's occupied_modes() reports, and what a per-gate generator's mode count also needs. +// what a spilled row's occupied_modes() reports, and what a per-gate generator's mode count also needs +// (see sparse_record_capacity in layer_build/TermProduct.h). [[nodiscard]] inline auto occupied_mode_count(const Bitset &mono) -> size_t { size_t n = 0; for_each_mode_slot(mono, [&n](size_t, unsigned int) { ++n; }); @@ -202,7 +203,7 @@ class SparseRowHasher { // Writes a sparse row's occupied slots into `mono`, which must already be at the row's width and // cleared -- a fresh Bitset(num_bits), or one a caller reset itself before refilling it. The shared body -// behind every dense materialization of a SparseRow below. +// behind every dense materialization of a SparseRow below and in layer_build/Common.h and TermProduct.h. inline auto fill_from_sparse_row(const SparseRow &row, Bitset &mono) -> void { const size_t n = row.num_slots(); for (size_t j = 0; j < n; ++j) { diff --git a/cpp/tests/InlineWidths.h b/cpp/tests/InlineWidths.h new file mode 100644 index 00000000..dab4da91 --- /dev/null +++ b/cpp/tests/InlineWidths.h @@ -0,0 +1,43 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include +#include + +#include "monoprop/Bitset.h" + +namespace test_utils { + +// Runs `body` once per compile-time width in [1, Bitset::kInlineWords] (32..256 inline modes). +// +// kNarrowKernelWords controls which widths pick a specialized scan kernel. The kernels themselves must +// work for the full inline range. We stop at kInlineWords because larger widths spill to heap storage, +// which word kernels do not allow. +// +// Keep this width sweep defined in one place. Multiple test files rely on the same range, and duplicating +// this fold could let one copy drift without notice. +template +auto for_each_inline_width(std::index_sequence, auto &&body) -> void { + // +1 because W == 0 is not a width any word kernel accepts. + (body(std::integral_constant{}), ...); +} + +auto for_each_inline_width(auto &&body) -> void { + for_each_inline_width(std::make_index_sequence{}, + std::forward(body)); +} + +} // namespace test_utils diff --git a/cpp/tests/README.md b/cpp/tests/README.md index daa1a631..30b113d1 100644 --- a/cpp/tests/README.md +++ b/cpp/tests/README.md @@ -14,9 +14,9 @@ discovers every Boost.Test case and registers it twice: registered when an MPI launcher is detected. Both `serial` and `mpi` also get a `sparse-rows` variant per case/rank count -(`monoprop_ROW_STORE=sparse` forced), because the sparse backend keys its rows -differently -- so it reaches the resolve path with a different insert order -- -and it is the backend wide systems resolve to. The MPI sparse-rows +(`monoprop_ROW_STORE=sparse` forced), because the sparse backend has its own +wire format (`query_payload_words_for`'s stride, the escape tail) with no dense +counterpart, and it is the backend wide systems resolve to. The MPI sparse-rows ranks are a separate list, `monoprop_MPI_SPARSE_ROWS_TEST_PROCS` (default `2`), kept independent of `monoprop_MPI_TEST_PROCS` so widening dense rank coverage does not silently multiply how many sparse-row `mpiexec` launches CI pays for. diff --git a/cpp/tests/bitset_tests.cpp b/cpp/tests/bitset_tests.cpp index e14c2135..f42d49c5 100644 --- a/cpp/tests/bitset_tests.cpp +++ b/cpp/tests/bitset_tests.cpp @@ -84,7 +84,7 @@ BOOST_AUTO_TEST_CASE(bitset_count_and_parity_and_cross_word) { } // fused_xor must agree with the composed operator^ / count_and it replaces in the hot path -// (Scan.h's emit_term_products) -- same operands, same three quantities, one pass instead of two. +// (TermProduct.h's emit_term_products) -- same operands, same three quantities, one pass instead of two. BOOST_AUTO_TEST_CASE(bitset_fused_xor_matches_composed_ops) { auto [a, ra] = make_pair<192>({1, 63, 64, 130, 191}); auto [b, rb] = make_pair<192>({63, 64, 65, 130}); diff --git a/cpp/tests/boostAddTests.cmake b/cpp/tests/boostAddTests.cmake index 4cb0f530..7faaae49 100644 --- a/cpp/tests/boostAddTests.cmake +++ b/cpp/tests/boostAddTests.cmake @@ -295,11 +295,12 @@ if(TEST_ENABLE_MPI_VARIANTS AND MPIEXEC_EXECUTABLE) ) endforeach() - # MPI counterpart of the "_sparse_rows" serial variant above: the sparse backend keys its rows - # differently, so it reaches the resolve path with a different row set and a different insert - # order, and it is the backend wide (MPI-scale) systems actually resolve to -- so it needs its - # own multi-rank coverage, not just the single-rank one above. Runs over _mpi_sparse_ranks, not - # _mpi_ranks, so it stays cheap by default regardless of how wide the dense rank list grows. + # MPI counterpart of the "_sparse_rows" serial variant above: the sparse backend is not just an + # alternate local layout, it changes the wire format (query_payload_words_for's per-backend + # stride, the escape tail, kOverflowLane, append_escape_tail have no dense counterpart), and it + # is the backend wide (MPI-scale) systems actually resolve to -- so it needs its own multi-rank + # coverage, not just the single-rank one above. Runs over _mpi_sparse_ranks, not _mpi_ranks, so + # it stays cheap by default regardless of how wide the dense rank list grows. foreach(_mpi_rank IN LISTS _mpi_sparse_ranks) set(mpi_cmd "${MPIEXEC_EXECUTABLE}") list( diff --git a/cpp/tests/evolution_detail_tests.cpp b/cpp/tests/evolution_detail_tests.cpp index ff4e6a35..db06c66c 100644 --- a/cpp/tests/evolution_detail_tests.cpp +++ b/cpp/tests/evolution_detail_tests.cpp @@ -146,6 +146,7 @@ BOOST_AUTO_TEST_CASE(self_resolve_mark_bounded_by_combined_size) { matched.begin_gate(op.size()); op.with_store([&](auto &store) { + const size_t capacity = 0; // no generator here, so no sparse record to size detail::LayerBuildEngine> eng( op, store, @@ -154,7 +155,8 @@ BOOST_AUTO_TEST_CASE(self_resolve_mark_bounded_by_combined_size) { /*my_rank_=*/0, matched, combined_size, - detail::query_payload_words_for(store), + detail::query_payload_words_for(store, capacity), + capacity, RecordingSink{}); eng.queries_r[0] = detail::query_buffer(); detail::query_push(eng.queries_r[0], terms[1], 1); diff --git a/cpp/tests/sparse_wire_tests.cpp b/cpp/tests/sparse_wire_tests.cpp new file mode 100644 index 00000000..b67a84b7 --- /dev/null +++ b/cpp/tests/sparse_wire_tests.cpp @@ -0,0 +1,534 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// The support-form query record: a row must survive the wire exactly, and the record must keep the dense +// record's shape so the stride arithmetic, alltoallv counts and phase/value readers work unchanged. +// +// Plus the two things that shape carries with it -- the header word every buffer opens with, and the tail a +// query too wide for any fixed-stride sparse record escapes to. + +#include + +#include +#include +#include +#include + +#include "monoprop/TypeAliases.h" +#include "monoprop/detail/evolution/layer_build/Common.h" +#include "monoprop/detail/operator/OperatorIndex.h" +#include "monoprop/detail/operator/SparseRowStore.h" + +using namespace monoprop; +using namespace monoprop::detail; + +namespace { + +struct OwnedRow { + std::vector lanes; + RowCodes codes = 0; + + [[nodiscard]] auto view() const -> SparseRow { return SparseRow{lanes.data(), codes}; } +}; + +// The four batch/codec cases below key their assertions off a per-record shape, so a default-constructed +// placeholder stands in for the escaped records' absent row. +static_assert(std::is_default_constructible_v); + +auto random_row(std::mt19937_64 &rng, size_t num_modes, size_t capacity) -> OwnedRow { + OwnedRow row{std::vector(capacity, 0), 0}; + // Ascending distinct modes, which is what a real row is; a random unsorted list would not be one. + std::vector modes; + for (size_t m = 0; m < num_modes && modes.size() < capacity; ++m) { + if ((rng() % 4) == 0) { + modes.push_back(m); + } + } + for (size_t j = 0; j < modes.size(); ++j) { + row.lanes[j] = static_cast(modes[j]); + row.codes |= static_cast(1U + (rng() % 3U)) << (2 * j); + } + return row; +} + +} // namespace + +// A record is lane words, then codes, then phase -- so the phase sits at offset `payload`, exactly where +// the dense reader expects it for a payload of that many words. That is the property that lets the record +// machinery stay untouched, so it is pinned rather than left implicit. +BOOST_AUTO_TEST_CASE(sparse_wire_record_keeps_the_dense_record_shape) { + BOOST_TEST(sparse_lane_words(1U) == 1U); + BOOST_TEST(sparse_lane_words(4U) == 1U); + BOOST_TEST(sparse_lane_words(5U) == 2U); + BOOST_TEST(sparse_lane_words(12U) == 3U); + BOOST_TEST(sparse_lane_words(32U) == 8U); + BOOST_TEST(sparse_payload_words(12U) == 4U); + BOOST_TEST(query_words(sparse_payload_words(12U)) == 5U); + BOOST_TEST(query_words_fused(sparse_payload_words(12U)) == 6U); + + // A 12-slot row rides 5 words where a 1024-mode monomial needs 33, and the two are equal at 128 modes + // -- which is roughly where the crossover puts the sparse backend in a wheel build anyway. + BOOST_TEST(query_words(sparse_payload_words(12U)) < query_words(1024U / 32U)); + BOOST_TEST(query_words(sparse_payload_words(12U)) == query_words(128U / 32U)); + + constexpr size_t kCapacity = 12; + VecZ buf = query_buffer(); + OwnedRow row{std::vector(kCapacity, 0), 0}; + row.lanes[0] = 7; + row.codes = 0b11ULL; + sparse_query_push(buf, row.view(), kCapacity, -1); + BOOST_REQUIRE(buf.size() == kQueryHeaderWords + query_words(sparse_payload_words(kCapacity))); + BOOST_TEST(query_record_count(buf) == 1U); + // query_phase reads the phase off the payload width alone, with no idea what the payload holds. + BOOST_TEST(query_phase(buf, 0, sparse_payload_words(kCapacity)) == -1); +} + +BOOST_AUTO_TEST_CASE(sparse_wire_round_trips_rows_exactly) { + std::mt19937_64 rng(20260812U); + size_t full_rows = 0; + size_t empty_rows = 0; + for (const size_t capacity : {1U, 4U, 5U, 12U, 32U}) { + for (const size_t num_modes : {32U, 64U, 1024U, 32000U}) { + const size_t stride = query_words(sparse_payload_words(capacity)); + VecZ buf = query_buffer(); + std::vector rows; + for (size_t t = 0; t < 40; ++t) { + // Every fifth record is the empty row: at these mode counts the random generator would + // essentially never produce one, and an empty row is the record whose lanes are all + // padding. + rows.push_back((t % 5) == 0 ? OwnedRow{std::vector(capacity, 0), 0} + : random_row(rng, num_modes, capacity)); + // Phases are +-1 on the real path; both must survive the unsigned round-trip. + sparse_query_push(buf, rows.back().view(), capacity, (t % 2) == 0 ? 1 : -1); + } + BOOST_REQUIRE(buf.size() == kQueryHeaderWords + (rows.size() * stride)); + BOOST_REQUIRE(query_record_count(buf) == rows.size()); + + for (size_t q = 0; q < rows.size(); ++q) { + std::vector lanes(capacity, 0xEEEE); // poisoned, so an unwritten lane shows up + RowCodes codes = 0; + int phase = 0; + sparse_query_read(buf, q, stride, capacity, lanes.data(), codes, phase); + BOOST_TEST(codes == rows[q].codes); + BOOST_TEST(phase == ((q % 2) == 0 ? 1 : -1)); + const size_t n = row_slot_count(codes); + BOOST_REQUIRE(n == rows[q].view().num_slots()); + for (size_t j = 0; j < n; ++j) { + BOOST_TEST(lanes[j] == rows[q].lanes[j]); + } + full_rows += n == capacity ? 1 : 0; + empty_rows += n == 0 ? 1 : 0; + } + } + } + // Both edges of the capacity have to have occurred, or the packing was never pushed to its bounds. + BOOST_TEST(full_rows > 0U); + BOOST_TEST(empty_rows > 0U); +} + +// Mode indices are packed four to a word, so a lane must not bleed into its neighbours. The largest mode +// a store admits is kMaxModes - 1, which is also the widest lane value. +BOOST_AUTO_TEST_CASE(sparse_wire_packs_lanes_without_bleeding) { + constexpr size_t kCapacity = 8; // two lane words, so the boundary between them is exercised + const size_t stride = query_words(sparse_payload_words(kCapacity)); + const auto top = static_cast(SparseRowStore::kMaxModes - 1); + + OwnedRow row{std::vector(kCapacity, 0), 0}; + // Ascending, and straddling the 4-lane word boundary with extreme values on both sides of it. + const std::vector modes{0, 1, top - 2, top - 1, top, 0, 0, 0}; + for (size_t j = 0; j < 5; ++j) { + row.lanes[j] = modes[j]; + row.codes |= RowCodes{0b11} << (2 * j); + } + + VecZ buf = query_buffer(); + sparse_query_push(buf, row.view(), kCapacity, 1); + std::vector lanes(kCapacity, 0); + RowCodes codes = 0; + int phase = 0; + sparse_query_read(buf, 0, stride, kCapacity, lanes.data(), codes, phase); + BOOST_TEST(codes == row.codes); + for (size_t j = 0; j < 5; ++j) { + BOOST_TEST(lanes[j] == modes[j]); + } +} + +// Records are read by position out of one flat buffer, so a short row must still occupy a full stride and +// must not be able to see the previous record's lanes. +BOOST_AUTO_TEST_CASE(sparse_wire_short_rows_keep_the_full_stride) { + constexpr size_t kCapacity = 12; + const size_t stride = query_words(sparse_payload_words(kCapacity)); + VecZ buf = query_buffer(); + + OwnedRow full{std::vector(kCapacity, 0), 0}; + for (size_t j = 0; j < kCapacity; ++j) { + full.lanes[j] = static_cast(100 + j); + full.codes |= RowCodes{0b11} << (2 * j); + } + OwnedRow empty{std::vector(kCapacity, 0), 0}; + + sparse_query_push(buf, full.view(), kCapacity, 1); + sparse_query_push(buf, empty.view(), kCapacity, -1); + BOOST_REQUIRE(buf.size() == kQueryHeaderWords + (2 * stride)); + + std::vector lanes(kCapacity, 0xEEEE); + RowCodes codes = 0xDEAD; + int phase = 0; + sparse_query_read(buf, 1, stride, kCapacity, lanes.data(), codes, phase); + BOOST_TEST(codes == 0U); + BOOST_TEST(phase == -1); + BOOST_TEST(row_slot_count(codes) == 0U); + // Nothing was written into lanes, so the poison is still there -- the empty row cannot have inherited + // the previous record's modes. + BOOST_TEST(lanes[0] == 0xEEEE); +} + +// --- the query key batches ------------------------------------------------------------------------- +// +// Both resolve paths fill one of these from wire records and hand it to find_batch contiguously. They are +// grow-only and reused across layers, which is the shape the measurements asked for, and that reuse is +// where the two ways to get it wrong live: a stale element read before being overwritten, and -- for the +// sparse batch, whose keys are views -- a view left pointing into storage that growth moved. + +// Each store is queried in the form it keys its rows by, so a resolve never converts one into the other. +static_assert(std::is_same_v::type, DenseQueryKeys>); +static_assert(std::is_same_v::type, SparseQueryKeys>); + +BOOST_AUTO_TEST_CASE(query_keys_dense_batch_round_trips_records) { + constexpr size_t kNumBits = 128; + const size_t stride = query_words(kNumBits / 64); + VecZ buf = query_buffer(); + std::vector monos; + for (size_t t = 0; t < 40; ++t) { + Bitset mono(kNumBits); + mono.set(t); + mono.set(kNumBits - 1 - t); + monos.push_back(mono); + query_push(buf, mono, (t % 2) == 0 ? 1 : -1); + } + + BOOST_TEST(query_record_count(buf) == monos.size()); + DenseQueryKeys keys; + keys.configure(kNumBits, /*capacity=*/0); + keys.ensure(monos.size()); + keys.begin_batch(); + for (size_t q = 0; q < monos.size(); ++q) { + BOOST_TEST(keys.read_record(buf, q, stride, q) == ((q % 2) == 0 ? 1 : -1)); + } + for (size_t q = 0; q < monos.size(); ++q) { + BOOST_TEST((keys[q] == monos[q])); + BOOST_TEST((keys.data()[q] == monos[q])); + } +} + +BOOST_AUTO_TEST_CASE(query_keys_sparse_batch_survives_growth) { + constexpr size_t kCapacity = 12; + const size_t stride = query_words(sparse_payload_words(kCapacity)); + std::mt19937_64 rng(20260812U); + VecZ buf = query_buffer(); + std::vector rows; + for (size_t t = 0; t < 64; ++t) { + rows.push_back(random_row(rng, 1024, kCapacity)); + sparse_query_push(buf, rows.back().view(), kCapacity, (t % 2) == 0 ? 1 : -1); + } + + SparseQueryKeys keys; + keys.configure(/*num_bits=*/2048, kCapacity); + + // Fill a small prefix, then grow: the lane array reallocates, so every view has to be rebuilt. If only + // the new tail were, the prefix read back below would be reading freed storage. + keys.ensure(8); + for (size_t q = 0; q < 8; ++q) { + BOOST_TEST(keys.read_record(buf, q, stride, q) == ((q % 2) == 0 ? 1 : -1)); + } + keys.ensure(rows.size()); + // Checked as a pointer invariant rather than by reading the prefix and hoping it looks wrong: a view + // left over from before the growth points into freed storage, which is undefined behaviour and might + // read back plausibly. Every view must address this batch's current lane array at its own stride. + for (size_t q = 0; q < rows.size(); ++q) { + BOOST_TEST(keys.data()[q].row.modes == keys.data()[0].row.modes + (q * kCapacity)); + } + for (size_t q = 8; q < rows.size(); ++q) { + BOOST_TEST(keys.read_record(buf, q, stride, q) == ((q % 2) == 0 ? 1 : -1)); + } + + for (size_t q = 0; q < 8; ++q) { + // Re-read the prefix records so the prefix slots are written again after the growth; what is being + // checked is that the view still addresses this batch's own lanes. + BOOST_TEST(keys.read_record(buf, q, stride, q) == ((q % 2) == 0 ? 1 : -1)); + } + for (size_t q = 0; q < rows.size(); ++q) { + const auto &key = keys.data()[q]; + BOOST_REQUIRE(!key.is_spilled()); + BOOST_TEST(key.row.codes == rows[q].codes); + const size_t n = key.row.num_slots(); + BOOST_REQUIRE(n == rows[q].view().num_slots()); + for (size_t j = 0; j < n; ++j) { + BOOST_TEST(key.row.mode(j) == rows[q].view().mode(j)); + } + } +} + +// configure() with a different extent must drop the storage: a thread servicing two propagators of +// different widths would otherwise write a wide record into a narrow element. +BOOST_AUTO_TEST_CASE(query_keys_reconfigure_resizes_the_elements) { + DenseQueryKeys dense; + dense.configure(64, 0); + dense.ensure(4); + BOOST_TEST(dense[0].size() == 64U); + dense.configure(256, 0); + dense.ensure(4); + BOOST_TEST(dense[0].size() == 256U); + + SparseQueryKeys sparse; + sparse.configure(64, 4); + sparse.ensure(4); + const auto *narrow_base = sparse.data()[0].row.modes; + const auto *narrow_next = sparse.data()[1].row.modes; + BOOST_TEST(narrow_next - narrow_base == 4); + sparse.configure(64, 12); + sparse.ensure(4); + BOOST_TEST(sparse.data()[1].row.modes - sparse.data()[0].row.modes == 12); +} + +// --- the escape tail --------------------------------------------------------------------------------- +// +// A query is M ⊕ G and a fully paired product escapes the cutoff, so a query's support is unbounded and no +// fixed-stride sparse record can hold every one. The escaped record keeps its place and its stride and +// names a tail entry instead, which is what leaves the engine's offsets, counts and compaction as plain +// arithmetic. + +namespace { + +// The scan's shape: records into one buffer, escape monomials into another, concatenated once pushing is +// done -- a record cannot be appended after the tail has started. +struct WireBuffers { + VecZ records = query_buffer(); + VecZ escapes; + + auto finish() -> VecZ { + VecZ out = records; + out.insert(out.end(), escapes.begin(), escapes.end()); + return out; + } +}; + +auto wide_monomial(size_t num_bits, size_t seed) -> Bitset { + Bitset mono(num_bits); + // Fully paired and far wider than any record capacity below: the shape that escapes the cutoff and so + // cannot be sized away. + for (size_t mode = seed % 3; mode < (num_bits / 2); mode += 3) { + mono.set(2 * mode); + mono.set((2 * mode) + 1); + } + return mono; +} + +} // namespace + +BOOST_AUTO_TEST_CASE(sparse_wire_escaped_records_carry_their_monomial_in_the_tail) { + constexpr size_t kCapacity = 8; + constexpr size_t kNumBits = 512; + const size_t stride = query_words(sparse_payload_words(kCapacity)); + std::mt19937_64 rng(20260813U); + + WireBuffers wire; + std::vector rows; + std::vector escaped; + std::vector is_escape; + for (size_t t = 0; t < 24; ++t) { + // Every third record escapes, so escaped and plain records interleave -- the ordering the tail + // indices have to survive. + if (t % 3 == 2) { + escaped.push_back(wide_monomial(kNumBits, t)); + sparse_query_push_escape(wire.records, wire.escapes, escaped.back(), kCapacity, (t % 2) == 0 ? 1 : -1); + rows.emplace_back(); + is_escape.push_back(true); + continue; + } + rows.push_back(random_row(rng, kNumBits / 2, kCapacity)); + sparse_query_push(wire.records, rows.back().view(), kCapacity, (t % 2) == 0 ? 1 : -1); + is_escape.push_back(false); + } + const VecZ buf = wire.finish(); + + // The header counts records, escaped ones included; the tail sits right after them. + BOOST_REQUIRE(query_record_count(buf) == is_escape.size()); + BOOST_REQUIRE(query_tail_offset(buf, stride) == kQueryHeaderWords + (is_escape.size() * stride)); + BOOST_REQUIRE(buf.size() == query_tail_offset(buf, stride) + (escaped.size() * (kNumBits / 64))); + + SparseQueryKeys keys; + keys.configure(kNumBits, kCapacity); + keys.ensure(is_escape.size()); + keys.begin_batch(); + for (size_t q = 0; q < is_escape.size(); ++q) { + BOOST_TEST(keys.read_record(buf, q, stride, q) == ((q % 2) == 0 ? 1 : -1)); + } + + size_t next_escape = 0; + for (size_t q = 0; q < is_escape.size(); ++q) { + const auto &key = keys[q]; + BOOST_REQUIRE(key.is_spilled() == is_escape[q]); + if (is_escape[q]) { + // The monomial has to arrive bit-for-bit: it is the only form this query exists in. + BOOST_TEST((*key.spilled == escaped[next_escape])); + ++next_escape; + continue; + } + BOOST_TEST(key.row.codes == rows[q].codes); + for (size_t j = 0; j < key.row.num_slots(); ++j) { + BOOST_TEST(key.row.mode(j) == rows[q].view().mode(j)); + } + } + BOOST_TEST(next_escape == escaped.size()); +} + +// The fused sink widens every record by a value word, which moves the tail. The escape indices survive +// because they name a position *within* the tail rather than an offset into the buffer -- so the same +// records read back at the fused stride. +BOOST_AUTO_TEST_CASE(sparse_wire_escape_indices_survive_the_fused_relayout) { + constexpr size_t kCapacity = 6; + constexpr size_t kNumBits = 256; + const size_t payload = sparse_payload_words(kCapacity); + std::mt19937_64 rng(20260814U); + + WireBuffers wire; + std::vector escaped; + std::vector values; + std::vector is_escape; + for (size_t t = 0; t < 12; ++t) { + if (t % 2 == 0) { + escaped.push_back(wide_monomial(kNumBits, t)); + sparse_query_push_escape(wire.records, wire.escapes, escaped.back(), kCapacity, 1); + is_escape.push_back(true); + } + else { + const OwnedRow row = random_row(rng, kNumBits / 2, kCapacity); + sparse_query_push(wire.records, row.view(), kCapacity, -1); + is_escape.push_back(false); + } + values.push_back(static_cast(t) * 0.5 - 3.0); + } + const VecZ plain = wire.finish(); + + VecZ fused; + build_fused_query_value(plain, values, fused, payload); + BOOST_REQUIRE(query_record_count(fused) == is_escape.size()); + BOOST_REQUIRE(fused.size() + == kQueryHeaderWords + (is_escape.size() * query_words_fused(payload)) + + (escaped.size() * (kNumBits / 64))); + + SparseQueryKeys keys; + keys.configure(kNumBits, kCapacity); + keys.ensure(is_escape.size()); + keys.begin_batch(); + size_t next_escape = 0; + for (size_t q = 0; q < is_escape.size(); ++q) { + BOOST_TEST(keys.read_record(fused, q, query_words_fused(payload), q) == (is_escape[q] ? 1 : -1)); + BOOST_TEST(query_value(fused, q, payload) == values[q]); + BOOST_REQUIRE(keys[q].is_spilled() == is_escape[q]); + if (is_escape[q]) { + BOOST_TEST((*keys[q].spilled == escaped[next_escape])); + ++next_escape; + } + } + BOOST_TEST(next_escape == escaped.size()); +} + +// Deferred self-misses are inserted after both resolve passes, by which time the batch has been refilled +// many times over -- so a key that must survive that has to be retained, and retaining it has to copy. +BOOST_AUTO_TEST_CASE(query_keys_retained_survive_a_refill) { + constexpr size_t kCapacity = 6; + constexpr size_t kNumBits = 256; + const size_t stride = query_words(sparse_payload_words(kCapacity)); + std::mt19937_64 rng(20260815U); + + WireBuffers wire; + std::vector rows; + std::vector escaped; + std::vector is_escape; + for (size_t t = 0; t < 8; ++t) { + if (t % 4 == 3) { + escaped.push_back(wide_monomial(kNumBits, t)); + sparse_query_push_escape(wire.records, wire.escapes, escaped.back(), kCapacity, 1); + rows.emplace_back(); + is_escape.push_back(true); + continue; + } + rows.push_back(random_row(rng, kNumBits / 2, kCapacity)); + sparse_query_push(wire.records, rows.back().view(), kCapacity, 1); + is_escape.push_back(false); + } + const VecZ buf = wire.finish(); + + SparseQueryKeys keys; + keys.configure(kNumBits, kCapacity); + keys.ensure(2); + + // Read the records two at a time into the same two slots, retaining every key -- the resolve path's + // batching, with kResolveBatch of 2. + std::vector handles; + for (size_t q = 0; q < is_escape.size(); q += 2) { + keys.begin_batch(); + for (size_t j = 0; j < 2; ++j) { + (void)keys.read_record(buf, q + j, stride, j); + } + for (size_t j = 0; j < 2; ++j) { + handles.push_back(keys.retain(j)); + } + } + + size_t next_escape = 0; + for (size_t q = 0; q < handles.size(); ++q) { + const auto key = keys.retained(handles[q]); + BOOST_REQUIRE(key.is_spilled() == is_escape[q]); + if (is_escape[q]) { + BOOST_TEST((*key.spilled == escaped[next_escape])); + ++next_escape; + continue; + } + BOOST_TEST(key.row.codes == rows[q].codes); + for (size_t j = 0; j < key.row.num_slots(); ++j) { + BOOST_TEST(key.row.mode(j) == rows[q].view().mode(j)); + } + } + BOOST_TEST(next_escape == escaped.size()); + + // The dense batch owes the same guarantee, and its keys are whole monomials. + VecZ dense_buf = query_buffer(); + std::vector monos; + for (size_t t = 0; t < 8; ++t) { + Bitset mono(kNumBits); + mono.set(t); + mono.set(kNumBits - 1 - t); + monos.push_back(mono); + query_push(dense_buf, mono, 1); + } + DenseQueryKeys dense; + dense.configure(kNumBits, 0); + dense.ensure(2); + std::vector dense_handles; + for (size_t q = 0; q < monos.size(); q += 2) { + dense.begin_batch(); + for (size_t j = 0; j < 2; ++j) { + (void)dense.read_record(dense_buf, q + j, query_words(kNumBits / 64), j); + } + for (size_t j = 0; j < 2; ++j) { + dense_handles.push_back(dense.retain(j)); + } + } + for (size_t q = 0; q < dense_handles.size(); ++q) { + BOOST_TEST((dense.retained(dense_handles[q]) == monos[q])); + } +} diff --git a/cpp/tests/term_product_tests.cpp b/cpp/tests/term_product_tests.cpp new file mode 100644 index 00000000..80865b58 --- /dev/null +++ b/cpp/tests/term_product_tests.cpp @@ -0,0 +1,636 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// The scan's per-term kernel, sparse against dense, through the interface the scan actually calls: +// product() -> passes() -> owner() -> push(). codes_algebra_tests.cpp and codes_product_tests.cpp already +// pin the pieces against their dense counterparts; what this adds is the *composition* -- that +// SparseTermProducts routes a term to the right one of them, and that its three fallbacks (a spilled store +// row, a product past the scratch capacity, a cutoff with no codes form) produce the dense answer rather +// than a wrong one. +// +// This is the gate on Stage 6's store swap: with MPOperator::Store still OperatorIndex the sparse kernel +// is unreachable from the library, so nothing else would instantiate it. +// +// The width-bound dense kernel, DenseTermProductsW, is the third answer to the same questions and +// is compared here for the same reason: it restates four of the five off raw words with the storage word +// count fixed at compile time, so nothing but a differential test can tell a restatement that agrees +// from one that merely runs. Its cases live at the bottom of the file; WordKernel's own primitives +// are pinned separately in word_kernel_tests.cpp. + +#include + +#include +#include +#include +#include +#include + +#include "monoprop/TypeAliases.h" +#include "monoprop/algebra/Algebra.h" +#include "monoprop/algebra/AlgebraCommon.h" +#include "monoprop/detail/evolution/layer_build/TermProduct.h" +#include "monoprop/detail/monomial_propagator/MonomialPropagatorCommon.h" +#include "monoprop/detail/operator/OperatorIndex.h" +#include "monoprop/detail/operator/SparseRowStore.h" + +#include "InlineWidths.h" +#include "RandomMonomial.h" +#include "TestData.h" +#include "TestUtilities.h" + +using namespace monoprop; +using namespace monoprop::detail; + +namespace { + +// Which branches a run exercised. Every case asserts on these: a fallback that silently swallowed every +// term would otherwise pass by comparing the dense kernel against itself. +struct Seen { + size_t sparse_terms = 0; + size_t fallback_terms = 0; + size_t cutoff_passed = 0; + size_t cutoff_failed = 0; + size_t sparse_cutoff_failed = 0; // the codes cutoff said no, not the dense one + size_t row_records = 0; // survivors pushed as a support-form record + size_t escaped_records = 0; // survivors pushed as an escape plus a tail entry +}; + +// A fully paired term of `modes` modes. These are the rows that make support unbounded -- xor_sum == 0 +// escapes the cutoff -- so they are what drives a store row to spill and a product to overflow. +auto paired_term(size_t num_modes, size_t modes) -> Bitset { + Bitset mono(2 * num_modes); + for (size_t m = 0; m < modes; ++m) { + mono.set(2 * m); + mono.set((2 * m) + 1); + } + return mono; +} + +// One term against one generator through both kernels, comparing every answer the scan reads. The two +// stores must hold the same monomial at index i, which the caller guarantees by inserting in lockstep. +template +auto check_term(DenseTermProducts &dense, + SparseTermProducts &sparse_kernel, + const OperatorIndex &packed, + const SparseRowStore &sparse_store, + size_t i, + size_t gen_pop, + Seen &seen) -> void { + const size_t mono_pop = packed.popcount(i); + BOOST_REQUIRE(sparse_store.popcount(i) == mono_pop); + + const auto reference = dense.product(packed, i); + const auto candidate = sparse_kernel.product(sparse_store, i); + BOOST_TEST(candidate.overlap == reference.overlap); + BOOST_TEST(candidate.phase_factor == reference.phase_factor); + + const size_t new_pop = mono_pop + gen_pop - (2 * reference.overlap); + const bool reference_passes = dense.passes(new_pop); + const bool candidate_passes = sparse_kernel.passes(new_pop); + BOOST_TEST(candidate_passes == reference_passes); + + // Both remaining answers are read only for a surviving term, so ask them only there. owner() is still + // the dense hash on both sides, because owner routing is monomial_hash everywhere. + if (reference_passes) { + for (const size_t rank_count : {2U, 3U, 8U}) { + BOOST_TEST(sparse_kernel.owner(rank_count) == dense.owner(rank_count)); + } + + // The two records no longer agree byte for byte -- that is what the support form is for -- so what + // is compared is what a resolver reads back out of them: the same key, and the same phase. + const size_t num_bits = packed.num_bits(); + const size_t capacity = sparse_kernel.record_capacity(); + BOOST_TEST(sparse_kernel.record_words() == query_words(sparse_payload_words(capacity))); + + VecZ unused_escapes; + VecZ reference_record = query_buffer(); + dense.push(QueryOut{reference_record, unused_escapes}, -1); + BOOST_TEST(unused_escapes.empty()); // a dense record has nowhere to escape to and never needs one + + VecZ candidate_record = query_buffer(); + VecZ candidate_escapes; + sparse_kernel.push(QueryOut{candidate_record, candidate_escapes}, -1); + const bool escaped = !candidate_escapes.empty(); + BOOST_TEST(escaped == sparse_kernel.fell_back()); + append_escape_tail(candidate_record, candidate_escapes); + seen.row_records += escaped ? 0 : 1; + seen.escaped_records += escaped ? 1 : 0; + + DenseQueryKeys dense_keys; + dense_keys.configure(num_bits, 0); + dense_keys.ensure(1); + dense_keys.begin_batch(); + BOOST_TEST(dense_keys.read_record(reference_record, 0, dense.record_words(), 0) == -1); + + SparseQueryKeys sparse_keys; + sparse_keys.configure(num_bits, capacity); + sparse_keys.ensure(1); + sparse_keys.begin_batch(); + BOOST_TEST(sparse_keys.read_record(candidate_record, 0, sparse_kernel.record_words(), 0) == -1); + BOOST_TEST(sparse_keys[0].is_spilled() == escaped); + BOOST_TEST((key_monomial(sparse_keys[0], num_bits) == dense_keys[0])); + // And the store agrees with the key either way, which is what the resolve's find_batch rests on. + BOOST_TEST(sparse_row_hash(sparse_keys[0]) == sparse_row_hash(dense_keys[0])); + } + + if (sparse_kernel.fell_back()) { + ++seen.fallback_terms; + } + else { + ++seen.sparse_terms; + seen.sparse_cutoff_failed += reference_passes ? 0 : 1; + } + seen.cutoff_passed += reference_passes ? 1 : 0; + seen.cutoff_failed += reference_passes ? 0 : 1; +} + +// Every term of `terms` against every generator of `gens`, over one algebra and one cutoff. +template +auto sweep(const std::vector &terms, + const std::vector &gens, + const CutoffFn &cutoff_fn, + size_t sparse_slots, + Seen &seen) -> void { + BOOST_REQUIRE(!terms.empty()); + const size_t num_bits = terms.front().size(); + OperatorIndex packed(num_bits); + SparseRowStore sparse_store(num_bits, sparse_slots); + for (const auto &mono : terms) { + packed.push_back(mono); + sparse_store.push_back(mono); + } + + const CutoffEvaluator cutoff_eval{cutoff_fn}; + for (const auto &gen : gens) { + DenseTermProducts dense(gen, cutoff_eval); + SparseTermProducts sparse_kernel(gen, cutoff_eval); + const size_t gen_pop = gen.count(); + for (size_t i = 0; i < terms.size(); ++i) { + check_term(dense, sparse_kernel, packed, sparse_store, i, gen_pop, seen); + } + } +} + +} // namespace + +// Randomized terms and generators at three widths, both algebras, both cutoff kinds. Slots are generous +// enough that no row spills, so what is under test is the sparse path itself. +BOOST_AUTO_TEST_CASE(term_product_sparse_kernel_matches_dense_on_randomized_rows) { + std::mt19937_64 rng(20260812U); + Seen seen; + for (const size_t num_modes : {32U, 64U, 300U}) { + std::vector terms; + for (size_t t = 0; t < 150; ++t) { + terms.push_back(test_utils::random_monomial(rng, num_modes, 6)); + } + std::vector gens; + for (size_t t = 0; t < 8; ++t) { + gens.push_back(test_utils::random_monomial(rng, num_modes, 4)); + } + for (const unsigned int cutoff : {4U, 8U}) { + const auto length = cutoff_function(CutoffType::Length, cutoff, num_modes, 2 * num_modes); + const auto support = cutoff_function(CutoffType::Support, cutoff, num_modes, 2 * num_modes); + sweep(terms, gens, length, SparseRowStore::kMaxSlots, seen); + sweep(terms, gens, support, SparseRowStore::kMaxSlots, seen); + sweep(terms, gens, support, SparseRowStore::kMaxSlots, seen); + sweep(terms, gens, length, SparseRowStore::kMaxSlots, seen); + } + } + // Most terms take the sparse path, and some do not: the terms are drawn up to 6 modes wide against a + // cutoff of 4 or 8, and a term above the bound overflows a capacity that is sized from that bound. + // Which is the real shape of the thing -- a stored row exceeds the bound whenever it is fully paired. + BOOST_TEST(seen.sparse_terms > 20000U); + BOOST_TEST(seen.fallback_terms > 0U); + BOOST_TEST(seen.cutoff_passed > 0U); + // The point of the codes cutoff is rejecting a product without materializing it, so a run where + // nothing was rejected would not have tested it. + BOOST_TEST(seen.sparse_cutoff_failed > 0U); + // Every survivor here was pushed as a row: an overflowing product is wider than the cutoff's bound, so + // unless it is fully paired the cutoff rejects it before a record is ever asked for. A *surviving* + // escape needs a fully paired product, which has a case of its own below. + BOOST_TEST(seen.row_records > 0U); + BOOST_TEST(seen.escaped_records == 0U); +} + +// Real terms and generators: the fixtures' Hamiltonian keys against their Majorana generator list, which +// is where the products and the ordering signs are the ones the engine actually computes. +BOOST_AUTO_TEST_CASE(term_product_sparse_kernel_matches_dense_on_fixture_generators) { + Seen seen; + for (const std::string name : {"random_exact.msgpack", "lih_fermionic_spin_exact.msgpack"}) { + const auto data = test_utils::load_case_data(name); + const size_t num_bits = 2 * data.num_modes; + const size_t max_index = 2 * data.num_modes; + + std::vector terms; + for (const auto &[inds, coeff] : data.hamiltonian) { + terms.push_back(indices_to_bitset_checked(inds, max_index, num_bits)); + } + std::vector gens; + for (const auto &inds : data.majoranas) { + gens.push_back(indices_to_bitset_checked(inds, max_index, num_bits)); + } + BOOST_REQUIRE(!terms.empty()); + BOOST_REQUIRE(!gens.empty()); + // A stride, not the whole cross product: the fixtures are large and the pairs are homogeneous. + std::vector sampled_gens; + const size_t stride = gens.size() > 12 ? (gens.size() / 12) + 1 : 1; + for (size_t g = 0; g < gens.size(); g += stride) { + sampled_gens.push_back(gens[g]); + } + + for (const unsigned int cutoff : {4U, 6U}) { + sweep(terms, + sampled_gens, + cutoff_function(CutoffType::Length, cutoff, data.num_modes, num_bits), + SparseRowStore::kMaxSlots, + seen); + sweep(terms, + sampled_gens, + cutoff_function(CutoffType::Support, cutoff, data.num_modes, num_bits), + SparseRowStore::kMaxSlots, + seen); + } + } + // No claim that nothing fell back: a wide fixture term against a wide generator legitimately overflows + // a capacity of cutoff + |G|, and that case is covered on its own below. + BOOST_TEST(seen.sparse_terms > 100U); + BOOST_TEST(seen.cutoff_passed > 0U); + BOOST_TEST(seen.sparse_cutoff_failed > 0U); +} + +// A store row with no view: the kernel must take the dense path for that term alone and keep taking the +// sparse one for the rest. Rows are a mix, so both happen in the same gate. +BOOST_AUTO_TEST_CASE(term_product_falls_back_on_a_spilled_row) { + std::mt19937_64 rng(99U); + constexpr size_t kNumModes = 64; + std::vector terms; + for (size_t t = 0; t < 60; ++t) { + terms.push_back(test_utils::random_monomial(rng, kNumModes, 3)); + terms.push_back(paired_term(kNumModes, 9)); // 9 modes > the 4 slots below + } + std::vector gens{test_utils::random_monomial(rng, kNumModes, 2), + test_utils::random_monomial(rng, kNumModes, 4)}; + + Seen seen; + sweep(terms, gens, cutoff_function(CutoffType::Length, 6, kNumModes, 2 * kNumModes), 4, seen); + sweep(terms, gens, cutoff_function(CutoffType::Support, 6, kNumModes, 2 * kNumModes), 4, seen); + BOOST_TEST(seen.fallback_terms > 0U); + BOOST_TEST(seen.sparse_terms > 0U); +} + +// A product past the scratch capacity, with the store row itself perfectly representable: capacity is +// max_mode_bound() + the generator's locality, so a fully paired row well above the cutoff overflows it. +// sparse_toggle reports that rather than truncating, and the kernel must then answer densely. +BOOST_AUTO_TEST_CASE(term_product_falls_back_on_a_capacity_overflow) { + std::mt19937_64 rng(1010U); + constexpr size_t kNumModes = 64; + std::vector terms; + for (size_t t = 0; t < 40; ++t) { + terms.push_back(test_utils::random_monomial(rng, kNumModes, 3)); + // 12 modes: inside the 20-slot store rows below, past a capacity of 4 + |G|. + terms.push_back(paired_term(kNumModes, 12)); + } + std::vector gens{test_utils::random_monomial(rng, kNumModes, 2)}; + + Seen seen; + sweep(terms, gens, cutoff_function(CutoffType::Length, 4, kNumModes, 2 * kNumModes), 20, seen); + sweep(terms, gens, cutoff_function(CutoffType::Support, 4, kNumModes, 2 * kNumModes), 20, seen); + BOOST_TEST(seen.fallback_terms > 0U); + BOOST_TEST(seen.sparse_terms > 0U); + BOOST_TEST(seen.cutoff_passed > 0U); +} + +// The escape's own case, and the reason no capacity can settle the question: a fully paired product is kept +// unconditionally (xor_sum == 0), so a product both wider than the record and kept by the cutoff exists by +// construction. The generator is fully paired on modes the term already holds, so the product is the term +// minus those modes -- still fully paired, still far wider than a capacity of bound + |G|. +BOOST_AUTO_TEST_CASE(term_product_escapes_a_surviving_product_no_record_can_hold) { + constexpr size_t kNumModes = 64; + const std::vector terms{paired_term(kNumModes, 12)}; + Bitset gen(2 * kNumModes); + for (const size_t m : {3U, 4U}) { + gen.set(2 * m); + gen.set((2 * m) + 1); + } + const std::vector gens{gen}; + + Seen seen; + sweep(terms, gens, cutoff_function(CutoffType::Length, 4, kNumModes, 2 * kNumModes), 20, seen); + sweep(terms, gens, cutoff_function(CutoffType::Support, 4, kNumModes, 2 * kNumModes), 20, seen); + BOOST_TEST(seen.fallback_terms == 2U); + BOOST_TEST(seen.sparse_terms == 0U); + BOOST_TEST(seen.cutoff_passed == 2U); + BOOST_TEST(seen.escaped_records == 2U); + BOOST_TEST(seen.row_records == 0U); +} + +// A logical width narrower than the storage width, which is what storage_modes_for() produces for any +// mode count that is not a whole 32-mode block. The inactive modes are the low ones, so the codes cutoff +// drops them as a prefix of the ascending slots where the dense one shifts the whole register -- and the +// terms here deliberately occupy modes on both sides of that boundary. +BOOST_AUTO_TEST_CASE(term_product_sparse_kernel_matches_dense_in_a_narrow_active_window) { + std::mt19937_64 rng(7070U); + constexpr size_t kStorageModes = 64; + constexpr size_t kLogicalModes = 50; // an inactive prefix of 14 modes + std::vector terms; + for (size_t t = 0; t < 120; ++t) { + terms.push_back(test_utils::random_monomial(rng, kStorageModes, 6)); + } + std::vector gens{test_utils::random_monomial(rng, kStorageModes, 3), + test_utils::random_monomial(rng, kStorageModes, 2)}; + + Seen seen; + sweep(terms, + gens, + cutoff_function(CutoffType::Length, 4, kLogicalModes, 2 * kStorageModes), + SparseRowStore::kMaxSlots, + seen); + sweep(terms, + gens, + cutoff_function(CutoffType::Support, 4, kLogicalModes, 2 * kStorageModes), + SparseRowStore::kMaxSlots, + seen); + BOOST_TEST(seen.sparse_terms > 0U); + BOOST_TEST(seen.cutoff_passed > 0U); + BOOST_TEST(seen.sparse_cutoff_failed > 0U); +} + +// A cutoff that is neither concrete functor -- the basis-change form, which is a lambda on purpose -- has +// no codes counterpart, so the whole gate falls back. Asserted because the alternative to falling back is +// answering with the wrong cutoff, which no other test would catch. +BOOST_AUTO_TEST_CASE(term_product_falls_back_when_the_cutoff_has_no_codes_form) { + std::mt19937_64 rng(2020U); + constexpr size_t kNumModes = 32; + std::vector terms; + for (size_t t = 0; t < 80; ++t) { + terms.push_back(test_utils::random_monomial(rng, kNumModes, 5)); + } + std::vector gens{test_utils::random_monomial(rng, kNumModes, 3)}; + + // The identity basis, so the predicate is an ordinary length cutoff wrapped in a lambda: the answers + // must still match, and every term must have gone the dense way to produce them. + MonomialList basis; + for (size_t b = 0; b < 2 * kNumModes; ++b) { + Bitset single(2 * kNumModes); + single.set(b); + basis.push_back(single); + } + const auto wrapped = cutoff_function_basis_change(CutoffType::Length, 4, basis, kNumModes); + BOOST_REQUIRE(CutoffEvaluator{wrapped}.length_cutoff() == nullptr); + + Seen seen; + sweep(terms, gens, wrapped, SparseRowStore::kMaxSlots, seen); + BOOST_TEST(seen.sparse_terms == 0U); + BOOST_TEST(seen.fallback_terms == terms.size()); +} + +// A generator wider than one codes word cannot be encoded as a row at all, so the gate falls back +// wholesale. Exotic, but the branch exists and an unencodable generator must not be silently truncated +// into a *different* generator. +BOOST_AUTO_TEST_CASE(term_product_falls_back_on_a_generator_past_one_codes_word) { + std::mt19937_64 rng(3030U); + constexpr size_t kNumModes = 128; + std::vector terms; + for (size_t t = 0; t < 40; ++t) { + terms.push_back(test_utils::random_monomial(rng, kNumModes, 4)); + } + std::vector gens{paired_term(kNumModes, SparseRowStore::kMaxSlots + 5)}; + + Seen seen; + sweep(terms, gens, cutoff_function(CutoffType::Length, 6, kNumModes, 2 * kNumModes), 40, seen); + BOOST_TEST(seen.sparse_terms == 0U); + BOOST_TEST(seen.fallback_terms == terms.size()); +} + +// --------------------------------------------------------------------------------------------------- +// The width-bound dense kernel against the width-agnostic one. +// --------------------------------------------------------------------------------------------------- + +namespace { + +// Which of the two cutoff paths a narrow-kernel run took, per term. Both have to occur across the +// cases below or one of them ships compared against nothing. +struct NarrowSeen { + size_t terms = 0; + size_t word_cutoff_terms = 0; // passes() answered off the words + size_t evaluator_terms = 0; // passes() went through the CutoffEvaluator + size_t passed = 0; + size_t failed = 0; + size_t paired_rescues = 0; // kept although longer than the cutoff, i.e. the fully-paired clause +}; + +// One term through both dense kernels, comparing every answer the scan reads plus the product itself. +// Unlike the sparse comparison the records must agree *byte for byte*: both push a dense monomial, so +// any difference here is a difference in the product. +template +auto check_narrow_term(DenseTermProducts &reference_kernel, + DenseTermProductsW &candidate_kernel, + const OperatorIndex &packed, + size_t i, + size_t gen_pop, + unsigned int cutoff, + NarrowSeen &seen) -> void { + const auto reference = reference_kernel.product(packed, i); + const auto candidate = candidate_kernel.product(packed, i); + BOOST_TEST(candidate.overlap == reference.overlap); + BOOST_TEST(candidate.phase_factor == reference.phase_factor); + BOOST_TEST((candidate_kernel.product_row() == reference_kernel.product_row())); + BOOST_TEST(candidate_kernel.record_words() == reference_kernel.record_words()); + + const size_t new_pop = packed.popcount(i) + gen_pop - (2 * reference.overlap); + const bool reference_passes = reference_kernel.passes(new_pop); + BOOST_TEST(candidate_kernel.passes(new_pop) == reference_passes); + + if (reference_passes) { + for (const size_t rank_count : {2U, 3U, 8U}) { + BOOST_TEST(candidate_kernel.owner(rank_count) == reference_kernel.owner(rank_count)); + } + + VecZ unused_escapes; + VecZ reference_record = query_buffer(); + reference_kernel.push(QueryOut{reference_record, unused_escapes}, -1); + VecZ candidate_record = query_buffer(); + candidate_kernel.push(QueryOut{candidate_record, unused_escapes}, -1); + BOOST_TEST(candidate_record == reference_record, boost::test_tools::per_element()); + BOOST_TEST(unused_escapes.empty()); // neither dense kernel has anywhere to escape to + seen.paired_rescues += new_pop > cutoff ? 1 : 0; + } + + ++seen.terms; + seen.word_cutoff_terms += candidate_kernel.uses_word_cutoff() ? 1 : 0; + seen.evaluator_terms += candidate_kernel.uses_word_cutoff() ? 0 : 1; + seen.passed += reference_passes ? 1 : 0; + seen.failed += reference_passes ? 0 : 1; +} + +// Every term of `terms` against every generator of `gens`, at the one width W the terms are built for. +template +auto sweep_narrow(const std::vector &terms, + const std::vector &gens, + const CutoffFn &cutoff_fn, + unsigned int cutoff, + NarrowSeen &seen) -> void { + BOOST_REQUIRE(!terms.empty()); + const size_t num_bits = terms.front().size(); + BOOST_REQUIRE(num_bits == W * Bitset::word_width); + OperatorIndex packed(num_bits); + for (const auto &mono : terms) { + packed.push_back(mono); + } + + const CutoffEvaluator cutoff_eval{cutoff_fn}; + for (const auto &gen : gens) { + DenseTermProducts reference_kernel(gen, cutoff_eval); + DenseTermProductsW candidate_kernel(gen, cutoff_eval); + const size_t gen_pop = gen.count(); + for (size_t i = 0; i < terms.size(); ++i) { + check_narrow_term(reference_kernel, candidate_kernel, packed, i, gen_pop, cutoff, seen); + } + } +} + +// The W the dispatch bound, read back out of the arm it selected -- the only thing about the seam that +// is observable, since every arm answers identically. +template +auto bound_kernel_width(size_t num_words) -> size_t { + return with_kernel_width(num_words, [](std::integral_constant) -> size_t { return W; }); +} + +// Terms wide enough to be rejected, narrow enough to survive, and some fully paired well above the +// cutoff -- which is the clause the word cutoff answers with its own fold rather than a popcount. +auto narrow_kernel_terms(std::mt19937_64 &rng, size_t num_modes) -> std::vector { + std::vector terms; + for (size_t t = 0; t < 60; ++t) { + terms.push_back(test_utils::random_monomial(rng, num_modes, 7)); + } + for (const size_t paired : {1U, 3U, 8U}) { + terms.push_back(paired_term(num_modes, paired)); + } + return terms; +} + +} // namespace + +// A length cutoff over the whole register: the one shape the word cutoff answers, so this is the case +// where the specialized path is actually under test. Asserted through uses_word_cutoff(), because a +// change that quietly stopped engaging it would leave every other assertion here still passing. +BOOST_AUTO_TEST_CASE(narrow_kernel_matches_dense_under_a_whole_register_length_cutoff) { + std::mt19937_64 rng(20260821U); + NarrowSeen seen; + test_utils::for_each_inline_width([&](std::integral_constant) { + const size_t num_modes = (W * Bitset::word_width) / 2; + const auto terms = narrow_kernel_terms(rng, num_modes); + const std::vector gens{test_utils::random_monomial(rng, num_modes, 2), + test_utils::random_monomial(rng, num_modes, 4), + paired_term(num_modes, 2)}; + for (const unsigned int cutoff : {4U, 8U}) { + const auto length = cutoff_function(CutoffType::Length, cutoff, num_modes, 2 * num_modes); + sweep_narrow(terms, gens, length, cutoff, seen); + // The Pauli algebra with a length cutoff: the kind and the basis are independent here even + // though the shipping models pair them, and the kernel's cutoff arm must not depend on the + // algebra it was instantiated with. + sweep_narrow(terms, gens, length, cutoff, seen); + } + }); + BOOST_TEST(seen.terms > 0U); + BOOST_TEST(seen.evaluator_terms == 0U); // every term took the word cutoff + BOOST_TEST(seen.word_cutoff_terms == seen.terms); + BOOST_TEST(seen.passed > 0U); + BOOST_TEST(seen.failed > 0U); + // Products kept although longer than the cutoff, i.e. the fully-paired fold answering yes. Without + // these the word cutoff would be tested as nothing but `new_pop <= cutoff`. + BOOST_TEST(seen.paired_rescues > 0U); +} + +// The two shapes the word cutoff declines: a support cutoff, and a length cutoff whose active window is +// narrower than the storage register -- which is what storage_modes_for() produces for any mode count +// that is not a whole 32-mode block. Both must fall through to the evaluator and still agree; the +// failure this guards against is answering a *different* cutoff, which changes which terms survive. +BOOST_AUTO_TEST_CASE(narrow_kernel_defers_to_the_evaluator_off_the_whole_register_length_cutoff) { + std::mt19937_64 rng(606U); + NarrowSeen seen; + test_utils::for_each_inline_width([&](std::integral_constant) { + const size_t num_modes = (W * Bitset::word_width) / 2; + const size_t logical_modes = num_modes - 5; // an inactive prefix of 5 modes + const auto terms = narrow_kernel_terms(rng, num_modes); + const std::vector gens{test_utils::random_monomial(rng, num_modes, 3), + test_utils::random_monomial(rng, num_modes, 2)}; + for (const unsigned int cutoff : {4U, 6U}) { + const auto support = cutoff_function(CutoffType::Support, cutoff, num_modes, 2 * num_modes); + const auto narrow = cutoff_function(CutoffType::Length, cutoff, logical_modes, 2 * num_modes); + sweep_narrow(terms, gens, support, cutoff, seen); + sweep_narrow(terms, gens, support, cutoff, seen); + sweep_narrow(terms, gens, narrow, cutoff, seen); + } + }); + BOOST_TEST(seen.terms > 0U); + BOOST_TEST(seen.word_cutoff_terms == 0U); + BOOST_TEST(seen.evaluator_terms == seen.terms); + BOOST_TEST(seen.passed > 0U); + BOOST_TEST(seen.failed > 0U); +} + +// A cutoff that is neither concrete functor, so CutoffEvaluator recovered nothing: the kernel keeps its +// word product and defers the whole predicate. Same case as the sparse kernel's, for the same reason. +BOOST_AUTO_TEST_CASE(narrow_kernel_defers_when_the_cutoff_has_no_concrete_functor) { + std::mt19937_64 rng(707U); + constexpr size_t kWords = 2; + constexpr size_t kNumModes = (kWords * Bitset::word_width) / 2; + const auto terms = narrow_kernel_terms(rng, kNumModes); + const std::vector gens{test_utils::random_monomial(rng, kNumModes, 3)}; + + MonomialList basis; + for (size_t b = 0; b < 2 * kNumModes; ++b) { + Bitset single(2 * kNumModes); + single.set(b); + basis.push_back(single); + } + const auto wrapped = cutoff_function_basis_change(CutoffType::Length, 4, basis, kNumModes); + BOOST_REQUIRE(CutoffEvaluator{wrapped}.length_cutoff() == nullptr); + + NarrowSeen seen; + sweep_narrow(terms, gens, wrapped, 4, seen); + BOOST_TEST(seen.terms == terms.size()); + BOOST_TEST(seen.word_cutoff_terms == 0U); + BOOST_TEST(seen.passed > 0U); + BOOST_TEST(seen.failed > 0U); +} + +// The seam itself: which W the scan binds for a given storage width and store. Pinned because it is the +// one thing above that no differential test can see -- every arm computes the same answers, so a +// dispatch that always chose 0 would leave the whole suite green and only the benchmark different. +BOOST_AUTO_TEST_CASE(with_kernel_width_binds_the_capped_storage_word_count) { + for (size_t words = 1; words <= Bitset::kInlineWords; ++words) { + BOOST_TEST(bound_kernel_width(words) == (words <= kNarrowKernelWords ? words : 0)); + // The sparse store is never specialized: its per-term work is O(slots), not O(storage words). + BOOST_TEST(bound_kernel_width(words) == 0U); + } + // Above the inline regime the words are on the heap, so there is no width to bind. + BOOST_TEST(bound_kernel_width(Bitset::kInlineWords + 1) == 0U); +} + +// with_kernel_width binds W from the store's row word count, independently of whatever gen a caller +// hands to the constructor it dispatches into -- so a mismatch is a real (if never-yet-observed) call +// site bug, not an internal-only invariant, and unlike the per-term word ops in Bitset.h (deliberately +// assert-only) this must not compile away under NDEBUG. +BOOST_AUTO_TEST_CASE(narrow_kernel_constructor_rejects_a_generator_of_the_wrong_width) { + constexpr size_t kWords = 2; + constexpr size_t kNumModes = (kWords * Bitset::word_width) / 2; + const Bitset narrower_gen(2 * (kNumModes - Bitset::word_width / 2)); + const Bitset wider_gen(2 * (kNumModes + Bitset::word_width / 2)); + + const CutoffEvaluator cutoff_eval{cutoff_function(CutoffType::Length, 4, kNumModes, 2 * kNumModes)}; + BOOST_CHECK_THROW((DenseTermProductsW(narrower_gen, cutoff_eval)), KernelWidthMismatch); + BOOST_CHECK_THROW((DenseTermProductsW(wider_gen, cutoff_eval)), KernelWidthMismatch); +} diff --git a/cpp/tests/word_kernel_tests.cpp b/cpp/tests/word_kernel_tests.cpp new file mode 100644 index 00000000..a14a3a55 --- /dev/null +++ b/cpp/tests/word_kernel_tests.cpp @@ -0,0 +1,204 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// The scan's bound-width word passes against the computations they restate: detail::WordKernel's +// four Bitset methods, and detail::fully_paired_words, whose oracle is cutoff_sums rather than a +// Bitset method (which is why it lives in AlgebraCommon.h and is tested here anyway -- one sweep over +// the inline regime, one set of word patterns). Every one of them is a second copy of an existing +// computation, so the only thing worth testing is that the two copies agree -- at every W in the +// inline regime, on the word patterns that distinguish a per-word fold from a whole-register one. +// +// splitmix carries the strongest obligation and gets the strictest test: that value is monomial_hash, +// which routes MPI ownership, so a divergence would move terms between ranks rather than merely run +// slower. It is asserted equal to SplitmixHash for every W, not merely well-distributed. + +#include + +#include +#include +#include +#include +#include +#include + +#include "monoprop/Bitset.h" +// AlgebraCommon.h before Monomial.h and without Utilities.h: Utilities.h is not self-contained (it is +// reached through TypeAliases.h, which pulls MPOperator.h in ahead of the free functions that header +// calls), so even_bits arrives transitively here the same way term_product_tests.cpp gets it. +#include "monoprop/algebra/AlgebraCommon.h" +#include "monoprop/core/Monomial.h" + +#include "InlineWidths.h" + +using monoprop::Bitset; +using monoprop::detail::WordKernel; +using test_utils::for_each_inline_width; + +namespace { + +// A bitset of exactly W words with the given words written straight in. Going through data() rather +// than set() because these tests are about the words: a pattern like "every odd bit of word 3" is a +// word, not a bit list. +template +auto from_words(const std::array &words) -> Bitset { + Bitset bs(W * Bitset::word_width); + BOOST_REQUIRE(bs.num_words() == W); + for (size_t w = 0; w < W; ++w) { + bs.data()[w] = words[w]; + } + return bs; +} + +// The patterns a per-word fold can get wrong where a whole-register one cannot: all-zero and all-ones +// (parity of an even count of set bits), the two single-mode halves (a fold that dropped a word would +// still see one of them), and the paired/unpaired even-odd patterns fully_paired keys on. Randomized +// words follow in every case; these are the ones worth naming. +template +auto interesting_words(std::mt19937_64 &rng) -> std::vector> { + std::vector> out; + const auto fill = [&](uint64_t v) { + std::array a{}; + a.fill(v); + out.push_back(a); + }; + fill(0); + fill(~uint64_t{0}); + fill(0x5555555555555555ULL); // every even bit: every occupied mode singly occupied + fill(0xAAAAAAAAAAAAAAAAULL); // every odd bit: likewise, the other Majorana + fill(0xFFFFFFFFFFFFFFFFULL); // every mode fully paired + // One word at a time set, so a fold that skipped word w fails on exactly one entry. + for (size_t w = 0; w < W; ++w) { + std::array a{}; + a[w] = 0x0123456789ABCDEFULL; + out.push_back(a); + std::array b{}; + b[w] = uint64_t{1} << (w % Bitset::word_width); + out.push_back(b); + } + for (size_t t = 0; t < 64; ++t) { + std::array a{}; + for (auto &word : a) { + word = rng(); + } + out.push_back(a); + // ...and the same words sparsified, since real monomials are sparse and a dense random word + // never exercises the "one set bit in the whole register" shapes. + std::array sparse{}; + for (size_t k = 0; k < 3; ++k) { + const size_t pos = rng() % (W * Bitset::word_width); + sparse[pos / Bitset::word_width] |= uint64_t{1} << (pos % Bitset::word_width); + } + out.push_back(sparse); + } + return out; +} + +} // namespace + +// The owner-routing value. Equality with SplitmixHash (and hence monomial_hash) at every W, including +// the W == 1 arm both sides special-case. +BOOST_AUTO_TEST_CASE(word_kernel_splitmix_is_the_owner_routing_hash) { + std::mt19937_64 rng(20260821U); + for_each_inline_width([&](std::integral_constant) { + for (const auto &words : interesting_words(rng)) { + const Bitset bs = from_words(words); + const size_t expected = monoprop::SplitmixHash{}(bs); + BOOST_TEST(WordKernel::splitmix(bs.data()) == expected); + // The name that marks the value as pinned, asserted separately from the hash functor so a + // future indirection between them cannot pass this test by tautology. + BOOST_TEST(WordKernel::splitmix(bs.data()) == monoprop::monomial_hash(bs)); + } + }); +} + +// The product and its two counts, in the same destination the scan writes. +BOOST_AUTO_TEST_CASE(word_kernel_fused_xor_into_matches_bitset) { + std::mt19937_64 rng(31337U); + for_each_inline_width([&](std::integral_constant) { + const auto lhs_words = interesting_words(rng); + const auto rhs_words = interesting_words(rng); + for (size_t i = 0; i < lhs_words.size(); ++i) { + const Bitset lhs = from_words(lhs_words[i]); + const Bitset rhs = from_words(rhs_words[i]); + + Bitset reference_out(W * Bitset::word_width); + const auto reference = lhs.fused_xor_into(rhs, reference_out); + + Bitset candidate_out(W * Bitset::word_width); + const auto candidate = WordKernel::fused_xor_into(lhs.data(), rhs.data(), candidate_out.data()); + + BOOST_TEST(candidate.overlap == reference.overlap); + BOOST_TEST(candidate.result_count == reference.result_count); + BOOST_TEST((candidate_out == reference_out)); + } + }); +} + +// The Majorana rotation sign's parity, which is of the whole AND and not of any per-word rounding -- +// so the all-ones patterns above matter: they make the per-word popcounts even and the total even too, +// where a wrong fold would still agree. +BOOST_AUTO_TEST_CASE(word_kernel_parity_and_matches_bitset) { + std::mt19937_64 rng(4242U); + for_each_inline_width([&](std::integral_constant) { + const auto lhs_words = interesting_words(rng); + const auto rhs_words = interesting_words(rng); + for (size_t i = 0; i < lhs_words.size(); ++i) { + const Bitset lhs = from_words(lhs_words[i]); + const Bitset rhs = from_words(rhs_words[i]); + BOOST_TEST(WordKernel::parity_and(lhs.data(), rhs.data()) == lhs.parity_and(rhs)); + } + }); +} + +// The cutoff's fully-paired clause. The oracle is cutoff_sums' xor_sum over the whole register, which +// is the only window the pass is allowed to answer for (a narrower one keeps going through the +// evaluator -- see DenseTermProductsW). It carries the even-bit pattern as a literal, so the mask built +// here is also the check that the literal is what even_bits would have produced. +BOOST_AUTO_TEST_CASE(fully_paired_words_matches_cutoff_sums) { + std::mt19937_64 rng(5150U); + for_each_inline_width([&](std::integral_constant) { + const size_t num_bits = W * Bitset::word_width; + const Bitset mask = monoprop::even_bits(num_bits); + for (size_t w = 0; w < W; ++w) { + BOOST_TEST(mask.word(w) == 0x5555555555555555ULL); + } + size_t paired = 0; + size_t unpaired = 0; + for (const auto &words : interesting_words(rng)) { + const Bitset bs = from_words(words); + const bool expected = monoprop::cutoff_sums(bs, num_bits / 2).xor_sum == 0; + BOOST_TEST(monoprop::detail::fully_paired_words(bs.data()) == expected); + paired += expected ? 1 : 0; + unpaired += expected ? 0 : 1; + } + // Both answers occur, so neither is passing by always returning the same one. + BOOST_TEST(paired > 0U); + BOOST_TEST(unpaired > 0U); + }); +} + +// clear() zeroes exactly W words. The word above is checked because the kernel's whole contract is +// "the caller bound the width": writing one word too many would corrupt an unrelated monomial's +// storage, and no other test reads that word. +BOOST_AUTO_TEST_CASE(word_kernel_clear_zeroes_exactly_its_width) { + for_each_inline_width([&](std::integral_constant) { + std::array buffer{}; + buffer.fill(~uint64_t{0}); + WordKernel::clear(buffer.data()); + for (size_t w = 0; w < W; ++w) { + BOOST_TEST(buffer[w] == 0U); + } + BOOST_TEST(buffer[W] == ~uint64_t{0}); + }); +}