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}); + }); +}