perf(cpp): ⚡ bind the per-term kernel and the support-form query record

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<Store, A, W>` 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<Store>` (`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<A, W>` and the `TermProductsFor<Store, A, W>` specialization
+  `DenseTermProductsW<A, W>` 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<NumBits>`
+  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<W>` (`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<W>`, 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<W>`, 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<Store>` (`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

cpp/monoprop/Bitset.h

@@ -60,9 +60,10 @@
     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
@@ -81,7 +82,7 @@
         overlap += static_cast<size_t>(std::popcount(a[i] & b[i]));
         result_count += static_cast<size_t>(std::popcount(n));
     }
     return {overlap, result_count};
 }
 
 // XOR-fold of a & b into one word. Folding first and popcounting once is what makes the caller's
@@ -118,7 +119,8 @@
 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;
@@ -227,7 +229,7 @@
         : nwords_(static_cast<uint32_t>(words_for(num_bits))),
           top_bits_(static_cast<uint32_t>(num_bits % word_width)) {
         if (spilled()) {
             s_.heap_ = new word_type[nwords_]{};
         }
         else {
             zero_inline();
@@ -245,7 +247,7 @@
 
     Bitset(const Bitset &o) noexcept : nwords_(o.nwords_), top_bits_(o.top_bits_) {
         if (spilled()) {
             s_.heap_ = new word_type[nwords_];
             std::memcpy(s_.heap_, o.s_.heap_, nwords_ * sizeof(word_type));
         }
         else {
@@ -289,7 +291,7 @@
         nwords_ = o.nwords_;
         top_bits_ = o.top_bits_;
         if (spilled()) {
             s_.heap_ = new word_type[nwords_];
             std::memcpy(s_.heap_, o.s_.heap_, nwords_ * sizeof(word_type));
         }
         else {
@@ -340,7 +342,7 @@
         if (nwords_ == 0) {
             return 0;
         }
         return top_bits_ != 0 ? (static_cast<size_t>(nwords_ - 1) * word_width + top_bits_)
                               : static_cast<size_t>(nwords_) * word_width;
     }
 
@@ -419,7 +421,7 @@
     // path wants this form: `out` is a scratch monomial that lives for the whole gate, so a term
     // costs the word pass alone -- no width recomputation, no spill test, no allocation, and none of
     // the copies a by-value result goes through on its way to the caller's variable.
     auto fused_xor_into(const Bitset &gen, Bitset &out) const noexcept -> FusedCounts;
 
     [[nodiscard]] auto fused_xor(const Bitset &gen) const noexcept -> FusedXor;
 
@@ -547,7 +549,7 @@
             }
             return static_cast<size_t>(-1);
         });
         return hit == static_cast<size_t>(-1) ? size() : hit;
     }
 
     [[nodiscard]] auto find_next(size_t pos) const noexcept -> size_t { // size() if none
@@ -597,7 +599,7 @@
 inline auto Bitset::fused_xor(const Bitset &gen) const noexcept -> FusedXor {
     Bitset result(size());
     const auto counts = fused_xor_into(gen, result);
     return {result, counts.overlap, counts.result_count};
 }
 
 // Bit-identical to the old per-width SplitmixHash<Bitset<NumBits>>: same mix(), same per-word fold
@@ -626,6 +628,68 @@
     }
 };
 
+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<NumBits> 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 <size_t W>
+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<size_t, W>{});
+    }
+
+    // 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<size_t, W>{})) & 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<size_t>(SplitmixHash::mix(a[0]));
+        }
+        else {
+            uint64_t h = 0;
+            for (size_t i = 0; i < W; ++i) {
+                h ^= SplitmixHash::mix(a[i] + static_cast<uint64_t>(i));
+            }
+            return static_cast<size_t>(h);
+        }
+    }
+};
+
+} // namespace detail
+
 } // namespace monoprop
 
 namespace std {

cpp/monoprop/algebra/Algebra.h

@@ -24,6 +24,7 @@
 #include <utility>
 
 #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 <size_t W>
+    static auto rotation_sign_words(const GenContext &ctx,
+                                    const Bitset::word_type *mono,
+                                    const Bitset::word_type * /*new_mono*/) -> int {
+        return detail::WordKernel<W>::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 <size_t W>
+    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;

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<LSb0> 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 <size_t W>
+[[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

cpp/monoprop/detail/evolution/layer_build/CMakeLists.txt

@@ -9,4 +9,5 @@ target_sources(
         "FusedApply.h"
         "Resolve.h"
         "Scan.h"
+        "TermProduct.h"
 )