From 00b8db96457f56dc38b1f85ab4b22642af808f7f Mon Sep 17 00:00:00 2001 From: Aaron Miller Date: Thu, 20 Aug 2026 20:59:11 +0100 Subject: [PATCH 1/5] =?UTF-8?q?perf(evolution):=20=E2=9A=A1=20carry=20cros?= =?UTF-8?q?s-rank=20queries=20as=20a=20width-adaptive=20record?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Five changes to the layer-build query path, squashed because the middle three share one data structure and the first is not measurable without them. The structural cutoff is decided from a (k, d) digest carried inline rather than from the wider comparison it replaced. Cross-rank queries then travel as a width-adaptive record, dropping the fixed stride the older CompactQuery design assumed. Queries whose owner is this rank resolve from positions instead of going through the cross-rank path at all. Two further changes ride with them and are named here because a reviewer cannot revert along a mechanism the message does not mention. Sink::kStride becomes incoming_layout()/querier_layout(): ContractSink's stride was the fused width while on_response_block was called with the plain buffer and took query_phase's plain-stride default, so the two agreed only by accident, and a wrong layout reads a neighbouring record's phase -- a coefficient sign flip, not a crash. And bulk_insert gains a group-prefetched hashed path, shared by the cross-rank and self inserts. Tests cover the paths that can disagree: a dense reference implementation the sparse path is differentially compared against, the digest's tie cases, the spill boundary where a row exceeds the inline width, and the prefetched insert. OperatorIndex gains a public overflow_size() so a test can assert the dense and positional insert paths spill identically -- the claim this change rests on, and otherwise unreachable from outside the class. Co-Authored-By: Claude Opus 5 (1M context) --- cpp/monoprop/algebra/AlgebraCommon.h | 60 +++ cpp/monoprop/core/SparseMonomial.h | 35 ++ .../detail/evolution/layer_build/Common.h | 61 +-- .../detail/evolution/layer_build/Engine.h | 147 ++++-- .../detail/evolution/layer_build/QueryCodec.h | 138 +++++ .../detail/evolution/layer_build/Resolve.h | 114 +++-- .../detail/evolution/layer_build/Scan.h | 67 ++- .../evolution/layer_build/SparseQuery.h | 401 +++++++++++++++ cpp/monoprop/detail/operator/OperatorIndex.h | 164 +++++- cpp/tests/README.md | 24 +- cpp/tests/bulk_insert_tests.cpp | 188 +++++++ cpp/tests/dense_query_reference.h | 72 +++ cpp/tests/digest_cutoff_tests.cpp | 143 ++++++ cpp/tests/fused_query_codec_tests.cpp | 125 ----- cpp/tests/mpi_utils_tests.cpp | 105 +++- cpp/tests/sparse_monomial_tests.cpp | 159 ++++++ cpp/tests/sparse_query_tests.cpp | 482 ++++++++++++++++++ cpp/tests/sparse_resolve_tests.cpp | 370 ++++++++++++++ 18 files changed, 2569 insertions(+), 286 deletions(-) create mode 100644 cpp/monoprop/core/SparseMonomial.h create mode 100644 cpp/monoprop/detail/evolution/layer_build/QueryCodec.h create mode 100644 cpp/monoprop/detail/evolution/layer_build/SparseQuery.h create mode 100644 cpp/tests/bulk_insert_tests.cpp create mode 100644 cpp/tests/dense_query_reference.h create mode 100644 cpp/tests/digest_cutoff_tests.cpp delete mode 100644 cpp/tests/fused_query_codec_tests.cpp create mode 100644 cpp/tests/sparse_monomial_tests.cpp create mode 100644 cpp/tests/sparse_query_tests.cpp create mode 100644 cpp/tests/sparse_resolve_tests.cpp diff --git a/cpp/monoprop/algebra/AlgebraCommon.h b/cpp/monoprop/algebra/AlgebraCommon.h index 80ef32b9..87da0c16 100644 --- a/cpp/monoprop/algebra/AlgebraCommon.h +++ b/cpp/monoprop/algebra/AlgebraCommon.h @@ -16,6 +16,7 @@ #include #include +#include #include #include #include @@ -24,6 +25,7 @@ #include "monoprop/TypeAliases.h" #include "monoprop/Utilities.h" +#include "monoprop/core/SparseMonomial.h" #include "monoprop/detail/operator/RowAccess.h" namespace monoprop { @@ -96,6 +98,8 @@ auto is_paired(const VecZ &mono) -> bool { return is_paired(indices_to_bitset(mono)); } +// The (k, d) digest form of the predicate above is is_paired(size_t, size_t) in SparseMonomial.h. + template auto is_fully_paired(const VecZ &inds, const Rows &op) -> VecZ { VecZ result; @@ -155,6 +159,27 @@ template return {(first_pair ^ second_pair).count(), active_mono.count(), (first_pair | second_pair).count()}; } +// The same sums from a (k, d) digest; no logical_num_modes because the masking above is inert for a +// well-formed monomial (every set bit at physical position >= 2 * (NumModes - logical_num_modes)). +[[nodiscard]] inline constexpr auto cutoff_sums(size_t k, size_t d) noexcept -> CutoffSums { + return {k - (2 * d), k, k - d}; +} + +// d alone: mode m owns bits (2m, 2m+1) LSb0, so `w & (w >> 1)` masked to even bits counts each +// doubly-occupied mode once. The shift is word-local because a carry would land on odd bit 63. +// Same well-formedness precondition as cutoff_sums(k, d); a monomial built below the active offset by +// hand (majorana_cutoff_tests.cpp:79,101) must keep the bitset overload, which stays the oracle. +template +[[gnu::always_inline]] [[nodiscard]] inline auto paired_mode_count(const Monomial &mono) noexcept -> size_t { + constexpr auto even = even_bits<2 * NumModes, LSb0>(); + size_t d = 0; + for (size_t w = 0; w < Monomial::num_words(); ++w) { + const uint64_t word = mono.word(w); + d += static_cast(std::popcount(word & (word >> 1) & even.word(w))); + } + return d; +} + // Both cutoffs below keep a fully paired monomial (xor_sum == 0) unconditionally: those are the only // terms contributing to an expectation value against a product reference state, so bounding them by // length or support would discard signal. @@ -170,6 +195,11 @@ auto length_cutoff(const Monomial &mono, unsigned int cutoff) -> bool return length_cutoff(mono, cutoff, NumModes); } +// Digest form, on the same precondition as cutoff_sums(k, d). Width-independent, hence not templated. +[[nodiscard]] inline constexpr auto length_cutoff(size_t k, size_t d, unsigned int cutoff) noexcept -> bool { + return length_keeps(k, d, cutoff); +} + template auto support_cutoff(const Monomial &mono, unsigned int cutoff, size_t logical_num_modes) -> bool { const auto sums = cutoff_sums(mono, logical_num_modes); @@ -181,6 +211,11 @@ auto support_cutoff(const Monomial &mono, unsigned int cutoff) -> bool return support_cutoff(mono, cutoff, NumModes); } +// Digest form, on the same precondition as cutoff_sums(k, d). Width-independent, hence not templated. +[[nodiscard]] inline constexpr auto support_cutoff(size_t k, size_t d, unsigned int cutoff) noexcept -> bool { + return support_keeps(k, d, cutoff); +} + namespace detail { template @@ -243,6 +278,20 @@ class CutoffEvaluator { return cutoff_fn_(mono); } + // Same decision without cutoff_sums: the caller already knows k, so only d is computed. No + // popcount early-out, deliberately -- the digest is cheaper than the branch. nullopt if opaque. + auto passes_from_dense(const Monomial &mono, size_t k) const -> std::optional { + // paired_mode_count has no active_mask, so it agrees with cutoff_sums(mono, L) only above it. + assert(mono.find_first() >= active_bit_offset_() && "monomial has a set bit below its active offset"); + if (length_cutoff_ != nullptr) { + return length_keeps(k, paired_mode_count(mono), length_cutoff_->cutoff); + } + if (support_cutoff_ != nullptr) { + return support_keeps(k, paired_mode_count(mono), support_cutoff_->cutoff); + } + return std::nullopt; + } + // Upper bound on the set bits (physical slots) a surviving term can carry, so the store can size // its packed inline rows. A length cutoff counts set bits directly; a support cutoff counts // modes/qubits, each spanning two slots, hence the x2. @@ -257,6 +306,17 @@ class CutoffEvaluator { } private: + // 2 * (NumModes - logical_num_modes) of whichever concrete cutoff is configured; assert-only. + [[nodiscard]] auto active_bit_offset_() const -> size_t { + if (length_cutoff_ != nullptr) { + return 2 * (NumModes - length_cutoff_->logical_num_modes); + } + if (support_cutoff_ != nullptr) { + return 2 * (NumModes - support_cutoff_->logical_num_modes); + } + return 0; + } + const CutoffFn &cutoff_fn_; const LengthCutoff *length_cutoff_; const SupportCutoff *support_cutoff_; diff --git a/cpp/monoprop/core/SparseMonomial.h b/cpp/monoprop/core/SparseMonomial.h new file mode 100644 index 00000000..20dd5c08 --- /dev/null +++ b/cpp/monoprop/core/SparseMonomial.h @@ -0,0 +1,35 @@ +// 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 structural cutoffs over a monomial's (k, d) digest: k = popcount, d = modes carrying BOTH +// Majoranas. Lets CutoffEvaluator decide from integers the emit site already has, without cutoff_sums. + +#include + +namespace monoprop { + +// xor_sum = k - 2d, popcount_sum = k, or_sum = k - d; a fully paired monomial is kept unconditionally. +[[nodiscard]] inline constexpr auto is_paired(size_t k, size_t d) noexcept -> bool { + return k == 2 * d; +} +[[nodiscard]] inline constexpr auto length_keeps(size_t k, size_t d, size_t cutoff) noexcept -> bool { + return k == 2 * d || k <= cutoff; +} +[[nodiscard]] inline constexpr auto support_keeps(size_t k, size_t d, size_t cutoff) noexcept -> bool { + return k == 2 * d || k - d <= cutoff; +} + +} // namespace monoprop diff --git a/cpp/monoprop/detail/evolution/layer_build/Common.h b/cpp/monoprop/detail/evolution/layer_build/Common.h index e3c203d3..8777ce24 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Common.h +++ b/cpp/monoprop/detail/evolution/layer_build/Common.h @@ -104,24 +104,8 @@ struct FusedContract { std::vector cross_half; // R>1: one half per cross-rank query (resolver +φ, querier −φ) }; -// Queries ride flat VecZ buffers: kQueryWords elements per query (W monomial words + one ±1 phase word). -// The source index is not in the payload — the resolver answers by position; the querier holds src_idx_r[r][q]. -template -inline constexpr size_t kQueryWords = mpi_detail::kWords + 1; - -// 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. -template -inline constexpr size_t kQueryWordsFused = kQueryWords + 1; - -// The unsigned-int intermediate normalizes the ±1 sign bit into a fixed 32-bit pattern so the round-trip -// is exact for any VecZ element width. Edit encode/decode as a pair. -inline auto encode_phase(int phase) -> size_t { - return static_cast(static_cast(phase)); -} -inline auto decode_phase(size_t word) -> int { - return static_cast(static_cast(word)); -} +// Queries ride flat VecZ buffers in one VARIABLE-WIDTH format (SparseQuery): no stride exists, so every +// offset comes from QueryCodec's walk. The source index is not on the wire; the querier holds src_idx_r. // bit_cast, not a conversion, so v_src arrives over the wire bit-identical. static_assert(sizeof(size_t) == sizeof(double), "fused query value word assumes 64-bit VecZ element"); @@ -132,45 +116,4 @@ inline auto decode_value(size_t word) -> double { return std::bit_cast(word); } -template -inline auto query_push(VecZ &buf, const Monomial &mono, int phase) -> void { - mpi_detail::append_monomial_words(mono, buf); - buf.push_back(encode_phase(phase)); -} - -// The mono + phase words occupy the same leading offsets in the plain and fused record, so readers differ -// only in the per-record stride QW (defaulted to the plain width). -template > -inline auto query_read(const VecZ &buf, size_t q, Monomial &mono_out, int &phase_out) -> void { - const size_t base = q * QW; - mono_out = mpi_detail::read_monomial_from_words(buf, base); - phase_out = decode_phase(buf[base + mpi_detail::kWords]); -} - -// No monomial reconstruction: process_responses needs only the phase. -template > -inline auto query_phase(const VecZ &buf, size_t q) -> int { - return decode_phase(buf[q * QW + mpi_detail::kWords]); -} - -template -inline auto query_value(const VecZ &buf, size_t q) -> double { - return decode_value(buf[q * kQueryWordsFused + mpi_detail::kWords + 1]); -} - -// Requires v.size() == q.size()/kQueryWords: exactly one value per query record. -template -inline auto build_fused_query_value(const VecZ &q, const std::vector &v, VecZ &out) -> void { - constexpr size_t W = kQueryWords; - const size_t nq = q.empty() ? 0 : q.size() / W; - out.clear(); - out.reserve(nq * kQueryWordsFused); - for (size_t i = 0; i < nq; ++i) { - out.insert(out.end(), - q.begin() + static_cast(i * W), - q.begin() + static_cast((i + 1) * W)); - out.push_back(encode_value(v[i])); - } -} - } // namespace monoprop::detail diff --git a/cpp/monoprop/detail/evolution/layer_build/Engine.h b/cpp/monoprop/detail/evolution/layer_build/Engine.h index d073f00c..88983362 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Engine.h +++ b/cpp/monoprop/detail/evolution/layer_build/Engine.h @@ -29,6 +29,7 @@ #include "monoprop/algebra/Algebra.h" #include "monoprop/detail/evolution/CutoffContext.h" #include "monoprop/detail/evolution/layer_build/Common.h" +#include "monoprop/detail/evolution/layer_build/QueryCodec.h" #include "monoprop/detail/evolution/layer_build/Resolve.h" #include "monoprop/detail/evolution/layer_build/Scan.h" #include "monoprop/detail/graph_encoding/MPGraphEncodingStorage.h" @@ -70,7 +71,10 @@ inline auto append_inserted_endpoints(CosMask &cos_all, size_t combined_size, co template struct GraphSink { static constexpr bool wants_values = false; - static constexpr size_t kStride = kQueryWords; + // Named apart: incoming_layout is what this rank RECEIVES, querier_layout its OWN send buffer. They + // coincide here only because GraphSink never fuses -- see ContractSink::querier_layout. + [[nodiscard]] auto incoming_layout() const -> QueryLayout { return {/*fused=*/false}; } + [[nodiscard]] auto querier_layout() const -> QueryLayout { return {/*fused=*/false}; } using Response = TermIndex; static auto init_response() -> Response { return std::numeric_limits::max(); } @@ -134,11 +138,16 @@ struct GraphSink { auto &out = acc[r].out_entries; const size_t base = out.size(); const size_t nq = resp.size(); + const QueryLayout layout = querier_layout(); out.resize(base + nq); + // Forward walk, not indexing by q: a compact query's width depends on its own popcount. + size_t off = 0; for (size_t q = 0; q < nq; ++q) { assert(resp[q] != std::numeric_limits::max() && "resolver must insert absent cross-rank terms"); - out[base + q] = {srcs[q], query_phase(qbuf, q)}; + out[base + q] = {srcs[q], QueryCodec::phase_at(qbuf, off)}; + off = QueryCodec::next_off(qbuf, layout, off); } + assert(off == qbuf.size() && "querier buffer does not hold exactly one query per response"); } // Drains the per-rank accumulators into the LayerCore's sin_send/sin_recv lists (layout derivation: @@ -186,7 +195,11 @@ struct GraphSink { template struct ContractSink { static constexpr bool wants_values = true; - static constexpr size_t kStride = kQueryWordsFused; + // This rank RECEIVES fused (query+value) records, but the buffer on_response_block is handed is its + // own queries_r, which is PLAIN (build_fused writes the fused form into combined_qv_). Reading the + // phase with the wrong layout takes a neighbouring record's, which is a silent coefficient sign flip. + [[nodiscard]] auto incoming_layout() const -> QueryLayout { return {/*fused=*/true}; } + [[nodiscard]] auto querier_layout() const -> QueryLayout { return {/*fused=*/false}; } using Response = double; static auto init_response() -> Response { return 0.0; } @@ -226,7 +239,7 @@ struct ContractSink { -> std::vector & { scratch.resize(queries.size()); for (size_t r = 0; r < queries.size(); ++r) { - build_fused_query_value(queries[r], vals[r], scratch[r]); + QueryCodec::build_fused(queries[r], vals[r], scratch[r]); } return scratch; } @@ -240,7 +253,7 @@ struct ContractSink { } auto on_resolved(size_t g, size_t s, - size_t q, + size_t /*q*/, size_t ip, const IncomingProbe &pr, const std::vector &incoming) -> Response { @@ -249,16 +262,19 @@ struct ContractSink { v_tgt = fused_scale ? op_coeffs[ip] * inv_cos : op_coeffs[ip]; } else if (schrodinger) { - v_tgt = - is_paired(pr.mono[g]) ? algebra_state_phase(basis, pr.mono[g], state_mask_) : 0.0; + // Through the probe's accessors: it holds position lists, and mono_at builds a bitset only + // for the fully paired minority that is_paired_at admits. + v_tgt = pr.is_paired_at(g) ? algebra_state_phase(basis, pr.mono_at(g), state_mask_) : 0.0; } else { v_tgt = 0.0; // Heisenberg fresh insert } - fc.cross_half[cross_base_ + g] = HalfRotationRec{ip, - query_value(incoming[s], q), - static_cast(pr.phase_of[g]), - /*is_insert=*/ip >= pr.base}; + // pr.off_of[g], not q: under the compact record a query ordinal does not name a buffer position. + fc.cross_half[cross_base_ + g] = + HalfRotationRec{ip, + QueryCodec::value_at(incoming[s], incoming_layout(), pr.off_of[g]), + static_cast(pr.phase_of[g]), + /*is_insert=*/ip >= pr.base}; return v_tgt; } auto process_reserve(const std::vector> &inc_r, size_t rank_count, size_t my_rank_) -> void { @@ -276,10 +292,14 @@ struct ContractSink { const std::vector &srcs, const VecZ &qbuf) -> void { const size_t nq = rval.size(); + const QueryLayout layout = querier_layout(); + size_t off = 0; for (size_t q = 0; q < nq; ++q) { - const auto nphase = static_cast(-query_phase(qbuf, q)); + const auto nphase = static_cast(-QueryCodec::phase_at(qbuf, off)); fc.cross_half.push_back(HalfRotationRec{srcs[q], rval[q], nphase, /*is_insert=*/false}); + off = QueryCodec::next_off(qbuf, layout, off); } + assert(off == qbuf.size() && "querier buffer does not hold exactly one query per response"); } // No LayerCore in the fused path → nullptr. Two-pass fused (k>0 / cos==0 fallback) appends inserted @@ -297,8 +317,14 @@ struct ContractSink { // Owns build_layer's machinery over a compile-time Sink policy. combined_size = the pre-layer operator size. template struct LayerBuildEngine { + // The store's position type, narrower than the wire's below 129 modes; decoded straight into. + using RowPosT = typename OperatorIndex::PosT; + + // A miss keeps its decoded positions (pos_at indexes deferred_pos_flat_) and the probe's hash. struct DeferredSelfMiss { - Monomial mono; + size_t pos_at; + uint32_t k; + uint32_t hash; size_t src; int phase; double v_src = 0.0; // ContractSink only: op_pre[src] captured at scan emit; 0 for GraphSink @@ -314,6 +340,10 @@ struct LayerBuildEngine { std::vector queries_r; std::vector> src_idx_r; std::vector deferred_self_misses; + // Deferred-miss positions, concatenated in miss order; parallel to deferred_self_misses. + std::vector deferred_pos_flat_; + // Per-batch decode scratch for resolve_range_; a member so one allocation serves every batch. + std::vector self_pos_flat_; // Scan-captured v_src per query (ContractSink only via Sink::wants_values; empty for GraphSink). std::vector> src_val_r; // Fused query+value send scratch (ContractSink, R>1): shared by a gate's two exchange passes. @@ -347,8 +377,10 @@ struct LayerBuildEngine { if constexpr (Sink::wants_values) { lv = &src_val_r[my_rank]; } - const size_t nq = lq.empty() ? 0 : lq.size() / kQueryWords; - resolve_range_(lq, ls, lv, 0, nq, is_leader_pass); + // One source per query, pushed by the scan, so the count needs neither a walk nor a division. + assert(ls.size() == QueryCodec::count_queries(lq, sink.querier_layout()) + && "the self query buffer does not hold exactly one query per source"); + resolve_range_(lq, ls, lv, ls.size(), is_leader_pass); lq.clear(); ls.clear(); if constexpr (Sink::wants_values) { @@ -389,7 +421,8 @@ struct LayerBuildEngine { // Followers a leader already matched must not be re-resolved over the wire, so compact them out. auto drop_matched_cross_rank_followers() -> void { - constexpr size_t W = kQueryWords; + using QC = QueryCodec; + const QueryLayout layout = sink.querier_layout(); for (size_t r = 0; r < R; ++r) { if (r == my_rank) { continue; @@ -403,22 +436,23 @@ struct LayerBuildEngine { } const size_t nq = s.size(); size_t kept = 0; + // Two cursors, since a dropped query has no fixed width; order is the accumulation order. + size_t src_off = 0; + size_t dst_off = 0; for (size_t k = 0; k < nq; ++k) { - if (matched.is_marked(s[k])) { - continue; - } - if (kept != k) { - std::copy(q.begin() + static_cast(k * W), - q.begin() + static_cast((k + 1) * W), - q.begin() + static_cast(kept * W)); - } - s[kept] = s[k]; - if (v != nullptr) { - (*v)[kept] = (*v)[k]; + const size_t next = QC::next_off(q, layout, src_off); + if (!matched.is_marked(s[k])) { + dst_off += QC::move_query(q, layout, src_off, dst_off); + s[kept] = s[k]; + if (v != nullptr) { + (*v)[kept] = (*v)[k]; + } + ++kept; } - ++kept; + src_off = next; } - q.resize(kept * W); + assert(src_off == q.size() && "follower compaction did not consume the whole query buffer"); + q.resize(dst_off); s.resize(kept); if (v != nullptr) { v->resize(kept); @@ -435,13 +469,19 @@ struct LayerBuildEngine { if (n_miss == 0) { return; } - auto key_at = [&](size_t k) -> const Monomial & { return deferred_self_misses[k].mono; }; sink.prepare_deferred(n_miss); - insert_absent_terms(local_op, n_miss, key_at, [&](size_t k, size_t base) { + // insert_absent_terms' three steps without its dense round trips, on the same ordering contract: + // miss k lands at base+k, in leader-then-follower order. + // insert_absent_terms is the dense reference this path is differentially tested against + // (sparse_resolve_tests.cpp), so it must not be deleted for having no library caller. + const size_t base = local_op.store->grow_rows_geometric(n_miss); + for (size_t k = 0; k < n_miss; ++k) { const auto &m = deferred_self_misses[k]; - assign_row(*local_op.store, base + k, m.mono); + local_op.store->set_positions(base + k, deferred_pos_flat_.data() + m.pos_at, m.k); sink.emit_deferred(k, base + k, m.src, m.phase, m.v_src); - }); + } + local_op.store->bulk_insert_hashed(n_miss, base, [&](size_t j) { return deferred_self_misses[j].hash; }); + local_op.reindex_after_growth(base, n_miss); } auto finish(CosMask &&cos_all, CosMask *out_cos = nullptr) -> std::shared_ptr { @@ -455,7 +495,10 @@ struct LayerBuildEngine { auto response_recv_counts() const -> std::vector { std::vector counts(R); for (size_t r = 0; r < R; ++r) { - counts[r] = static_cast(queries_r[r].size() / kQueryWords); + // One response per QUERY, and src_idx_r[r] holds one source per query: no walk, no division. + assert(src_idx_r[r].size() == QueryCodec::count_queries(queries_r[r], sink.querier_layout()) + && "a querier buffer does not hold exactly one query per source"); + counts[r] = static_cast(src_idx_r[r].size()); } return counts; } @@ -466,24 +509,38 @@ struct LayerBuildEngine { auto resolve_range_(VecZ &lq, std::vector &ls, [[maybe_unused]] std::vector *lv, - size_t lo, size_t hi, bool is_leader_pass) -> void { const size_t op_size = local_op.store->size(); - std::array, kResolveBatch> keys; + // pos_off/k_of index self_pos_flat_, rebuilt per batch but keeping its capacity; no dense keys. + std::array pos_off; + std::array k_of; + std::array hashes; std::array phases; std::array srcs; std::array vals; std::array found; - size_t q = lo; + using QC = QueryCodec; + const QueryLayout layout = sink.querier_layout(); + // A cursor, not an ordinal, and it must advance even when the query is skipped. + size_t off = 0; + size_t q = 0; while (q < hi) { size_t m = 0; + self_pos_flat_.clear(); for (; q < hi && m < kResolveBatch; ++q) { const size_t src = ls[q]; + const size_t this_off = off; if (!is_leader_pass && matched.is_marked(src)) { + off = QC::next_off(lq, layout, this_off); continue; // follower already matched by a leader → not an independent rotation } - query_read(lq, q, keys[m], phases[m]); + const size_t k = QC::k_at(lq, this_off); + const size_t at = self_pos_flat_.size(); + self_pos_flat_.resize(at + k); // default-init grow: read_positions writes every element + off = QC::read_positions(lq, layout, this_off, self_pos_flat_.data() + at, phases[m]); + pos_off[m] = at; + k_of[m] = static_cast(k); srcs[m] = src; if constexpr (Sink::wants_values) { vals[m] = (*lv)[q]; @@ -493,7 +550,13 @@ struct LayerBuildEngine { if (m == 0) { break; } - local_op.store->find_batch(keys.data(), m, found.data()); + // The hashes come back because a miss needs one at insert, folded from these same positions. + local_op.store->find_batch_positions(self_pos_flat_.data(), + pos_off.data(), + k_of.data(), + m, + found.data(), + hashes.data()); for (size_t j = 0; j < m; ++j) { double v_src = 0.0; if constexpr (Sink::wants_values) { @@ -508,7 +571,11 @@ struct LayerBuildEngine { sink.self_hit(srcs[j], found[j], phases[j], v_src); } else { - deferred_self_misses.push_back({keys[j], srcs[j], phases[j], v_src}); + // self_pos_flat_ is cleared by the next batch, so copy now, into one gate-long buffer. + const size_t at = deferred_pos_flat_.size(); + const auto *const first = self_pos_flat_.data() + pos_off[j]; + deferred_pos_flat_.insert(deferred_pos_flat_.end(), first, first + k_of[j]); + deferred_self_misses.push_back({at, k_of[j], hashes[j], srcs[j], phases[j], v_src}); } } } diff --git a/cpp/monoprop/detail/evolution/layer_build/QueryCodec.h b/cpp/monoprop/detail/evolution/layer_build/QueryCodec.h new file mode 100644 index 00000000..a7439a30 --- /dev/null +++ b/cpp/monoprop/detail/evolution/layer_build/QueryCodec.h @@ -0,0 +1,138 @@ +// 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 + +#include "monoprop/detail/evolution/layer_build/Common.h" +#include "monoprop/detail/evolution/layer_build/SparseQuery.h" + +namespace monoprop::detail { + +// The one interface every site that walks a query buffer is written against: the record is variable +// width, so no caller may hold a stride. + +// `fused` is a property of the BUFFER, not the process: queries_r is always plain, while the send +// scratch and what a ContractSink resolver receives are fused. A named field, not a bare bool, so +// `next_off(buf, true, off)` cannot read as plausibly-correct-either-way. +struct QueryLayout { + bool fused = false; // one value word follows each query +}; + +// One alias, so the tests and the codec name the same record type. +template +using QueryRecord = SparseQuery; + +template +struct QueryCodec { + using CQ = QueryRecord; + using PosT = typename CQ::PosT; + + // Words the QUERY at `off` occupies, NOT counting a trailing fused value word. Asked of the record + // rather than derived from k, which does not determine the width once the mode and gw can vary. + [[nodiscard]] static auto query_words(const VecZ &buf, size_t off) -> size_t { return CQ::words_at(buf, off); } + + // Complete mode pairs among ascending positions, exposed so no caller names a concrete record type. + template + [[nodiscard]] static auto pair_count(const OutT *pos, size_t k) noexcept -> size_t { + return CQ::pair_count(pos, k); + } + + // Reserve hints, not correctness: sized from the measured mean of 5.33 positions per query. + static constexpr size_t kReservePositionsPerQuery = 6; + static constexpr size_t kReserveWordsPerQuery = 2; + + // Offset of the next query; `off` always names the START of one, and the rest is derived. + [[nodiscard]] static auto next_off(const VecZ &buf, QueryLayout layout, size_t off) -> size_t { + return off + query_words(buf, off) + (layout.fused ? 1U : 0U); + } + + // Returns the WORDS written, which is not a constant: byte accounting must not assume a width. + static auto push(VecZ &buf, const Monomial &mono, int phase) -> size_t { + return CQ::push_mono(buf, mono, phase); + } + + // Identical in both formats: the value is one bit_cast word after the query's words. + static auto push_value(VecZ &buf, double v) -> void { buf.push_back(encode_value(v)); } + + // Inflates the record back into a dense Monomial, for callers that cannot consume positions. + static auto read_mono(const VecZ &buf, size_t off, Monomial &mono_out, int &phase_out) -> void { + (void)CQ::read_mono(buf, off, mono_out, phase_out); + } + + // The query's popcount, straight out of the record's header field. + [[nodiscard]] static auto k_at(const VecZ &buf, size_t off) -> size_t { return CQ::k_at(buf, off); } + + // Positions plus phase, into the CALLER's element type: the store's PosT is narrower below 129 modes. + template + static auto read_positions(const VecZ &buf, QueryLayout layout, size_t off, OutT *out, int &phase_out) -> size_t { + phase_out = CQ::phase_at(buf, off); + return CQ::read_positions(buf, off, out) + (layout.fused ? 1U : 0U); + } + + [[nodiscard]] static auto phase_at(const VecZ &buf, size_t off) -> int { return CQ::phase_at(buf, off); } + + [[nodiscard]] static auto value_at(const VecZ &buf, [[maybe_unused]] QueryLayout layout, size_t off) -> double { + assert(layout.fused && "there is no value word in a plain query buffer"); + return decode_value(buf[off + query_words(buf, off)]); + } + + // The number of QUERIES. Genuinely a walk: records vary in width, so there is no stride to divide by. + [[nodiscard]] static auto count_queries(const VecZ &buf, QueryLayout layout) -> size_t { + size_t off = 0; + size_t n = 0; + while (off < buf.size()) { + off = next_off(buf, layout, off); + ++n; + } + assert(off == buf.size() && "a compact query ran past the end of the buffer"); + return n; + } + + // Interleave a plain query stream with its parallel v_src array; a size mismatch shifts every coeff. + static auto build_fused(const VecZ &queries, const std::vector &vals, VecZ &out) -> void { + out.clear(); + out.reserve(queries.size() + vals.size()); + size_t off = 0; + size_t i = 0; + while (off < queries.size()) { + const size_t n = query_words(queries, off); + out.insert(out.end(), + queries.begin() + static_cast(off), + queries.begin() + static_cast(off + n)); + assert(i < vals.size() && "fused build needs exactly one value per query"); + out.push_back(encode_value(vals[i])); + off += n; + ++i; + } + assert(i == vals.size() && "fused build needs exactly one value per query"); + } + + // Copy the query at `src_off` (with its value word, if fused) to `dst_off`; returns words written. + static auto move_query(VecZ &buf, QueryLayout layout, size_t src_off, size_t dst_off) -> size_t { + const size_t n = query_words(buf, src_off) + (layout.fused ? 1U : 0U); + if (src_off != dst_off) { + assert(dst_off < src_off && "compaction only ever moves a query earlier"); + std::copy(buf.begin() + static_cast(src_off), + buf.begin() + static_cast(src_off + n), + buf.begin() + static_cast(dst_off)); + } + return n; + } +}; + +} // namespace monoprop::detail diff --git a/cpp/monoprop/detail/evolution/layer_build/Resolve.h b/cpp/monoprop/detail/evolution/layer_build/Resolve.h index ef90d91f..0c49706f 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Resolve.h +++ b/cpp/monoprop/detail/evolution/layer_build/Resolve.h @@ -14,6 +14,7 @@ #pragma once +#include #include #include #include @@ -22,6 +23,7 @@ #include "monoprop/algebra/Algebra.h" #include "monoprop/detail/evolution/CutoffContext.h" #include "monoprop/detail/evolution/layer_build/Common.h" +#include "monoprop/detail/evolution/layer_build/QueryCodec.h" #include "monoprop/detail/operator/MPOperator.h" #include "monoprop/detail/operator/RowAccess.h" @@ -33,28 +35,58 @@ namespace monoprop::detail { // pairwise distinct ⇒ misses distinct and absent. template struct IncomingProbe { - std::vector goff; // rank_count+1 flat offsets: g = goff[s] + q - DefaultInitVector sender_of; // g → sender rank - DefaultInitVector> mono; // g → deserialized query monomial - DefaultInitVector phase_of; // g → query phase - DefaultInitVector idx_of; // g → resolved index (hit: < base; miss: base+j) - std::vector miss_g; // j → the g that became miss j (Phase 4 reads mono[miss_g[j]]) - size_t base = 0; // op size before the miss inserts (the miss-index base) + // The STORE's position width, not the wire's: these positions exist to become rows. + using PosT = typename OperatorIndex::PosT; + + std::vector goff; // rank_count+1 flat offsets: g = goff[s] + q + DefaultInitVector sender_of; // g → sender rank + DefaultInitVector phase_of; // g → query phase + // g → WORD offset of that query inside incoming[sender_of[g]]; a query ordinal names no position. + DefaultInitVector off_of; + DefaultInitVector idx_of; // g → resolved index (hit: < base; miss: base+j) + std::vector miss_g; // j → the g that became miss j (Phase 4 reads the key of miss_g[j]) + size_t base = 0; // op size before the miss inserts (the miss-index base) size_t nq_total = 0; + + // The queries as they arrived, flat: query g owns pos_flat[pos_off[g] .. pos_off[g] + k_of[g]). + DefaultInitVector pos_flat; + DefaultInitVector pos_off; + DefaultInitVector k_of; + // g → fold_hash of the query key, folded by the probe and reused by the insert. + DefaultInitVector hash_of; + + // BUILDS a bitset, so cold consumers only -- the fully paired minority, never anything per-term. + [[nodiscard]] auto mono_at(size_t g) const -> Monomial { + Monomial m; + const PosT *p = pos_flat.data() + pos_off[g]; + for (size_t j = 0; j < k_of[g]; ++j) { + m.set(static_cast(p[j])); + } + return m; + } + + // is_paired from the positions' (k, d) digest, no bitset built. + [[nodiscard]] auto is_paired_at(size_t g) const -> bool { + const PosT *p = pos_flat.data() + pos_off[g]; + const size_t k = k_of[g]; + return monoprop::is_paired(k, QueryCodec::pair_count(p, k)); + } }; -// Phases 1-2, read-only w.r.t. operator contents. QW = per-record stride: the plain query width, or -// kQueryWordsFused for the fused resolver. The caller runs Phase 3, then insert_incoming_misses. -template > +// Phases 1-2, read-only w.r.t. operator contents. `layout` describes the records this rank RECEIVES: +// fused for the ContractSink resolver, plain for GraphSink. The caller runs Phase 3, then +// insert_incoming_misses. Counts and offsets come from the decode walk; there is no record stride. +template auto probe_incoming_queries(const std::vector &incoming, // serialized, one VecZ per sender MPOperator &op, - size_t rank_count) -> IncomingProbe { - constexpr size_t W = QW; + size_t rank_count, + QueryLayout layout) -> IncomingProbe { + using QC = QueryCodec; IncomingProbe pr; pr.goff.assign(rank_count + 1, 0); for (size_t s = 0; s < rank_count; ++s) { - const size_t nq = incoming[s].empty() ? 0 : incoming[s].size() / W; + const size_t nq = QC::count_queries(incoming[s], layout); pr.goff[s + 1] = pr.goff[s] + nq; } pr.nq_total = pr.goff[rank_count]; @@ -69,22 +101,39 @@ auto probe_incoming_queries(const std::vector &incoming, // serialized, on static_cast(s)); } - // Phase 1 (read-only): deserialize, then probe with the group-prefetch batch find. - pr.mono.resize(pr.nq_total); + // Phase 1 (read-only): deserialize, then probe with the group-prefetch batch find. One walk per sender. pr.phase_of.resize(pr.nq_total); + pr.off_of.resize(pr.nq_total); pr.idx_of.resize(pr.nq_total); - for (size_t g = 0; g < pr.nq_total; ++g) { - const size_t s = pr.sender_of[g]; - const size_t q = g - pr.goff[s]; - Monomial m; - int ph = 0; - query_read(incoming[s], q, m, ph); - pr.mono[g] = m; - pr.phase_of[g] = ph; + pr.pos_off.resize(pr.nq_total); + pr.k_of.resize(pr.nq_total); + pr.hash_of.resize(pr.nq_total); + pr.pos_flat.clear(); + // A hint only: the measured mean is 5.33 positions, so this is one allocation but for an outlier. + pr.pos_flat.reserve(pr.nq_total * QueryCodec::kReservePositionsPerQuery); + for (size_t s = 0; s < rank_count; ++s) { + size_t off = 0; + for (size_t g = pr.goff[s]; g < pr.goff[s + 1]; ++g) { + int ph = 0; + const size_t k = QC::k_at(incoming[s], off); + const size_t at = pr.pos_flat.size(); + pr.pos_flat.resize(at + k); // default-init grow: read_positions writes every element + pr.pos_off[g] = at; + pr.k_of[g] = static_cast(k); + pr.off_of[g] = off; + off = QC::read_positions(incoming[s], layout, off, pr.pos_flat.data() + at, ph); + pr.phase_of[g] = ph; + } + assert(off == incoming[s].size() && "the query walk did not consume the sender's whole buffer"); } { const size_t op_size = op.store->size(); - op.store->find_batch(pr.mono.data(), pr.nq_total, pr.idx_of.data()); + op.store->find_batch_positions(pr.pos_flat.data(), + pr.pos_off.data(), + pr.k_of.data(), + pr.nq_total, + pr.idx_of.data(), + pr.hash_of.data()); for (size_t g = 0; g < pr.nq_total; ++g) { if (pr.idx_of[g] >= op_size) { // kNotFound is size_t max → also lands here pr.idx_of[g] = kMissingIndex; @@ -111,11 +160,15 @@ auto insert_incoming_misses(MPOperator &op, const IncomingProbe( - op, - n_miss, - [&](size_t j) -> const Monomial & { return pr.mono[pr.miss_g[j]]; }, - [&](size_t j, size_t base) { assign_row(*op.store, base + j, pr.mono[pr.miss_g[j]]); }); + // insert_absent_terms' three steps without its two dense round-trips, and on the same ordering + // contract, which is what matters: slot j lands at base+j, in miss order = (sender, record) order. + const size_t base = op.store->grow_rows_geometric(n_miss); + for (size_t j = 0; j < n_miss; ++j) { + const size_t g = pr.miss_g[j]; + op.store->set_positions(base + j, pr.pos_flat.data() + pr.pos_off[g], pr.k_of[g]); + } + op.store->bulk_insert_hashed(n_miss, base, [&](size_t j) { return pr.hash_of[pr.miss_g[j]]; }); + op.reindex_after_growth(base, n_miss); } // resolve_incoming / process_responses are the picture-independent cross-rank exchange skeletons; what @@ -134,7 +187,8 @@ 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 IncomingProbe pr = probe_incoming_queries(incoming, op, rank_count); + const IncomingProbe pr = + probe_incoming_queries(incoming, op, rank_count, sink.incoming_layout()); 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 d5b8a77e..59c1a030 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Scan.h +++ b/cpp/monoprop/detail/evolution/layer_build/Scan.h @@ -21,6 +21,7 @@ #include #include #include +#include #include #include "monoprop/TypeAliases.h" @@ -29,6 +30,7 @@ #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/QueryCodec.h" #include "monoprop/detail/graph_encoding/MPGraphEncodingTypes.h" #include "monoprop/detail/mpi/MPIUtils.h" #include "monoprop/detail/operator/InvertedIndex.h" @@ -268,6 +270,40 @@ auto fused_find_and_collect(const MPOperator &op, auto &fs = res.follower_src; auto &fv = res.follower_val; + const OperatorIndex &ham = *op.store; + + // Everything after a term survives the structural cutoff, from the dense partner emit built. + auto push = [&](const Monomial &dense, + size_t mono_pop, + size_t overlap, + int phase_factor, + size_t i, + double v_src, + bool is_follower) { + 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. + // Must be the SAME function find_rank computes (MPIUtils.h) or a term is placed and queried + // on different ranks, which duplicates a row silently; mpi_utils_tests.cpp asserts it. + size_t r_prime = my_rank; + if (rank_count != 1) { + r_prime = monomial_hash(dense) % rank_count; + } + if (is_follower) { + QueryCodec::push(fq[r_prime], dense, phase); + fs[r_prime].push_back(i); + if (capture_values) { + fv[r_prime].push_back(v_src); + } + } + else { + QueryCodec::push(lq[r_prime], dense, phase); + ls[r_prime].push_back(i); + if (capture_values) { + lv[r_prime].push_back(v_src); + } + } + }; + // The dynamic gate runs before emit_term_products, so a gate-rejected term computes no products. // abs_c/v_src come from the caller's coeff read, not re-read. auto emit = [&](size_t mono_pop, size_t i, double abs_c, double v_src, bool is_follower) { @@ -277,31 +313,16 @@ auto fused_find_and_collect(const MPOperator &op, Monomial new_mono; size_t overlap = 0; int phase_factor = 0; - emit_term_products(*op.store, i, ectx, new_mono, overlap, phase_factor); + emit_term_products(ham, i, ectx, new_mono, overlap, phase_factor); // 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; - const bool struct_pass = cutoff_eval.passes_with_popcount(new_mono, new_pop); + // nullopt only for an opaque cutoff_fn_, which has no (k, d) form and must be invoked. + const auto keep = cutoff_eval.passes_from_dense(new_mono, new_pop); + const bool struct_pass = keep.value_or(false) || (!keep.has_value() && cutoff_eval(new_mono)); if (!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 source = i; - if (is_follower) { - query_push(fq[r_prime], new_mono, 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); - ls[r_prime].push_back(source); - if (capture_values) { - lv[r_prime].push_back(v_src); - } - } + push(new_mono, mono_pop, overlap, phase_factor, i, v_src, is_follower); }; // Pass 1 and pass 2 stay fused over `nz`: splitting them regressed measurably, as `nz` spills L1 @@ -327,9 +348,11 @@ auto fused_find_and_collect(const MPOperator &op, n_foll); } if (rank_count == 1) { - lq[my_rank].reserve((n_anti - n_foll) * kQueryWords); + // A hint only: 2 words covers k <= 6; wider terms grow the buffer rather than be reserved for. + const size_t qw = QueryCodec::kReserveWordsPerQuery; + lq[my_rank].reserve((n_anti - n_foll) * qw); ls[my_rank].reserve(n_anti - n_foll); - fq[my_rank].reserve(n_foll * kQueryWords); + fq[my_rank].reserve(n_foll * qw); fs[my_rank].reserve(n_foll); } auto derive_coeff = [&](size_t i) -> std::pair { diff --git a/cpp/monoprop/detail/evolution/layer_build/SparseQuery.h b/cpp/monoprop/detail/evolution/layer_build/SparseQuery.h new file mode 100644 index 00000000..330a182c --- /dev/null +++ b/cpp/monoprop/detail/evolution/layer_build/SparseQuery.h @@ -0,0 +1,401 @@ +// 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 +#include +#include + +#include "monoprop/TypeAliases.h" + +namespace monoprop::detail { + +// A variable-width query record: one word-aligned record per term holding the term's ascending set-bit +// positions. push() picks the argmin of three closed forms -- FIXED (k lanes of kPosBits), GAP (first +// position raw, then k-1 gaps of gw bits) and BITMAP (a raw kBits mask) -- so the record is never larger +// than any of the three, which is also what bounds Writer's kMaxWords. Record order is preserved +// everywhere, because it is the floating-point accumulation order (Resolve.h mints misses in it). +// +// Header in the low bits of word 0, then the payload LSB-first, both by explicit shift, never punning: +// [0..1] mode (FIXED/GAP/BITMAP) [2..3] phase+1 (emit_phase is TERNARY) [4..9] k, 63 escaping to a +// following 16-bit k then [4 bits] gw, in GAP mode only then payload. +// Header width is a design variable, not overhead: it decides which side of a 64-bit line a term falls +// on. +template +struct SparseQuery { + using PosT = uint16_t; + + static constexpr size_t kBits = 2 * NumModes; + static_assert(kBits <= 65535, "a physical bit position and the popcount must both fit a uint16_t"); + + //: Bits for one raw position in [0, 2*NumModes); compile-time, so the lane width is free. + static constexpr size_t kPosBits = static_cast(std::bit_width(kBits - 1)); + + static constexpr uint64_t kModeFixed = 0; + static constexpr uint64_t kModeGap = 1; + static constexpr uint64_t kModeBitmap = 2; + + static constexpr size_t kModeBits = 2; + static constexpr size_t kPhaseBits = 2; + static constexpr size_t kKBits = 6; + static constexpr size_t kLongKBits = 16; + static constexpr size_t kGwBits = 4; + static constexpr size_t kKEscape = (1U << kKBits) - 1U; + static constexpr size_t kBaseHeaderBits = kModeBits + kPhaseBits + kKBits; + static_assert(kPosBits <= (1U << kGwBits) - 1U, "gw <= kPosBits must fit the header's 4-bit gap-width field"); + + static constexpr size_t kMaxPositions = 65535; + //: Words in a full occupancy mask. ceil, NOT kBits/64: a C++ caller may pick any width (LiH: 24). + static constexpr size_t kMaskWords = (kBits + 63U) / 64U; + //: BITMAP bounds every mode from above, because the encoder takes the minimum of the three. + static constexpr size_t kMaxWords = (kBits + kBaseHeaderBits + kLongKBits + kGwBits + 63U) / 64U + 1U; + + // ---- bit stream ------------------------------------------------------------------------------- + + struct Writer { + uint64_t w[kMaxWords] = {}; + size_t nbits = 0; + + constexpr auto put(uint64_t v, size_t width) noexcept -> void { + if (width == 0) { + assert(v == 0 && "a zero-width field cannot carry a value"); + return; + } + // Assert BEFORE masking: masking alone turns an overflow into a different well-formed record. + assert((width >= 64 || (v >> width) == 0) && "field value does not fit its width"); + if (width < 64) { + v &= (uint64_t{1} << width) - 1U; + } + const size_t word = nbits >> 6U; + const size_t off = nbits & 63U; + assert(word < kMaxWords && "record overran its worst-case word bound"); + w[word] |= v << off; + // off > 0 is implied here (width <= 64), so the shift is in [1, 63]; `v >> 64` would be UB. + if (off + width > 64) { + w[word + 1] |= v >> (64U - off); + } + nbits += width; + } + }; + + struct Reader { + const VecZ &buf; + size_t base; // word offset of the record start + size_t nbits = 0; + + [[nodiscard]] auto get(size_t width) noexcept -> uint64_t { + if (width == 0) { + return 0; + } + const size_t word = nbits >> 6U; + const size_t off = nbits & 63U; + uint64_t v = static_cast(buf[base + word]) >> off; + if (off + width > 64) { + v |= static_cast(buf[base + word + 1]) << (64U - off); + } + nbits += width; + return (width < 64) ? (v & ((uint64_t{1} << width) - 1U)) : v; + } + }; + + // ---- header ----------------------------------------------------------------------------------- + + struct Header { + uint64_t mode = 0; + int phase = 0; + size_t k = 0; + size_t gw = 0; + size_t bits = 0; // header width, i.e. where the payload begins + }; + + // One word load and a few masks; deliberately does NOT touch the payload -- the cursor walks call it. + [[nodiscard]] static auto header_at(const VecZ &buf, size_t off) noexcept -> Header { + const auto w0 = static_cast(buf[off]); + Header h; + h.mode = w0 & 0x3U; + h.phase = static_cast((w0 >> kModeBits) & 0x3U) - 1; + h.k = static_cast((w0 >> (kModeBits + kPhaseBits)) & kKEscape); + h.bits = kBaseHeaderBits; + if (h.k == kKEscape) { + h.k = static_cast((w0 >> h.bits) & 0xFFFFU); + h.bits += kLongKBits; + } + if (h.mode == kModeGap) { + h.gw = static_cast((w0 >> h.bits) & 0xFU); + h.bits += kGwBits; + } + return h; + } + + [[nodiscard]] static constexpr auto header_bits_for(size_t k, uint64_t mode) noexcept -> size_t { + return kBaseHeaderBits + ((k >= kKEscape) ? kLongKBits : 0U) + ((mode == kModeGap) ? kGwBits : 0U); + } + + [[nodiscard]] static constexpr auto fixed_bits(size_t k) noexcept -> size_t { + return header_bits_for(k, kModeFixed) + k * kPosBits; + } + [[nodiscard]] static constexpr auto gap_bits(size_t k, size_t gw) noexcept -> size_t { + return header_bits_for(k, kModeGap) + ((k == 0) ? 0U : kPosBits + (k - 1U) * gw); + } + [[nodiscard]] static constexpr auto bitmap_bits(size_t k) noexcept -> size_t { + return header_bits_for(k, kModeBitmap) + kBits; + } + [[nodiscard]] static constexpr auto words_of(size_t bits) noexcept -> size_t { return (bits + 63U) / 64U; } + + //: The record's word count, from the header alone: k does not determine it, mode and gw do too. + [[nodiscard]] static constexpr auto words_of_header(const Header &h) noexcept -> size_t { + switch (h.mode) { + case kModeGap: + return words_of(gap_bits(h.k, h.gw)); + case kModeBitmap: + return words_of(bitmap_bits(h.k)); + default: + return words_of(fixed_bits(h.k)); + } + } + + [[nodiscard]] static auto words_at(const VecZ &buf, size_t off) noexcept -> size_t { + return words_of_header(header_at(buf, off)); + } + + [[nodiscard]] static auto k_at(const VecZ &buf, size_t off) noexcept -> size_t { return header_at(buf, off).k; } + [[nodiscard]] static auto phase_at(const VecZ &buf, size_t off) noexcept -> int { + return header_at(buf, off).phase; + } + + // ---- encode ----------------------------------------------------------------------------------- + + //: gw = bit_width(max gap). Folded into the caller's single pass in push(), never a second walk. + template + [[nodiscard]] static auto gap_width(const PosU *pos, size_t k) noexcept -> size_t { + size_t g = 0; + for (size_t j = 1; j < k; ++j) { + const size_t d = static_cast(pos[j] - pos[j - 1] - 1U); + const auto b = static_cast(std::bit_width(d)); + g = (b > g) ? b : g; + } + return g; + } + + // Precondition: k STRICTLY ASCENDING physical bit positions in [0, kBits). A violation is silent in + // release -- gap coding is meaningless without it and an out-of-range position decodes to a different + // valid-looking monomial. Returns the WORDS written; PosU is generic because the store's position + // type is narrower than the wire's below 129 modes, and the encoding does not depend on it. + template + static auto push(VecZ &buf, const PosU *pos, size_t k, int phase) -> size_t { + assert(k <= kMaxPositions && "term has more positions than the record's k field can hold"); + assert(phase >= -1 && phase <= 1 && "emit_phase is ternary: rotation_sign, or REAL_PARTS entry"); + for (size_t j = 1; j < k; ++j) { + assert(pos[j] > pos[j - 1] && "positions must be strictly ascending"); + } + + const size_t gw = gap_width(pos, k); + const size_t wf = words_of(fixed_bits(k)); + const size_t wg = words_of(gap_bits(k, gw)); + const size_t wb = words_of(bitmap_bits(k)); + + uint64_t mode = kModeFixed; + size_t want = wf; + if (wg < want) { + mode = kModeGap; + want = wg; + } + if (wb < want) { + mode = kModeBitmap; + want = wb; + } + + Writer w; + w.put(mode, kModeBits); + w.put(static_cast(phase + 1), kPhaseBits); + if (k >= kKEscape) { + w.put(kKEscape, kKBits); + w.put(static_cast(k), kLongKBits); + } + else { + w.put(static_cast(k), kKBits); + } + + if (mode == kModeGap) { + w.put(static_cast(gw), kGwBits); + if (k != 0) { + w.put(static_cast(pos[0]), kPosBits); + for (size_t j = 1; j < k; ++j) { + w.put(static_cast(pos[j] - pos[j - 1] - 1U), gw); + } + } + } + else if (mode == kModeBitmap) { + // The trailing partial word is NOT optional: at kBits 24 a `kBits / 64` loop writes nothing + // while bitmap_bits() still charges 24 bits, and the decoder then recovers no positions. + uint64_t mask[kMaskWords] = {}; + for (size_t j = 0; j < k; ++j) { + const auto p = static_cast(pos[j]); + mask[p >> 6U] |= uint64_t{1} << (p & 63U); + } + for (size_t done = 0, i = 0; done < kBits; ++i) { + const size_t chunk = (kBits - done < 64U) ? (kBits - done) : 64U; + w.put(mask[i], chunk); + done += chunk; + } + } + else { + for (size_t j = 0; j < k; ++j) { + w.put(static_cast(pos[j]), kPosBits); + } + } + + assert(words_of(w.nbits) == want && "encoder wrote a different width than it costed"); + for (size_t i = 0; i < want; ++i) { + buf.push_back(static_cast(w.w[i])); + } + return want; + } + + // ---- decode ----------------------------------------------------------------------------------- + + // OutT is generic so the resolve path decodes straight into the store's (narrower) position width. + template + static auto read_positions(const VecZ &buf, size_t off, OutT *out) -> size_t { + const Header h = header_at(buf, off); + Reader r{buf, off, h.bits}; + if (h.mode == kModeGap) { + if (h.k != 0) { + auto prev = static_cast(r.get(kPosBits)); + out[0] = static_cast(prev); + for (size_t j = 1; j < h.k; ++j) { + prev += static_cast(r.get(h.gw)) + 1U; + out[j] = static_cast(prev); + } + } + } + else if (h.mode == kModeBitmap) { + // Symmetric with the encoder, trailing partial word included -- see the note there. + size_t n = 0; + for (size_t done = 0, i = 0; done < kBits; ++i) { + const size_t chunk = (kBits - done < 64U) ? (kBits - done) : 64U; + uint64_t word = r.get(chunk); + while (word != 0) { + const auto b = static_cast(std::countr_zero(word)); + out[n++] = static_cast(done + b); + word &= word - 1U; + } + done += chunk; + } + assert(n == h.k && "bitmap popcount disagrees with the record's k"); + } + else { + for (size_t j = 0; j < h.k; ++j) { + out[j] = static_cast(r.get(kPosBits)); + } + } + const size_t next = off + words_of_header(h); + assert(check_header(buf, off, out) && "record header is inconsistent with its own positions"); + return next; + } + + // Debug-only: every wire field must be checkable from the rest of the record, or it rots. + template + [[nodiscard]] static auto check_header(const VecZ &buf, size_t off, const OutT *pos) -> bool { + const Header h = header_at(buf, off); + if (h.mode > kModeBitmap) { + return false; // an unknown mode is a format the reader does not understand + } + if (h.phase < -1 || h.phase > 1) { + return false; + } + for (size_t j = 0; j + 1 < h.k; ++j) { + if (static_cast(pos[j]) >= static_cast(pos[j + 1])) { + return false; // positions must arrive strictly ascending + } + } + if (h.k != 0 && static_cast(pos[h.k - 1]) >= kBits) { + return false; + } + if (h.mode == kModeGap && h.k > 1) { + // gw is the MAXIMUM gap width: too small truncates a gap silently, too large wastes bits. + size_t g = 0; + for (size_t j = 1; j < h.k; ++j) { + const auto b = static_cast(std::bit_width(static_cast(pos[j] - pos[j - 1] - 1))); + g = (b > g) ? b : g; + } + if (g != h.gw) { + return false; + } + } + return true; + } + + // d, recomputed rather than carried: ascending order makes a pair an even position then its successor. + template + [[nodiscard]] static auto pair_count(const OutT *pos, size_t k) noexcept -> size_t { + size_t d = 0; + for (size_t j = 0; j + 1 < k; ++j) { + if ((pos[j] % 2 == 0) && (pos[j + 1] == pos[j] + 1)) { + ++d; + } + } + return d; + } + + static constexpr size_t kStackPositions = 64; + + static auto read_mono(const VecZ &buf, size_t off, Monomial &mono_out, int &phase_out) -> size_t { + const Header h = header_at(buf, off); + phase_out = h.phase; + mono_out = Monomial{}; + if (h.k <= kStackPositions) { + PosT scratch[kStackPositions]; + const size_t next = read_positions(buf, off, scratch); + for (size_t j = 0; j < h.k; ++j) { + mono_out.set(static_cast(scratch[j])); + } + assert(mono_out.count() == h.k && "decoded popcount disagrees with the record's k"); + return next; + } + std::vector scratch(h.k); + const size_t next = read_positions(buf, off, scratch.data()); + for (size_t j = 0; j < h.k; ++j) { + mono_out.set(static_cast(scratch[j])); + } + assert(mono_out.count() == h.k && "decoded popcount disagrees with the record's k"); + return next; + } + + // Encode from a dense monomial, which is what the scan holds: the partner is built densely anyway. + static auto push_mono(VecZ &buf, const Monomial &mono, int phase) -> size_t { + const size_t k = mono.count(); + if (k <= kStackPositions) { + PosT scratch[kStackPositions]; + size_t j = 0; + for (size_t b = mono.find_first(); b < mono.size(); b = mono.find_next(b)) { + scratch[j++] = static_cast(b); + } + assert(j == k && "find_first/find_next walk disagrees with count()"); + return push(buf, scratch, k, phase); + } + std::vector scratch(k); + size_t j = 0; + for (size_t b = mono.find_first(); b < mono.size(); b = mono.find_next(b)) { + scratch[j++] = static_cast(b); + } + assert(j == k && "find_first/find_next walk disagrees with count()"); + return push(buf, scratch.data(), k, phase); + } +}; + +} // namespace monoprop::detail diff --git a/cpp/monoprop/detail/operator/OperatorIndex.h b/cpp/monoprop/detail/operator/OperatorIndex.h index 8a3c29b2..51323e22 100644 --- a/cpp/monoprop/detail/operator/OperatorIndex.h +++ b/cpp/monoprop/detail/operator/OperatorIndex.h @@ -17,8 +17,10 @@ #include #include #include +#include #include #include +#include #include #include #include @@ -96,6 +98,9 @@ class OperatorIndex { [[nodiscard]] auto size() const -> size_t { return size_; } + // Rows that exceeded inline_width_ and spilled; observable so a test can compare the two insert paths. + [[nodiscard]] auto overflow_size() const -> size_t { return overflow_.size(); } + auto reserve(size_t n) -> void { reserve_rows(n); reserve_index(n); @@ -137,6 +142,33 @@ class OperatorIndex { } } + // set() from the row's own form: a row IS an ascending position list. Same postcondition as set(), + // including the dropped stale overflow entry. + // + // Precondition: `pos` strictly ascending, every entry < 2*NumModes. A violation is silent in release + // -- an unsorted row simply never matches, and an out-of-range one decodes to a different term. + auto set_positions(size_t i, const PosT *pos, size_t count) -> void { + assert(std::adjacent_find(pos, pos + count, std::greater_equal{}) == pos + count + && "row positions must be strictly ascending"); + assert((count == 0 || static_cast(pos[count - 1]) < 2 * NumModes) && "row position out of range"); + PosT *row = &rows_[i * stride_]; + if (count > inline_width_) { + // The spill path has no position array, so build the dense form -- only here. + row[0] = kOverflowMarker; + value_type mono; + for (size_t j = 0; j < count; ++j) { + mono.set(pos[j]); + } + overflow_[i] = mono; + return; + } + if (!overflow_.empty()) { + overflow_.erase(i); + } + row[0] = static_cast(count); + std::copy_n(pos, count, row + 1); + } + [[nodiscard]] auto row(size_t i) const -> value_type { const PosT c = rows_[i * stride_]; if (c == kOverflowMarker) { @@ -170,6 +202,19 @@ class OperatorIndex { } return overflow_.at(i).count(); } + // The row's stored ascending positions; (nullptr, 0) for a spilled row, and invalidated by any insert. + struct RowPositions { + const PosT *pos; + size_t count; + [[nodiscard]] auto inlined() const -> bool { return pos != nullptr; } + }; + [[nodiscard]] auto row_positions(size_t i) const -> RowPositions { + const PosT c = rows_[i * stride_]; + if (c == kOverflowMarker) { + return {nullptr, 0}; + } + return {&rows_[(i * stride_) + 1], static_cast(c)}; + } [[nodiscard]] auto memory_bytes() const -> size_t { size_t total = rows_.capacity() * sizeof(PosT); total += overflow_.size() * (sizeof(value_type) + sizeof(size_t) + 24); @@ -233,6 +278,66 @@ class OperatorIndex { } } + // find_batch over ascending position lists: query q is pos_flat[pos_off[q] .. pos_off[q] + k_of[q]). + // Identical results to find_batch on the monomials those positions describe. Same three-stage + // prefetch pipeline, so the positions stay the currency without giving up find_batch's shape. + auto find_batch_positions(const PosT *pos_flat, + const size_t *pos_off, + const uint32_t *k_of, + size_t n, + size_t *out, + uint32_t *hash_out = nullptr) const -> void { + static constexpr size_t G = 16; + std::array hh; + std::array sp; + std::array cand; + for (size_t base = 0; base < n; base += G) { + const size_t g = std::min(G, n - base); + for (size_t j = 0; j < g; ++j) { + hh[j] = fold_hash_positions(pos_flat + pos_off[base + j], k_of[base + j]); + sp[j] = spread(hh[j]); + __builtin_prefetch(&table_.slots[sp[j] & table_.mask], 0, 0); + } + if (hash_out != nullptr) { + std::copy_n(hh.begin(), g, hash_out + base); + } + for (size_t j = 0; j < g; ++j) { + cand[j] = kEmptySlot; + if (table_.count == 0) { + continue; + } + cand[j] = probe_hash_match_(hh[j], sp[j] & table_.mask); + if (cand[j] != kEmptySlot) { + __builtin_prefetch(&rows_[static_cast(cand[j]) * stride_], 0, 0); + } + } + for (size_t j = 0; j < g; ++j) { + const size_t q = base + j; + const PosT *qpos = pos_flat + pos_off[q]; + const size_t qk = k_of[q]; + if (cand[j] == kEmptySlot) { + out[q] = kNotFound; + } + else if (row_eq_positions(static_cast(cand[j]), qpos, qk)) { + out[q] = static_cast(cand[j]); + } + else { + // A 32-bit collision: rare enough to walk the chain from the top rather than resume it. + out[q] = find_positions_(hh[j], qpos, qk); + } + } + } + } + + // fold_hash of the monomial `pos` describes, through the same fold, so it is equal by construction. + [[nodiscard]] static auto fold_hash_positions(const PosT *pos, size_t count) noexcept -> uint32_t { + key_type mono; + for (size_t j = 0; j < count; ++j) { + mono.set(pos[j]); + } + return fold_hash(mono); + } + // Insert-or-no-op. Row at `value` must already be written (the confirm reads dense rows). auto emplace(const key_type &key, mapped_type value) -> void { check_index_fits(value); @@ -251,12 +356,35 @@ class OperatorIndex { // Insert n distinct rows with consecutive indices [base, base+n). Rows must already be written. template auto bulk_insert(size_t n, mapped_type base, KeyFn &&key_at) -> void { + if (n == 0) { + return; + } + // Delegating means both entry points share one insert loop, prefetch pipeline included. + bulk_insert_hashed(n, base, [&](size_t k) { return fold_hash(key_at(k)); }); + } + // bulk_insert with the hashes already in hand: same precondition (n distinct rows, already written, + // at consecutive indices) and the same slot assignment. `hashes[k]` MUST be fold_hash of the key of + // row base+k -- a wrong one leaves the row unfindable, which surfaces later as a duplicate insert. + // + // Group-prefetched like find_batch: correctness does not depend on it (a prefetch is a hint and the + // insert re-reads the slot), but hash_at is called exactly ONCE per element and buffered. + template + auto bulk_insert_hashed(size_t n, mapped_type base, HashFn &&hash_at) -> void { if (n == 0) { return; } check_index_fits(base + n - 1); - for (size_t k = 0; k < n; ++k) { - insert_slot_(static_cast(base + k), fold_hash(key_at(k))); + static constexpr size_t G = 16; // same group width as find_batch, for the same reason + std::array hh; + for (size_t b = 0; b < n; b += G) { + const size_t g = std::min(G, n - b); + for (size_t j = 0; j < g; ++j) { + hh[j] = hash_at(b + j); + __builtin_prefetch(&table_.slots[spread(hh[j]) & table_.mask], /*rw=*/1, /*locality=*/0); + } + for (size_t j = 0; j < g; ++j) { + insert_slot_(static_cast(base + b + j), hh[j]); + } } } template @@ -386,6 +514,38 @@ class OperatorIndex { return true; } + // Compare row i against an ascending position list; a spilled row falls back to a dense compare. + [[nodiscard]] auto row_eq_positions(size_t i, const PosT *q, size_t qk) const -> bool { + const PosT c = rows_[i * stride_]; + if (c == kOverflowMarker) { + key_type mono; + for (size_t j = 0; j < qk; ++j) { + mono.set(q[j]); + } + return overflow_.at(i) == mono; + } + if (qk != static_cast(c)) { + return false; + } + return std::equal(q, q + qk, &rows_[(i * stride_) + 1]); + } + + // find()'s chain walk for a position-list key, hash already folded; only the collision arm reaches it. + [[nodiscard]] auto find_positions_(uint32_t h, const PosT *q, size_t qk) const -> size_t { + if (table_.count == 0) { + return kNotFound; + } + for (size_t s = spread(h) & table_.mask;; s = (s + 1) & table_.mask) { + const Slot &e = table_.slots[s]; + if (e.idx == kEmptySlot) { + return kNotFound; + } + if (e.h == h && row_eq_positions(static_cast(e.idx), q, qk)) { + return static_cast(e.idx); + } + } + } + static auto check_index_fits(size_t value) -> void { if (value >= kIndexCeiling) { throw TermIndexCeilingReached("OperatorIndex: operator index reached the TermIndex ceiling; rebuild with " diff --git a/cpp/tests/README.md b/cpp/tests/README.md index 6fc34490..6cd35461 100644 --- a/cpp/tests/README.md +++ b/cpp/tests/README.md @@ -73,6 +73,9 @@ name and cannot address suite-nested cases, tests use flat - **`ExchangeLayoutOracle.h`**: `build_layer_exchange_layout` — the independent reference for a layer's exchange counts and displacements, which `derive_exchange_layout` in the library is checked against. +- **`dense_query_reference.h`**: the retired dense query record, frozen as the + independent oracle for `sparse_query_tests.cpp`. Test-only, and not kept in + sync with the wire format. - **`TestData.{h,cpp}`**: the `CaseData` struct and msgpack fixture loader. - **`boost-test.cmake` / `boostAddTests.cmake`**: CMake test discovery. @@ -83,17 +86,26 @@ name and cannot address suite-nested cases, tests use flat vs a std::bitset oracle), `mpfunctions.cpp` (MP utilities + bit-flip helpers), `pauli_algebra_tests.cpp`, `majorana_cutoff_tests.cpp` (length/support cutoff, CutoffEvaluator, interleave phase, coeff encode/decode), `validation_tests.cpp` - (parameter validators), `mpi_utils_tests.cpp` (find_rank + word serialization), - `evolution_detail_tests.cpp` (MatchedEpochSet + CutoffContext), - `row_accessor_tests.cpp` (dense vs OperatorIndex row accessors). + (parameter validators), `mpi_utils_tests.cpp` (find_rank, word serialization, + scan routing agreement), `evolution_detail_tests.cpp` (MatchedEpochSet + + CutoffContext), + `row_accessor_tests.cpp` (dense vs OperatorIndex row accessors), + `sparse_monomial_tests.cpp` (the `(k, d)` cutoff predicates vs their bitset + forms). - **Operator store**: `operator_index_tests.cpp`, `inverted_index_tests.cpp`, `mp_operator_tests.cpp` (MPOperator get_state Pauli/Majorana scoring, get_operator init-map drain, update_initial_operator picture branches, - insert_absent_terms, inverted-index sync, memory estimate, deep copy). + insert_absent_terms, inverted-index sync, memory estimate, deep copy), + `bulk_insert_tests.cpp` (the grouped-prefetch insert vs a one-key-at-a-time + reference: table state and enumeration order). - **Layer build / evolution**: `build_graph_tests.cpp`, `pauli_build_layer_tests.cpp`, `fused_cos_sweep_tests.cpp`, - `fused_query_codec_tests.cpp`, `combined_recompute_equivalence.cpp` - (recompute equivalence + snapshot invariance), `exact_upper_atol_rescue.cpp`, + `sparse_query_tests.cpp` (the SparseQuery wire record against the frozen dense + oracle, plus the fused value channel), `sparse_resolve_tests.cpp` (probe and + insert from wire positions vs the dense Monomial-keyed path), + `digest_cutoff_tests.cpp` (paired_mode_count and the digest cutoff predicate + vs cutoff_sums), `combined_recompute_equivalence.cpp` (recompute equivalence + + snapshot invariance), `exact_upper_atol_rescue.cpp`, `large_cosine_storage_tests.cpp`, `gate_boundaries.cpp`. - **Graph encoding / packing**: `graph_encoding_tests.cpp` (CosineWordBuilder coalescer, checked_* overflow guards, packed-phase storage + int8 read, diff --git a/cpp/tests/bulk_insert_tests.cpp b/cpp/tests/bulk_insert_tests.cpp new file mode 100644 index 00000000..2221b340 --- /dev/null +++ b/cpp/tests/bulk_insert_tests.cpp @@ -0,0 +1,188 @@ +// 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. + +// bulk_insert prefetches 16 slot addresses; being a pure hint it must leave the table in EXACTLY the +// state an unpipelined loop leaves it in, slot order included -- for_each makes it Python-visible. + +#include + +#include +#include +#include +#include +#include +#include + +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/operator/OperatorIndex.h" + +using namespace monoprop; + +namespace { + +constexpr size_t kN = 250; +using Index = detail::OperatorIndex; + +auto draw_distinct(std::mt19937_64 &rng, size_t n) -> std::vector> { + std::vector> out; + std::set> seen; + std::uniform_int_distribution bit(0, Monomial::size() - 1); + std::uniform_int_distribution pop(0, 12); + while (out.size() < n) { + Monomial m; + const size_t k = pop(rng); + for (size_t placed = 0; placed < k;) { + const size_t b = bit(rng); + if (!m.test(b)) { + m.set(b); + ++placed; + } + } + std::vector key; + for (size_t w = 0; w < Monomial::num_words(); ++w) { + key.push_back(m.word(w)); + } + if (seen.insert(key).second) { + out.push_back(m); + } + } + return out; +} + +// Deliberately NOT reserved: rehash_if_needed firing mid-group frees the table already-issued +// addresses point into, the interesting case for a prefetch. +auto build(const std::vector> &terms) -> std::unique_ptr { + auto idx = std::make_unique(); + const size_t base = idx->grow_rows_geometric(terms.size()); + for (size_t k = 0; k < terms.size(); ++k) { + idx->set(base + k, terms[k]); + } + idx->bulk_insert(terms.size(), base, [&](size_t k) -> const Monomial & { return terms[k]; }); + return idx; +} + +// The oracle: N one-key calls, so every group is of size one and none of the grouped loop's boundary +// arithmetic runs. It shares insert_slot_ but not the GROUPING, which is the thing under test. +auto build_reference(const std::vector> &terms) -> std::unique_ptr { + auto idx = std::make_unique(); + const size_t base = idx->grow_rows_geometric(terms.size()); + for (size_t k = 0; k < terms.size(); ++k) { + idx->set(base + k, terms[k]); + idx->bulk_insert(1, base + k, [&](size_t) -> const Monomial & { return terms[k]; }); + } + return idx; +} + +} // namespace + +BOOST_AUTO_TEST_CASE(bulk_insert_finds_every_key) { + std::mt19937_64 rng(20260814); + const auto terms = draw_distinct(rng, 4000); + const auto idx = build(terms); + BOOST_REQUIRE_EQUAL(idx->size(), terms.size()); + for (size_t i = 0; i < terms.size(); ++i) { + const auto found = idx->find(terms[i]); + BOOST_REQUIRE(found.has_value()); + BOOST_TEST(*found == i); + } +} + +BOOST_AUTO_TEST_CASE(bulk_insert_batch_find_agrees_with_one_key_at_a_time) { + std::mt19937_64 rng(20260815); + const auto terms = draw_distinct(rng, 4000); + const auto absent = draw_distinct(rng, 500); + const auto grouped = build(terms); + const auto one_by_one = build_reference(terms); + + std::vector a(terms.size(), 0); + std::vector b(terms.size(), 0); + one_by_one->find_batch(terms.data(), terms.size(), a.data()); + grouped->find_batch(terms.data(), terms.size(), b.data()); + for (size_t i = 0; i < terms.size(); ++i) { + BOOST_TEST(a[i] == b[i]); + } + + // A table that answered everything would satisfy the loop above and prove nothing. + std::vector ma(absent.size(), 0); + std::vector mb(absent.size(), 0); + one_by_one->find_batch(absent.data(), absent.size(), ma.data()); + grouped->find_batch(absent.data(), absent.size(), mb.data()); + size_t genuinely_absent = 0; + for (size_t i = 0; i < absent.size(); ++i) { + BOOST_TEST(ma[i] == mb[i]); + if (!one_by_one->find(absent[i]).has_value()) { + ++genuinely_absent; + } + } + BOOST_TEST(genuinely_absent > 0U); +} + +BOOST_AUTO_TEST_CASE(bulk_insert_preserves_the_enumeration_order) { + // Identical, not merely equal as a set: this sequence is what the Python API returns terms in. + std::mt19937_64 rng(20260816); + const auto terms = draw_distinct(rng, 4000); + const auto grouped = build(terms); + const auto one_by_one = build_reference(terms); + + std::vector order_ref; + std::vector order_grouped; + one_by_one->for_each([&](const Monomial &, size_t i) { order_ref.push_back(i); }); + grouped->for_each([&](const Monomial &, size_t i) { order_grouped.push_back(i); }); + + BOOST_REQUIRE_EQUAL(order_grouped.size(), order_ref.size()); + BOOST_REQUIRE_EQUAL(order_ref.size(), terms.size()); + // Both coming out in index order would make the comparison blind to a reshuffle preserving it. + bool is_sorted_by_index = true; + for (size_t i = 1; i < order_ref.size(); ++i) { + if (order_ref[i] < order_ref[i - 1]) { + is_sorted_by_index = false; + break; + } + } + BOOST_TEST(!is_sorted_by_index); + for (size_t i = 0; i < order_ref.size(); ++i) { + BOOST_TEST(order_grouped[i] == order_ref[i]); + } +} + +BOOST_AUTO_TEST_CASE(bulk_insert_handles_a_partial_final_group) { + // The group width is 16, so these sizes straddle the boundary, including fewer than one group. + std::mt19937_64 rng(20260817); + for (const size_t n : {size_t{1}, size_t{15}, size_t{16}, size_t{17}, size_t{31}, size_t{33}}) { + const auto terms = draw_distinct(rng, n); + const auto grouped = build(terms); + const auto one_by_one = build_reference(terms); + BOOST_REQUIRE_EQUAL(grouped->size(), n); + for (size_t i = 0; i < n; ++i) { + const auto fa = one_by_one->find(terms[i]); + const auto fb = grouped->find(terms[i]); + BOOST_REQUIRE(fa.has_value()); + BOOST_REQUIRE(fb.has_value()); + BOOST_TEST(*fb == *fa); + } + } +} + +BOOST_AUTO_TEST_CASE(bulk_insert_of_nothing_is_a_no_op) { + std::mt19937_64 rng(20260818); + const auto terms = draw_distinct(rng, 100); + auto idx = build(terms); + const size_t before = idx->size(); + const auto key_at = [&](size_t k) -> const Monomial & { return terms[k]; }; + idx->bulk_insert(0, 0, key_at); + BOOST_TEST(idx->size() == before); + for (size_t i = 0; i < terms.size(); ++i) { + BOOST_REQUIRE(idx->find(terms[i]).has_value()); + } +} diff --git a/cpp/tests/dense_query_reference.h b/cpp/tests/dense_query_reference.h new file mode 100644 index 00000000..f190ce26 --- /dev/null +++ b/cpp/tests/dense_query_reference.h @@ -0,0 +1,72 @@ +// 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 retired dense query record, frozen here as a deliberately independent oracle for +// sparse_query_tests.cpp's differential: test-only, and it does NOT follow the wire format. + +#pragma once + +#include +#include + +#include "monoprop/TypeAliases.h" +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/evolution/layer_build/Common.h" +#include "monoprop/detail/mpi/MPIUtils.h" + +namespace monoprop::test_ref { + +using monoprop::VecZ; + +// W monomial words + one ±1 phase word. +template +inline constexpr size_t kQueryWords = mpi_detail::kWords + 1; + +template +inline constexpr size_t kQueryWordsFused = kQueryWords + 1; + +template +inline auto query_push(VecZ &buf, const Monomial &mono, int phase) -> void { + mpi_detail::append_monomial_words(mono, buf); + buf.push_back(static_cast(static_cast(phase))); +} + +template > +inline auto query_read(const VecZ &buf, size_t q, Monomial &mono_out, int &phase_out) -> void { + const size_t base = q * QW; + mono_out = mpi_detail::read_monomial_from_words(buf, base); + phase_out = static_cast(static_cast(buf[base + mpi_detail::kWords])); +} + +template +inline auto query_value(const VecZ &buf, size_t q) -> double { + return detail::decode_value(buf[q * kQueryWordsFused + mpi_detail::kWords + 1]); +} + +// Requires v.size() == q.size()/kQueryWords: exactly one value per query record. +template +inline auto build_fused_query_value(const VecZ &q, const std::vector &v, VecZ &out) -> void { + constexpr size_t W = kQueryWords; + const size_t nq = q.empty() ? 0 : q.size() / W; + out.clear(); + out.reserve(nq * kQueryWordsFused); + for (size_t i = 0; i < nq; ++i) { + out.insert(out.end(), + q.begin() + static_cast(i * W), + q.begin() + static_cast((i + 1) * W)); + out.push_back(detail::encode_value(v[i])); + } +} + +} // namespace monoprop::test_ref diff --git a/cpp/tests/digest_cutoff_tests.cpp b/cpp/tests/digest_cutoff_tests.cpp new file mode 100644 index 00000000..b34a6114 --- /dev/null +++ b/cpp/tests/digest_cutoff_tests.cpp @@ -0,0 +1,143 @@ +// 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. + +// paired_mode_count's d, differentially against the bitset cutoff_sums: the IDENTITY +// (d == popcount_sum - or_sum) and the PREDICATE built on it are separable, so they are separate cases. + +#include + +#include +#include +#include +#include +#include + +#include "monoprop/algebra/Algebra.h" +#include "monoprop/algebra/AlgebraCommon.h" +#include "monoprop/detail/evolution/CutoffContext.h" +#include "monoprop/detail/evolution/layer_build/Scan.h" +#include "monoprop/detail/operator/MPOperator.h" + +using namespace monoprop; + +namespace { + +// indices_to_bitset places every bit at or above the active offset: the precondition inherited here. +template +auto draw_well_formed(std::mt19937_64 &rng, size_t logical, size_t weight) -> Monomial { + VecZ idx; + std::uniform_int_distribution dist(0, (2 * logical) - 1); + while (idx.size() < weight) { + const size_t v = dist(rng); + bool dup = false; + for (const auto x : idx) { + dup = dup || (x == v); + } + if (!dup) { + idx.push_back(v); + } + } + return indices_to_bitset(idx); +} + +// Weight spans both extremes: d == 0 and k == 2d are the cases the cutoffs branch on. +template +auto check_identity(std::mt19937_64 &rng, size_t logical, size_t &checked, size_t &paired_seen) -> void { + for (size_t weight = 1; weight <= 2 * logical && weight <= 24; ++weight) { + for (int rep = 0; rep < 40; ++rep) { + const auto mono = draw_well_formed(rng, logical, weight); + const auto sums = cutoff_sums(mono, logical); + const size_t d = paired_mode_count(mono); + BOOST_REQUIRE_EQUAL(d, sums.popcount_sum - sums.or_sum); + const auto rebuilt = cutoff_sums(sums.popcount_sum, d); + BOOST_REQUIRE_EQUAL(rebuilt.xor_sum, sums.xor_sum); + BOOST_REQUIRE_EQUAL(rebuilt.or_sum, sums.or_sum); + BOOST_REQUIRE_EQUAL(rebuilt.popcount_sum, sums.popcount_sum); + paired_seen += static_cast(sums.xor_sum == 0); + ++checked; + } + } +} + +} // namespace + +BOOST_AUTO_TEST_CASE(paired_mode_count_matches_cutoff_sums_across_widths) { + std::mt19937_64 rng(0xD16E57U); + size_t checked = 0; + size_t paired_seen = 0; + + check_identity<32>(rng, 32, checked, paired_seen); // W = 64, one word, no active offset + check_identity<32>(rng, 30, checked, paired_seen); // W = 64, active_bit_offset = 4 + check_identity<48>(rng, 45, checked, paired_seen); // W = 96 -- not a multiple of 64 + check_identity<64>(rng, 64, checked, paired_seen); // W = 128, exactly two words + check_identity<128>(rng, 120, checked, paired_seen); + check_identity<256>(rng, 250, checked, paired_seen); // the production shape + + BOOST_TEST(checked > 3000U); + // The identity is only interesting on the fully-paired branch, so the draw must reach it. + BOOST_TEST(paired_seen > 0U); +} + +BOOST_AUTO_TEST_CASE(paired_mode_count_exhaustive_at_small_width) { + constexpr size_t kN = 5; // W = 10 + size_t paired = 0; + for (uint64_t bits = 0; bits < (uint64_t{1} << (2 * kN)); ++bits) { + Monomial mono; + for (size_t b = 0; b < 2 * kN; ++b) { + if ((bits >> b) & 1U) { + mono.set(b); + } + } + const auto sums = cutoff_sums(mono, kN); + BOOST_REQUIRE_EQUAL(paired_mode_count(mono), sums.popcount_sum - sums.or_sum); + paired += static_cast(sums.xor_sum == 0); + } + BOOST_TEST(paired == 32U); // 2^5: each mode independently empty or doubly occupied +} + +// At the PREDICATE level, not the scan level: cutoff_sums is the independent form to compare against. +BOOST_AUTO_TEST_CASE(digest_predicate_matches_cutoff_sums_predicate) { + std::mt19937_64 rng(0xC0FFEEU); + size_t checked = 0; + size_t kept = 0; + size_t rejected = 0; + + // The popcount <= cutoff early-out is the asymmetry between the two, so cutoffs straddle it. + for (const unsigned int cutoff : {1U, 2U, 4U, 6U, 10U, 20U}) { + for (const bool support : {false, true}) { + constexpr size_t kN = 32; + constexpr size_t kLogical = 30; + const CutoffFn fn = support ? CutoffFn{detail::SupportCutoff{cutoff, kLogical}} + : CutoffFn{detail::LengthCutoff{cutoff, kLogical}}; + const detail::CutoffEvaluator eval(fn); + for (size_t w = 1; w <= 12; ++w) { + for (int rep = 0; rep < 40; ++rep) { + const auto mono = draw_well_formed(rng, kLogical, w); + const size_t k = mono.count(); + const auto digest = eval.passes_from_dense(mono, k); + BOOST_REQUIRE(digest.has_value()); // a concrete cutoff must always decide + const bool reference = eval.passes_with_popcount(mono, k); + BOOST_REQUIRE_EQUAL(*digest, reference); + ++checked; + kept += static_cast(*digest); + rejected += static_cast(!*digest); + } + } + } + } + BOOST_TEST(checked > 5000U); + // A sweep that only ever kept would agree with any predicate that returns true. + BOOST_TEST(kept > 0U); + BOOST_TEST(rejected > 0U); +} diff --git a/cpp/tests/fused_query_codec_tests.cpp b/cpp/tests/fused_query_codec_tests.cpp deleted file mode 100644 index 72f76763..00000000 --- a/cpp/tests/fused_query_codec_tests.cpp +++ /dev/null @@ -1,125 +0,0 @@ -// 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. - -#include - -#include -#include -#include - -#include "monoprop/TypeAliases.h" -#include "monoprop/core/Monomial.h" -#include "monoprop/detail/evolution/layer_build/Common.h" - -// Query+value fusion codec: the fused R>1 exchange rides the source coefficient (v_src) on each query -// record as a trailing bit-cast word so one alltoallv carries both streams. Round-tripping must be -// byte-for-byte, including the FP corner cases a lossy value channel would mangle. - -namespace { - -using namespace monoprop; -using monoprop::detail::build_fused_query_value; -using monoprop::detail::kQueryWords; -using monoprop::detail::kQueryWordsFused; -using monoprop::detail::query_push; -using monoprop::detail::query_read; -using monoprop::detail::query_value; - -constexpr size_t kModes = 8; // 2*kModes = 16 majorana bits, one 64-bit word - -// A deterministic, distinct majorana bit pattern per record index. -auto make_mono(size_t r) -> Monomial { - Monomial m; - for (size_t b = 0; b < 2 * kModes; ++b) { - if (((r * 2654435761u + b * 40503u) & 3u) == 0u) { - m.set(b); - } - } - return m; -} - -BOOST_AUTO_TEST_CASE(fused_record_roundtrip_exact) { - const std::vector phases = {1, -1, 1, -1, 1, 1, -1}; - const std::vector values = { - 0.0, - -0.0, - 1.0, - -1.0, - 3.141592653589793, - -2.718281828459045e-300, // near-denormal magnitude - std::numeric_limits::min(), // smallest normal - }; - const size_t nq = values.size(); - - VecZ plain; - std::vector> monos(nq); - for (size_t r = 0; r < nq; ++r) { - monos[r] = make_mono(r); - query_push(plain, monos[r], phases[r]); - } - BOOST_REQUIRE_EQUAL(plain.size(), nq * kQueryWords); - - VecZ fused; - build_fused_query_value(plain, values, fused); - BOOST_REQUIRE_EQUAL(fused.size(), nq * kQueryWordsFused); - - for (size_t q = 0; q < nq; ++q) { - Monomial m_out; - int ph_out = 0; - query_read>(fused, q, m_out, ph_out); - BOOST_CHECK(m_out == monos[q]); - BOOST_CHECK_EQUAL(ph_out, phases[q]); - // Compare the raw payload, so -0.0 and denormals stay distinguished from 0.0. - const double v_out = query_value(fused, q); - BOOST_CHECK(std::memcmp(&v_out, &values[q], sizeof(double)) == 0); - } -} - -// Reusing `out` across calls must leak no stale words: capacity is a high-water mark, size is exact. -// This is the reuse pattern LayerBuildEngine::combined_qv_ relies on gate to gate. -BOOST_AUTO_TEST_CASE(fused_buffer_reuse_shrinks_logical_size) { - VecZ plain_big; - std::vector vbig; - for (size_t r = 0; r < 32; ++r) { - query_push(plain_big, make_mono(r), (r % 2 == 0) ? 1 : -1); - vbig.push_back(static_cast(r) * 1.5 - 7.0); - } - VecZ out; - build_fused_query_value(plain_big, vbig, out); - const size_t cap_after_big = out.capacity(); - - VecZ plain_small; - std::vector vsmall = {42.0, -42.0, 0.25}; - for (size_t r = 0; r < vsmall.size(); ++r) { - query_push(plain_small, make_mono(100 + r), 1); - } - build_fused_query_value(plain_small, vsmall, out); - BOOST_CHECK_EQUAL(out.size(), vsmall.size() * kQueryWordsFused); - BOOST_CHECK_GE(out.capacity(), cap_after_big); - for (size_t q = 0; q < vsmall.size(); ++q) { - const double v_out = query_value(out, q); - BOOST_CHECK(std::memcmp(&v_out, &vsmall[q], sizeof(double)) == 0); - } -} - -// Empty input arises for the self slot, which resolve_self_queries clears before the exchange. -BOOST_AUTO_TEST_CASE(fused_empty_input) { - VecZ empty; - std::vector no_values; - VecZ out{1, 2, 3}; // pre-dirtied; build must clear it - build_fused_query_value(empty, no_values, out); - BOOST_CHECK(out.empty()); -} - -} // namespace diff --git a/cpp/tests/mpi_utils_tests.cpp b/cpp/tests/mpi_utils_tests.cpp index 8372e88f..0f296507 100644 --- a/cpp/tests/mpi_utils_tests.cpp +++ b/cpp/tests/mpi_utils_tests.cpp @@ -12,18 +12,29 @@ // See the License for the specific language governing permissions and // limitations under the License. -// The pure MPIUtils.h primitives (term->owner mapping, wire word packing), driven without a comm. +// The pure MPIUtils.h primitives (term->owner mapping, wire word packing), driven without a comm -- +// plus the routing agreement between find_rank and the scan, a property of neither call site alone. #include +#include +#include +#include +#include #include #include #include "monoprop/algebra/MajoranaAlgebra.h" +#include "monoprop/detail/evolution/CutoffContext.h" +#include "monoprop/detail/evolution/layer_build/Scan.h" #include "monoprop/detail/mpi/MPIUtils.h" +#include "monoprop/detail/operator/MPOperator.h" +#include "monoprop/detail/operator/OperatorIndex.h" using namespace monoprop; +// find_rank is splitmix over the dense words modulo the rank count, and nothing else, so the oracle +// is asserted unconditionally rather than as one of several permitted hashes. BOOST_AUTO_TEST_CASE(mpi_utils_find_rank_range_and_hash_mod) { constexpr size_t N = 32; std::mt19937_64 rng(0x9E3779B9ULL); @@ -36,8 +47,8 @@ BOOST_AUTO_TEST_CASE(mpi_utils_find_rank_range_and_hash_mod) { const auto mono = indices_to_bitset(inds); for (size_t n_ranks : {size_t{1}, size_t{2}, size_t{3}, size_t{7}}) { const size_t r = find_rank(mono, n_ranks); - BOOST_TEST(r < n_ranks); BOOST_TEST(r == monomial_hash(mono) % n_ranks); + BOOST_TEST(r < n_ranks); BOOST_TEST(r == find_rank(mono, n_ranks)); // deterministic } } @@ -73,3 +84,93 @@ BOOST_AUTO_TEST_CASE(mpi_utils_monomial_words_roundtrip) { BOOST_REQUIRE(sbuf.size() == mpi_detail::kWords); BOOST_TEST((mpi_detail::read_monomial_from_words(sbuf, 0) == d)); } + +namespace { + +template +auto draw_well_formed(std::mt19937_64 &rng, size_t logical, size_t weight) -> Monomial { + VecZ idx; + std::uniform_int_distribution dist(0, (2 * logical) - 1); + while (idx.size() < weight) { + const size_t v = dist(rng); + if (std::find(idx.begin(), idx.end(), v) == idx.end()) { + idx.push_back(v); + } + } + return indices_to_bitset(idx); +} + +auto build_op(const std::vector> &terms) -> detail::MPOperator<32> { + detail::MPOperator<32> op; + op.basis = Basis::Majorana; + detail::insert_absent_terms<32>( + op, + terms.size(), + [&](size_t k) -> const Monomial<32> & { return terms[k]; }, + [&](size_t k, size_t base) { assign_row<32>(*op.store, base + k, terms[k]); }); + return op; +} + +auto check_bucket_ownership(const std::vector &buckets, size_t ranks, size_t &checked) -> void { + // Every offset comes from the codec's walk: the record is VARIABLE WIDTH, so a hardcoded stride + // would compare a monomial decoded at the wrong offset against the wrong rank. + using QC = detail::QueryCodec<32>; + const detail::QueryLayout layout{/*fused=*/false}; + for (size_t r = 0; r < buckets.size(); ++r) { + size_t off = 0; + while (off < buckets[r].size()) { + Monomial<32> mono; + int phase = 0; + QC::read_mono(buckets[r], off, mono, phase); + BOOST_REQUIRE_EQUAL(find_rank<32>(mono, ranks), r); + off = QC::next_off(buckets[r], layout, off); + ++checked; + } + BOOST_REQUIRE_EQUAL(off, buckets[r].size()); + } +} + +} // namespace + +// The scan hashes the partner it just built; find_rank hashes what the resolve side decoded off the +// wire. Nothing downstream notices if they diverge -- the term simply exists twice. +BOOST_AUTO_TEST_CASE(mpi_utils_scan_routing_agrees_with_find_rank) { + constexpr size_t kN = 32; + constexpr size_t kLogical = 30; + std::mt19937_64 rng(0xB0B1E5U); + + std::vector> terms; + for (size_t i = 0; i < 2000; ++i) { + terms.push_back(draw_well_formed(rng, kLogical, 1 + (rng() % 6))); + } + auto op = build_op(terms); + const Monomial gen = draw_well_formed(rng, kLogical, 4); + VecD coeffs(op.store->size(), 1.0); + + const CutoffFn fn = detail::LengthCutoff{10, kLogical}; + const detail::CutoffEvaluator eval(fn); + const auto cut = detail::build_majorana_evolution_cutoff_state(std::nullopt, + std::cref(coeffs), + std::nullopt, + std::optional{0.3}); + + size_t checked = 0; + for (const size_t ranks : {2U, 4U, 8U}) { + const auto res = detail::fused_find_and_collect>(op, + gen, + eval, + cut, + coeffs, + std::nullopt, + ranks, + 0, + false, + nullptr, + 1.0); + BOOST_REQUIRE_EQUAL(res.leader_queries.size(), ranks); + check_bucket_ownership(res.leader_queries, ranks, checked); + check_bucket_ownership(res.follower_queries, ranks, checked); + } + // Without this the loop above passes trivially if the scan emitted nothing at all. + BOOST_TEST(checked > 1000U); +} diff --git a/cpp/tests/sparse_monomial_tests.cpp b/cpp/tests/sparse_monomial_tests.cpp new file mode 100644 index 00000000..006fac66 --- /dev/null +++ b/cpp/tests/sparse_monomial_tests.cpp @@ -0,0 +1,159 @@ +// 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 (k, d) integer predicates differentially against the dense bitset forms they displace; the emit +// path calls only these, so a disagreement is a silently wrong keep/reject, never a crash. Both +// populations are drawn on purpose: uniform draws land on the fully-paired branch 11 times in 28500. + +#include + +#include +#include +#include +#include + +#include "monoprop/algebra/Algebra.h" +#include "monoprop/algebra/AlgebraCommon.h" +#include "monoprop/algebra/MajoranaAlgebra.h" +#include "monoprop/core/SparseMonomial.h" + +using namespace monoprop; + +namespace { + +// indices_to_bitset is the only constructor user input reaches; any other draw is unreachable state. +template +auto draw(std::mt19937_64 &rng, size_t logical, size_t weight) -> Monomial { + VecZ idx; + std::uniform_int_distribution dist(0, (2 * logical) - 1); + while (idx.size() < weight) { + const size_t v = dist(rng); + bool dup = false; + for (const auto x : idx) { + dup = dup || (x == v); + } + if (!dup) { + idx.push_back(v); + } + } + return indices_to_bitset(idx); +} + +template +auto draw_paired(std::mt19937_64 &rng, size_t logical, size_t modes) -> Monomial { + VecZ idx; + std::uniform_int_distribution dist(0, logical - 1); + std::vector chosen; + while (chosen.size() < modes) { + const size_t q = dist(rng); + bool dup = false; + for (const auto x : chosen) { + dup = dup || (x == q); + } + if (!dup) { + chosen.push_back(q); + idx.push_back(2 * q); + idx.push_back((2 * q) + 1); + } + } + return indices_to_bitset(idx); +} + +struct Tally { + size_t comparisons = 0; + size_t mismatches = 0; + size_t paired_out = 0; // samples that are fully paired (the unconditionally-kept branch) + size_t kept = 0; + size_t rejected = 0; +}; + +enum class Population : uint8_t { Uniform, Paired }; + +template +auto check_width(std::mt19937_64 &rng, size_t logical, Population pop, Tally &t) -> void { + const bool paired_pop = pop == Population::Paired; + for (int rep = 0; rep < 400; ++rep) { + const size_t kw = 1 + (rng() % 12); + const auto x = paired_pop ? draw_paired(rng, logical, 1 + (kw % 5)) : draw(rng, logical, kw); + + const size_t k = x.count(); + const size_t d = paired_mode_count(x); + + const auto ref = cutoff_sums(x, logical); + const auto got = cutoff_sums(k, d); + bool ok = got.xor_sum == ref.xor_sum && got.popcount_sum == ref.popcount_sum && got.or_sum == ref.or_sum + && is_paired(k, d) == is_paired(x); + t.comparisons += 4; + if (is_paired(k, d)) { + ++t.paired_out; + } + + // 0 rejects all but the paired branch, 12 keeps everything, the rest straddle the weights. + for (const unsigned int c : {0U, 1U, 4U, 6U, 12U}) { + const bool len = length_cutoff(k, d, c); + const bool sup = support_cutoff(k, d, c); + ok = ok && len == length_cutoff(x, c, logical) && sup == support_cutoff(x, c, logical); + // Assert the forwarding too, or a wrapper that dropped a term would hide behind itself. + ok = ok && len == length_keeps(k, d, c) && sup == support_keeps(k, d, c); + t.comparisons += 4; + t.kept += static_cast(len); + t.rejected += static_cast(!len); + } + + if (!ok) { + ++t.mismatches; + } + } +} + +} // namespace + +BOOST_AUTO_TEST_CASE(sparse_predicates_match_bitset_forms_across_widths) { + std::mt19937_64 rng(0x5A5E0DDULL); + Tally t; + + for (const auto pop : {Population::Uniform, Population::Paired}) { + check_width<32>(rng, 32, pop, t); // W = 64, one word, no active offset + check_width<32>(rng, 30, pop, t); // W = 64, active_bit_offset = 4 + check_width<48>(rng, 45, pop, t); // W = 96 -- not a multiple of 64 + check_width<64>(rng, 64, pop, t); // W = 128, exactly two words + check_width<128>(rng, 120, pop, t); + check_width<256>(rng, 250, pop, t); // the production shape + } + + BOOST_TEST(t.mismatches == 0U); + BOOST_TEST(t.comparisons > 40000U); // a loop that never ran would report zero mismatches too + BOOST_TEST(t.paired_out > 0U); // the unconditionally-kept branch must be reached + BOOST_TEST(t.kept > 0U); + BOOST_TEST(t.rejected > 0U); +} + +// The boundaries as literals: length compares k, support compares k - d (a paired mode spans two). +BOOST_AUTO_TEST_CASE(sparse_predicates_pin_their_boundaries) { + BOOST_TEST(is_paired(0U, 0U)); // the identity is fully paired by this definition + BOOST_TEST(is_paired(4U, 2U)); + BOOST_TEST(!is_paired(3U, 1U)); + + // Fully paired: kept at cutoff 0, which rejects everything else. + BOOST_TEST(length_keeps(4U, 2U, 0U)); + BOOST_TEST(support_keeps(4U, 2U, 0U)); + BOOST_TEST(!length_keeps(1U, 0U, 0U)); + BOOST_TEST(!support_keeps(1U, 0U, 0U)); + + // Unpaired k=5, d=1: length sees 5, support sees k - d = 4. + BOOST_TEST(!length_keeps(5U, 1U, 4U)); + BOOST_TEST(length_keeps(5U, 1U, 5U)); + BOOST_TEST(support_keeps(5U, 1U, 4U)); + BOOST_TEST(!support_keeps(5U, 1U, 3U)); +} diff --git a/cpp/tests/sparse_query_tests.cpp b/cpp/tests/sparse_query_tests.cpp new file mode 100644 index 00000000..a0a8f544 --- /dev/null +++ b/cpp/tests/sparse_query_tests.cpp @@ -0,0 +1,482 @@ +// 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. + +// Differential against the frozen dense oracle in dense_query_reference.h; cases chosen, not sampled. + +#include + +#include +#include +#include +#include +#include + +#include "monoprop/detail/evolution/layer_build/Common.h" +#include "monoprop/detail/evolution/layer_build/QueryCodec.h" +#include "monoprop/detail/evolution/layer_build/SparseQuery.h" + +#include "dense_query_reference.h" + +using namespace monoprop; +using namespace monoprop::detail; + +namespace { + +template +auto differential(const std::vector &pos, int phase) -> size_t { + using SQ = SparseQuery; + const size_t k = pos.size(); + + Monomial want; + for (const auto p : pos) { + want.set(p); + } + BOOST_REQUIRE_EQUAL(want.count(), k); // the caller must not hand us duplicates + VecZ dbuf; + test_ref::query_push(dbuf, want, phase); + Monomial dmono; + int dphase = 99; + test_ref::query_read(dbuf, 0, dmono, dphase); + BOOST_REQUIRE((dmono == want)); + BOOST_REQUIRE_EQUAL(dphase, phase); + + VecZ sbuf; + const size_t sw = SQ::push(sbuf, pos.data(), k, phase); + BOOST_REQUIRE_EQUAL(sbuf.size(), sw); + + BOOST_TEST(SQ::words_at(sbuf, 0) == sw); + BOOST_TEST(SQ::k_at(sbuf, 0) == k); + BOOST_TEST(SQ::phase_at(sbuf, 0) == phase); + + std::vector sout(k == 0 ? 1 : k); + const size_t snext = SQ::read_positions(sbuf, 0, sout.data()); + BOOST_TEST(snext == sw); + sout.resize(k); + BOOST_TEST(sout == pos, boost::test_tools::per_element()); + + Monomial sm; + int sp = 99; + (void)SQ::read_mono(sbuf, 0, sm, sp); + BOOST_TEST(sm.count() == k); + BOOST_TEST((sm == dmono)); + BOOST_TEST(sp == dphase); + + VecZ mbuf; + const size_t mw = SQ::push_mono(mbuf, want, phase); + BOOST_TEST(mw == sw); + BOOST_TEST(mbuf == sbuf, boost::test_tools::per_element()); + + return sw; +} + +auto strided(size_t k, size_t start, size_t step, size_t universe) -> std::vector { + std::vector v; + for (size_t j = 0; j < k; ++j) { + const size_t p = start + j * step; + if (p >= universe) { + break; + } + v.push_back(static_cast(p)); + } + return v; +} + +// Uniform draws are what can reach bitmap mode -- gap coding wins every regular pattern -- but see +// sparse_record_actually_exercises_its_bitmap_mode: at narrow widths they cannot reach it either. +auto scattered(size_t k, size_t universe, std::mt19937_64 &rng) -> std::vector { + std::vector pool(universe); + for (size_t j = 0; j < universe; ++j) { + pool[j] = static_cast(j); + } + std::shuffle(pool.begin(), pool.end(), rng); + pool.resize(std::min(k, universe)); + std::sort(pool.begin(), pool.end()); + return pool; +} + +// The reference d, written against the definition (mode m owns bits 2m and 2m+1), not the codec. +auto reference_pair_count(const std::vector &pos) -> size_t { + size_t d = 0; + for (size_t j = 0; j + 1 < pos.size(); ++j) { + if ((pos[j] % 2 == 0) && (pos[j + 1] == pos[j] + 1)) { + ++d; + } + } + return d; +} + +} // namespace + +// Flat names: boostAddTests.cmake strips the indentation encoding suite nesting, so a suite-wrapped +// case errors at setup having asserted nothing. + +BOOST_AUTO_TEST_CASE(sparse_record_agrees_with_the_dense_oracle_across_widths) { + for (const int phase : {-1, 0, 1}) { + for (const size_t k : {size_t{0}, size_t{1}, size_t{2}, size_t{5}, size_t{6}, size_t{7}, size_t{15}}) { + differential<32>(strided(k, 0, 2, 64), phase); + differential<128>(strided(k, 3, 7, 256), phase); + differential<250>(strided(k, 11, 23, 500), phase); + differential<512>(strided(k, 1, 41, 1024), phase); + differential<1024>(strided(k, 5, 97, 2048), phase); + } + } +} + +BOOST_AUTO_TEST_CASE(sparse_record_handles_widths_that_are_not_whole_words) { + // Widths with no whole word: 12 modes (LiH) is kBits=24, so the bitmap payload is a partial word. + std::mt19937_64 rng(0xB17U); + for (const int phase : {-1, 0, 1}) { + for (const size_t k : {size_t{0}, size_t{1}, size_t{5}, size_t{9}, size_t{16}, size_t{24}}) { + differential<12>(strided(k, 0, 1, 24), phase); // kBits=24: no whole word at all + differential<12>(strided(k, 0, 2, 24), phase); + differential<12>(scattered(k, 24, rng), phase); + } + for (const size_t k : {size_t{1}, size_t{7}, size_t{20}, size_t{40}, size_t{70}}) { + differential<50>(strided(k, 0, 1, 100), phase); // kBits=100: one whole word plus 36 bits + differential<33>(strided(k, 0, 1, 66), phase); // kBits=66: one whole word plus 2 bits + differential<250>(strided(k, 0, 1, 500), phase); // kBits=500: seven words plus 52 bits + differential<50>(scattered(k, 100, rng), phase); + differential<33>(scattered(k, 66, rng), phase); + differential<250>(scattered(k, 500, rng), phase); + } + } + for (const size_t bits : {size_t{24}, size_t{100}, size_t{66}}) { + std::vector all(bits); + for (size_t j = 0; j < bits; ++j) { + all[j] = static_cast(j); + } + if (bits == 24) { + differential<12>(all, 1); + } + else if (bits == 100) { + differential<50>(all, 1); + } + else { + differential<33>(all, 1); + } + } +} + +BOOST_AUTO_TEST_CASE(sparse_record_actually_exercises_its_bitmap_mode) { + // Two generators because uniform draws CANNOT reach bitmap at kBits=24: bitmap is one word there, so + // gap must need two (19 + (k-1)*gw > 64), which only a dense run plus one far outlier forces. + const auto count_bitmap = [](auto tag, size_t universe) { + using SQ = SparseQuery; + size_t used = 0; + size_t bad = 0; + const auto tally = [&](const std::vector &pos) { + if (pos.size() < 2 || pos.size() > SQ::kMaxPositions) { + return; + } + VecZ buf; + (void)SQ::push(buf, pos.data(), pos.size(), 1); + if (SQ::header_at(buf, 0).mode != SQ::kModeBitmap) { + return; + } + ++used; + std::vector back(pos.size()); + SQ::read_positions(buf, 0, back.data()); + bad += static_cast(back != pos || SQ::k_at(buf, 0) != pos.size()); + }; + std::mt19937_64 rng(0xB1747U ^ universe); + for (size_t trial = 0; trial < 600; ++trial) { + tally(scattered(1 + (rng() % universe), universe, rng)); + } + for (size_t run = 1; run < universe; ++run) { + for (size_t outlier = run; outlier < universe; ++outlier) { + std::vector pos; + pos.reserve(run + 1); + for (size_t j = 0; j < run; ++j) { + pos.push_back(static_cast(j)); + } + pos.push_back(static_cast(outlier)); + tally(pos); + } + } + return std::pair{used, bad}; + }; + const auto narrow = count_bitmap(std::integral_constant{}, 24); + const auto partial = count_bitmap(std::integral_constant{}, 500); + const auto bucket = count_bitmap(std::integral_constant{}, 256); + BOOST_TEST(narrow.first > 0U); + BOOST_TEST(partial.first > 0U); + BOOST_TEST(bucket.first > 0U); + BOOST_TEST(narrow.second == 0U); + BOOST_TEST(partial.second == 0U); + BOOST_TEST(bucket.second == 0U); +} + +BOOST_AUTO_TEST_CASE(sparse_record_survives_the_six_bit_k_escape) { + for (const size_t k : {size_t{62}, size_t{63}, size_t{64}, size_t{200}}) { + const auto pos = strided(k, 0, 3, 2048); + BOOST_REQUIRE_EQUAL(pos.size(), k); + differential<1024>(pos, 1); + } +} + +BOOST_AUTO_TEST_CASE(sparse_record_bounds_the_fully_paired_term) { + // Every bit set: 514 words in FIXED at NumModes=1024, one word in BITMAP. Hence the argmin. + std::vector all(2048); + for (size_t j = 0; j < all.size(); ++j) { + all[j] = static_cast(j); + } + const size_t sw = differential<1024>(all, 1); + BOOST_TEST(sw == 1U); +} + +BOOST_AUTO_TEST_CASE(sparse_record_never_exceeds_the_dense_words_it_replaces) { + const auto check = [](auto tag, size_t universe, size_t dense_words) { + using SQ = SparseQuery; + std::mt19937_64 rng(0xC0FFEE ^ universe); + for (size_t k = 0; k <= universe; k += std::max(1, universe / 37)) { + const auto pos = scattered(k, universe, rng); + VecZ buf; + const size_t w = SQ::push(buf, pos.data(), pos.size(), 1); + BOOST_TEST(w <= dense_words + 1, + "k=" << k << " at U=" << universe << " took " << w << " words vs dense " << dense_words); + } + }; + check(std::integral_constant{}, 64, 2); + check(std::integral_constant{}, 256, 5); + check(std::integral_constant{}, 500, 9); +} + +BOOST_AUTO_TEST_CASE(sparse_record_walks_a_multi_query_buffer_exactly) { + // Mixed width and mixed mode, which is the case a hardcoded stride gets wrong. + using SQ = SparseQuery<128>; + using QC = QueryCodec<128>; + const QueryLayout layout{/*fused=*/false}; + VecZ buf; + std::vector offs; + size_t off = 0; + const std::vector> terms = { + strided(3, 0, 1, 256), // consecutive -> gap width 0 + strided(6, 10, 40, 256), // wide gaps -> gap width near the raw position width + {}, // empty + strided(40, 0, 6, 256), // wide enough that bitmap becomes competitive + strided(1, 255, 1, 256), // single position at the very top + strided(20, 7, 2, 256), // uniform stride 2 + }; + for (const auto &t : terms) { + offs.push_back(off); + off += SQ::push(buf, t.data(), t.size(), 1); + } + BOOST_TEST(QC::count_queries(buf, layout) == terms.size()); + + off = 0; + for (size_t i = 0; i < terms.size(); ++i) { + BOOST_TEST(off == offs[i]); + BOOST_TEST(SQ::k_at(buf, off) == terms[i].size()); + std::vector out(terms[i].size() + 1); + (void)SQ::read_positions(buf, off, out.data()); + out.resize(terms[i].size()); + BOOST_TEST(out == terms[i], boost::test_tools::per_element()); + off = QC::next_off(buf, layout, off); + } + BOOST_TEST(off == buf.size()); +} + +BOOST_AUTO_TEST_CASE(sparse_record_picks_the_smallest_of_its_three_modes) { + // Load-bearing: gap ALONE is 8.56 B/term against fixed lanes' 8.00 on uniform draws at 250 modes. + using SQ = SparseQuery<128>; + std::mt19937_64 rng(12345); + for (size_t trial = 0; trial < 400; ++trial) { + const size_t k = rng() % 60; + const auto pos = scattered(k, 256, rng); + VecZ buf; + const size_t w = SQ::push(buf, pos.data(), pos.size(), 1); + const size_t gwid = SQ::gap_width(pos.data(), pos.size()); + const size_t best = std::min({SQ::words_of(SQ::fixed_bits(pos.size())), + SQ::words_of(SQ::gap_bits(pos.size(), gwid)), + SQ::words_of(SQ::bitmap_bits(pos.size()))}); + BOOST_TEST(w == best, "k=" << k << " wrote " << w << " words, best was " << best); + } +} + +BOOST_AUTO_TEST_CASE(sparse_record_position_width_is_the_compile_time_bucket) { + BOOST_TEST(SparseQuery<32>::kPosBits == 6U); // U=64 + BOOST_TEST(SparseQuery<128>::kPosBits == 8U); // U=256, both lattice models + BOOST_TEST(SparseQuery<250>::kPosBits == 9U); // U=500 + BOOST_TEST(SparseQuery<512>::kPosBits == 10U); // U=1024 + BOOST_TEST(SparseQuery<1024>::kPosBits == 11U); // U=2048 +} + +BOOST_AUTO_TEST_CASE(sparse_record_carries_the_extreme_bit_positions) { + // MSb0 ordering puts logical index 0 at the TOP, so bit 2N-1 is the common case, not a rare one. + differential<32>({0}, 1); + differential<32>({63}, 1); + differential<32>({0, 63}, -1); + differential<128>({0, 255}, 1); + differential<250>({0, 499}, 1); + differential<1024>({0, 2047}, -1); + differential<12>({0, 23}, 1); // the non-word-multiple width, at both extremes +} + +BOOST_AUTO_TEST_CASE(sparse_record_encoding_is_deterministic) { + std::mt19937_64 rng(0xDE7ULL); + for (size_t trial = 0; trial < 200; ++trial) { + const auto pos = scattered(rng() % 40, 256, rng); + VecZ a; + VecZ b; + const size_t wa = SparseQuery<128>::push(a, pos.data(), pos.size(), 1); + const size_t wb = SparseQuery<128>::push(b, pos.data(), pos.size(), 1); + BOOST_TEST(wa == wb); + BOOST_TEST(a == b, boost::test_tools::per_element()); + } +} + +BOOST_AUTO_TEST_CASE(sparse_record_pair_count_recomputes_d_from_positions) { + // d is recomputed from positions, so it must agree with the definition: mode m owns bits 2m, 2m+1. + std::mt19937_64 rng(0xD1D1ULL); + for (size_t trial = 0; trial < 100; ++trial) { + const auto pos = scattered(rng() % 30, 256, rng); + BOOST_TEST(SparseQuery<128>::pair_count(pos.data(), pos.size()) == reference_pair_count(pos)); + } + const std::vector straddle{1, 2, 5, 6}; + BOOST_TEST(SparseQuery<128>::pair_count(straddle.data(), straddle.size()) == 0U); + const std::vector real{2, 3, 6, 7}; + BOOST_TEST(SparseQuery<128>::pair_count(real.data(), real.size()) == 2U); + std::vector paired; + for (uint16_t m = 0; m < 16; ++m) { + paired.push_back(static_cast(2 * m)); + paired.push_back(static_cast(2 * m + 1)); + } + BOOST_TEST(SparseQuery<128>::pair_count(paired.data(), paired.size()) == paired.size() / 2); +} + +BOOST_AUTO_TEST_CASE(sparse_record_fused_stream_interleaves_values_and_stays_walkable) { + using QC = QueryCodec<128>; + std::mt19937_64 rng(0xF5EDULL); + std::vector> terms; + std::vector vals; + VecZ plain; + for (size_t i = 0; i < 24; ++i) { + terms.push_back(scattered(rng() % 45, 256, rng)); + vals.push_back(static_cast(i) * 0.5 - 3.25); + (void)QC::push( + plain, + [&] { + Monomial<128> m; + for (const auto p : terms.back()) { + m.set(p); + } + return m; + }(), + 1); + } + VecZ fused; + QC::build_fused(plain, vals, fused); + + const QueryLayout layout{.fused = true}; + BOOST_TEST(QC::count_queries(fused, layout) == terms.size()); + size_t off = 0; + for (size_t i = 0; i < terms.size(); ++i) { + BOOST_TEST(QC::k_at(fused, off) == terms[i].size()); + BOOST_TEST(QC::value_at(fused, layout, off) == vals[i]); + std::vector out(terms[i].size() + 1); + int phase = 0; + const size_t next = QC::read_positions(fused, layout, off, out.data(), phase); + out.resize(terms[i].size()); + BOOST_TEST(out == terms[i], boost::test_tools::per_element()); + off = QC::next_off(fused, layout, off); + BOOST_TEST(next == off); + } + BOOST_TEST(off == fused.size()); +} + +// A separate case because `-0.0 == 0.0` is TRUE, so bit-exactness needs its own assertion. +BOOST_AUTO_TEST_CASE(sparse_record_fused_value_channel_is_bit_exact_and_reusable) { + using QC = QueryCodec<128>; + const QueryLayout layout{.fused = true}; + + auto push_terms = [](VecZ &buf, const std::vector> &terms) { + for (const auto &t : terms) { + Monomial<128> m; + for (const auto p : t) { + m.set(p); + } + (void)QC::push(buf, m, 1); + } + }; + + // 1. BIT-EXACTNESS via memcmp, so -0.0 stays distinguished from 0.0. Widths differ per term. + const std::vector values = { + 0.0, + -0.0, + 1.0, + -1.0, + 3.141592653589793, + -2.718281828459045e-300, // near-denormal magnitude + std::numeric_limits::min(), // smallest normal + std::numeric_limits::denorm_min(), + }; + const std::vector> terms = { + {3}, + {0, 255}, + {1, 2, 3, 4, 5, 6, 7}, + {9, 40}, + {2, 3}, + {5, 60, 61, 200}, + {17}, + {0, 1, 2}, + }; + BOOST_REQUIRE_EQUAL(terms.size(), values.size()); + + VecZ plain; + push_terms(plain, terms); + VecZ fused; + QC::build_fused(plain, values, fused); + + size_t off = 0; + for (size_t i = 0; i < values.size(); ++i) { + const double v_out = QC::value_at(fused, layout, off); + BOOST_CHECK(std::memcmp(&v_out, &values[i], sizeof(double)) == 0); + off = QC::next_off(fused, layout, off); + } + BOOST_TEST(off == fused.size()); + + // 2. BUFFER REUSE: size must be exact, or a shorter gate reads the previous gate's trailing words. + VecZ plain_big; + std::vector vbig; + std::vector> big; + for (size_t r = 0; r < 32; ++r) { + big.push_back({static_cast(r), static_cast(r + 60)}); + vbig.push_back(static_cast(r) * 1.5 - 7.0); + } + push_terms(plain_big, big); + VecZ out; + QC::build_fused(plain_big, vbig, out); + const size_t cap_after_big = out.capacity(); + + VecZ plain_small; + const std::vector vsmall = {42.0, -42.0, 0.25}; + const std::vector> small = {{1}, {2, 3}, {4, 5, 6}}; + push_terms(plain_small, small); + QC::build_fused(plain_small, vsmall, out); + BOOST_TEST(QC::count_queries(out, layout) == vsmall.size()); + BOOST_CHECK_GE(out.capacity(), cap_after_big); + off = 0; + for (size_t i = 0; i < vsmall.size(); ++i) { + const double v_out = QC::value_at(out, layout, off); + BOOST_CHECK(std::memcmp(&v_out, &vsmall[i], sizeof(double)) == 0); + off = QC::next_off(out, layout, off); + } + BOOST_TEST(off == out.size()); + + // 3. EMPTY INPUT: the self slot is cleared before the exchange, into a buffer holding stale words. + VecZ empty; + VecZ dirty{1, 2, 3}; + QC::build_fused(empty, {}, dirty); + BOOST_TEST(dirty.empty()); +} diff --git a/cpp/tests/sparse_resolve_tests.cpp b/cpp/tests/sparse_resolve_tests.cpp new file mode 100644 index 00000000..37b70873 --- /dev/null +++ b/cpp/tests/sparse_resolve_tests.cpp @@ -0,0 +1,370 @@ +// 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 position-form resolve path, differentially against the queries the caller built and the dense +// Monomial-keyed insert path, neither of which is the code under test. + +#include + +#include +#include +#include +#include +#include +#include + +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/evolution/layer_build/QueryCodec.h" +#include "monoprop/detail/evolution/layer_build/Resolve.h" +#include "monoprop/detail/operator/MPOperator.h" +#include "monoprop/detail/operator/OperatorIndex.h" + +using namespace monoprop; + +namespace { + +template +auto random_monomial(std::mt19937_64 &rng, size_t k) -> Monomial { + Monomial m; + std::uniform_int_distribution bit(0, Monomial::size() - 1); + size_t placed = 0; + while (placed < k) { + const size_t b = bit(rng); + if (!m.test(b)) { + m.set(b); + ++placed; + } + } + return m; +} + +// Fully paired terms are the only source of wide records: 94 in 20.9M in production, so drawn here. +template +auto random_paired_monomial(std::mt19937_64 &rng, size_t d) -> Monomial { + Monomial m; + std::uniform_int_distribution mode(0, NumModes - 1); + size_t placed = 0; + while (placed < d) { + const size_t mo = mode(rng); + if (!m.test(2 * mo)) { + m.set(2 * mo); + m.set((2 * mo) + 1); + ++placed; + } + } + return m; +} + +// 0 and 1 for the degenerate records, up to 20 for multi-word ones, 14 for the overflow spill. +const std::vector kPopcounts = {0, 1, 2, 4, 5, 6, 7, 8, 11, 12, 14, 20}; + +template +auto make_op(const std::vector> &terms) -> detail::MPOperator { + detail::MPOperator op; + op.basis = Basis::Majorana; + if (terms.empty()) { + return op; + } + detail::insert_absent_terms( + op, + terms.size(), + [&](size_t k) -> const Monomial & { return terms[k]; }, + [&](size_t k, size_t base) { assign_row(*op.store, base + k, terms[k]); }); + return op; +} + +template +auto draw_distinct(std::mt19937_64 &rng, size_t n) -> std::vector> { + std::vector> out; + std::set> seen; + std::uniform_int_distribution pick(0, kPopcounts.size() - 1); + while (out.size() < n) { + const size_t k = kPopcounts[pick(rng)]; + const auto m = + ((rng() & 1U) != 0U) ? random_paired_monomial(rng, k / 2) : random_monomial(rng, k); + std::vector key; + key.reserve(Monomial::num_words()); + for (size_t w = 0; w < Monomial::num_words(); ++w) { + key.push_back(m.word(w)); + } + if (seen.insert(key).second) { + out.push_back(m); + } + } + return out; +} + +template +auto serialize(const std::vector>> &queries, bool fused) -> std::vector { + std::vector incoming(queries.size()); + for (size_t s = 0; s < queries.size(); ++s) { + for (size_t q = 0; q < queries[s].size(); ++q) { + const int phase = ((q % 2) == 0) ? 1 : -1; + detail::QueryCodec::push(incoming[s], queries[s][q], phase); + if (fused) { + detail::QueryCodec::push_value(incoming[s], 0.5 + static_cast(q)); + } + } + } + return incoming; +} + +template +auto check_probe_matches_the_queries(std::mt19937_64 &rng, size_t n_seed, size_t n_query, size_t rank_count, bool fused) + -> void { + const auto seed_terms = draw_distinct(rng, n_seed); + const auto fresh_terms = draw_distinct(rng, n_query); + + // Hits matter even though the production hit rate is ~0: only they exercise the confirm. + std::vector>> queries(rank_count); + std::set> queried; + size_t hits_planned = 0; + size_t misses_planned = 0; + for (size_t i = 0; i < n_query; ++i) { + const bool want_hit = (i % 3) == 0 && !seed_terms.empty(); + const auto m = want_hit ? seed_terms[i % seed_terms.size()] : fresh_terms[i]; + std::vector key; + for (size_t w = 0; w < Monomial::num_words(); ++w) { + key.push_back(m.word(w)); + } + // A repeat would violate bulk_insert's precondition; the engine gets distinctness from ^G. + if (!queried.insert(key).second) { + continue; + } + (want_hit ? hits_planned : misses_planned) += 1; + queries[i % rank_count].push_back(m); + } + BOOST_REQUIRE(hits_planned > 0); + BOOST_REQUIRE(misses_planned > 0); + + std::vector> expect_mono; + std::vector expect_phase; + std::vector expect_sender; + for (size_t s = 0; s < rank_count; ++s) { + for (size_t q = 0; q < queries[s].size(); ++q) { + expect_mono.push_back(queries[s][q]); + expect_phase.push_back(((q % 2) == 0) ? 1 : -1); + expect_sender.push_back(s); + } + } + + const auto incoming = serialize(queries, fused); + const detail::QueryLayout layout{fused}; + + auto op = make_op(seed_terms); + const auto pr = detail::probe_incoming_queries(incoming, op, rank_count, layout); + + BOOST_REQUIRE_EQUAL(pr.nq_total, expect_mono.size()); + BOOST_REQUIRE(pr.nq_total > 0); + BOOST_REQUIRE_EQUAL(pr.pos_off.size(), pr.nq_total); + + std::set> seeded; + for (const auto &m : seed_terms) { + std::vector key; + for (size_t w = 0; w < Monomial::num_words(); ++w) { + key.push_back(m.word(w)); + } + seeded.insert(key); + } + + size_t hits_seen = 0; + size_t wide_seen = 0; + std::vector> expected_misses; + for (size_t g = 0; g < pr.nq_total; ++g) { + const Monomial &want = expect_mono[g]; + BOOST_TEST((pr.mono_at(g) == want)); + BOOST_TEST(pr.k_of[g] == want.count()); + BOOST_TEST(pr.phase_of[g] == expect_phase[g]); + BOOST_TEST(pr.sender_of[g] == expect_sender[g]); + BOOST_TEST(pr.is_paired_at(g) == monoprop::is_paired(want)); + + std::vector key; + for (size_t w = 0; w < Monomial::num_words(); ++w) { + key.push_back(want.word(w)); + } + const bool want_hit = seeded.count(key) != 0; + BOOST_TEST((pr.idx_of[g] < pr.base) == want_hit); + if (want_hit) { + ++hits_seen; + } + else { + expected_misses.push_back(want); + } + VecZ scratch; + if (detail::QueryCodec::push(scratch, want, expect_phase[g]) > 1U) { + ++wide_seen; + } + } + // Vacuous-pass guards: no hit means the confirm never ran, no wide record means the cursor didn't. + BOOST_TEST(hits_seen > 0); + BOOST_TEST(wide_seen > 0); + + BOOST_REQUIRE_EQUAL(pr.miss_g.size(), expected_misses.size()); + for (size_t j = 0; j < pr.miss_g.size(); ++j) { + BOOST_TEST((expect_mono[pr.miss_g[j]] == expected_misses[j])); + BOOST_TEST(pr.idx_of[pr.miss_g[j]] == pr.base + j); + } + + detail::insert_incoming_misses(op, pr); + + // The second implementation: the dense Monomial-keyed path, sharing no code with set_positions. + auto ref = make_op(seed_terms); + detail::insert_absent_terms( + ref, + expected_misses.size(), + [&](size_t j) -> const Monomial & { return expected_misses[j]; }, + [&](size_t j, size_t base) { assign_row(*ref.store, base + j, expected_misses[j]); }); + + BOOST_REQUIRE_EQUAL(op.store->size(), ref.store->size()); + BOOST_TEST(op.store->size() > pr.base); + size_t overflow_seen = 0; + for (size_t i = 0; i < ref.store->size(); ++i) { + BOOST_TEST((op.store->row(i) == ref.store->row(i))); + BOOST_TEST(op.store->popcount(i) == ref.store->popcount(i)); + if (!ref.store->row_positions(i).inlined()) { + ++overflow_seen; + } + } + BOOST_TEST(overflow_seen > 0); + + // The index, not just the rows: a wrong hash leaves the row correct and unfindable. + for (size_t i = 0; i < ref.store->size(); ++i) { + const auto key = ref.store->row(i); + const auto in_op = op.store->find(key); + const auto in_ref = ref.store->find(key); + BOOST_REQUIRE(in_ref.has_value()); + BOOST_REQUIRE(in_op.has_value()); + BOOST_TEST(*in_op == *in_ref); + BOOST_TEST(*in_ref == i); + } +} + +} // namespace + +/* ── The check, across both position widths and both buffer layouts ── */ + +BOOST_AUTO_TEST_CASE(sparse_resolve_probe_matches_narrow_positions) { + std::mt19937_64 rng(20260814); + static_assert(sizeof(detail::OperatorIndex<32>::PosT) == 1, "this case exists to cover the narrowing decode"); + check_probe_matches_the_queries<32>(rng, /*n_seed=*/40, /*n_query=*/90, /*rank_count=*/3, /*fused=*/false); +} + +BOOST_AUTO_TEST_CASE(sparse_resolve_probe_matches_wide_positions) { + std::mt19937_64 rng(20260815); + static_assert(sizeof(detail::OperatorIndex<250>::PosT) == 2, "this case exists to cover the wide store"); + check_probe_matches_the_queries<250>(rng, /*n_seed=*/60, /*n_query=*/140, /*rank_count=*/4, /*fused=*/false); +} + +BOOST_AUTO_TEST_CASE(sparse_resolve_probe_matches_fused_layout) { + std::mt19937_64 rng(20260816); + check_probe_matches_the_queries<250>(rng, /*n_seed=*/50, /*n_query=*/120, /*rank_count=*/2, /*fused=*/true); +} + +BOOST_AUTO_TEST_CASE(sparse_resolve_probe_matches_single_sender) { + std::mt19937_64 rng(20260817); + check_probe_matches_the_queries<32>(rng, /*n_seed=*/25, /*n_query=*/60, /*rank_count=*/1, /*fused=*/true); +} + +/* ── The pieces, pinned individually ──────────────────────────────────────── */ + +BOOST_AUTO_TEST_CASE(sparse_resolve_set_positions_matches_set) { + constexpr size_t kN = 250; + constexpr size_t kInlineWidth = 11; + std::mt19937_64 rng(20260818); + const auto terms = draw_distinct(rng, 200); + + detail::OperatorIndex from_mono(kInlineWidth); + detail::OperatorIndex from_pos(kInlineWidth); + from_mono.grow_rows_geometric(terms.size()); + from_pos.grow_rows_geometric(terms.size()); + + size_t spilled = 0; + for (size_t i = 0; i < terms.size(); ++i) { + from_mono.set(i, terms[i]); + std::vector::PosT> pos; + for (size_t b = terms[i].find_first(); b < terms[i].size(); b = terms[i].find_next(b)) { + pos.push_back(static_cast::PosT>(b)); + } + from_pos.set_positions(i, pos.data(), pos.size()); + if (pos.size() > kInlineWidth) { + ++spilled; + } + } + BOOST_TEST(spilled > 0); + for (size_t i = 0; i < terms.size(); ++i) { + BOOST_TEST((from_pos.row(i) == from_mono.row(i))); + BOOST_TEST((from_pos.row(i) == terms[i])); + BOOST_TEST(from_pos.popcount(i) == from_mono.popcount(i)); + BOOST_TEST(from_pos.row_positions(i).inlined() == from_mono.row_positions(i).inlined()); + } + BOOST_TEST(from_pos.overflow_size() == from_mono.overflow_size()); +} + +BOOST_AUTO_TEST_CASE(sparse_resolve_finds_dense_inserted_keys) { + // The hash identity, isolated: fold_hash_positions differing from fold_hash misses, and legally. + constexpr size_t kN = 250; + std::mt19937_64 rng(20260819); + const auto terms = draw_distinct(rng, 300); + auto op = make_op(terms); + + std::vector::PosT> flat; + std::vector off; + std::vector kk; + for (const auto &m : terms) { + off.push_back(flat.size()); + size_t k = 0; + for (size_t b = m.find_first(); b < m.size(); b = m.find_next(b)) { + flat.push_back(static_cast::PosT>(b)); + ++k; + } + kk.push_back(static_cast(k)); + } + std::vector out(terms.size(), 0); + std::vector hashes(terms.size(), 0); + op.store->find_batch_positions(flat.data(), off.data(), kk.data(), terms.size(), out.data(), hashes.data()); + + std::vector out_dense(terms.size(), 0); + op.store->find_batch(terms.data(), terms.size(), out_dense.data()); + for (size_t i = 0; i < terms.size(); ++i) { + BOOST_REQUIRE(out[i] != detail::OperatorIndex::kNotFound); + BOOST_TEST(out[i] == i); + BOOST_TEST(out[i] == out_dense[i]); + BOOST_TEST(hashes[i] == detail::OperatorIndex::fold_hash_positions(flat.data() + off[i], kk[i])); + } + + const auto absent = draw_distinct(rng, 50); + std::vector::PosT> aflat; + std::vector aoff; + std::vector akk; + for (const auto &m : absent) { + aoff.push_back(aflat.size()); + size_t k = 0; + for (size_t b = m.find_first(); b < m.size(); b = m.find_next(b)) { + aflat.push_back(static_cast::PosT>(b)); + ++k; + } + akk.push_back(static_cast(k)); + } + std::vector aout(absent.size(), 0); + op.store->find_batch_positions(aflat.data(), aoff.data(), akk.data(), absent.size(), aout.data(), nullptr); + size_t genuinely_absent = 0; + for (size_t i = 0; i < absent.size(); ++i) { + // draw_distinct may re-draw a seeded term; only genuinely absent ones are evidence. + if (!op.store->find(absent[i]).has_value()) { + BOOST_TEST(aout[i] == detail::OperatorIndex::kNotFound); + ++genuinely_absent; + } + } + BOOST_TEST(genuinely_absent > 0); +} From 058bcd0012458a6eac18de9afa0a3143a0f98029 Mon Sep 17 00:00:00 2001 From: Aaron Miller Date: Mon, 24 Aug 2026 22:01:06 +0100 Subject: [PATCH 2/5] =?UTF-8?q?perf(evolution):=20=E2=9A=A1=20stage=20self?= =?UTF-8?q?-owned=20queries=20as=20positions?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The self leg encoded a query record and decoded it straight back on the same rank: Scan pushed through QueryCodec for every destination including my_rank, and resolve_range_ inflated it again into self_pos_flat_ before calling find_batch_positions -- whose signature is already the flat-positions form the encoder had just destroyed. At one partition that round trip is 100% of queries against exchange_s = 0.0000. SelfQueryStage carries them as positions instead, and merge_partner_positions supplies the positions, k, the cutoff digest's d and the overlap in the one merge that replaces the walk. Four measured shapes were rejected on the way and are recorded at their sites so they are not tried again: a lambda in the merge (15,279,191 non-inlined calls), thread_local for the position buffer (__tls_get_addr, 54.7M), a byte loop for row_eq_positions (glibc's AVX2 memcmp wins at ~5 bytes), and walking the dense partner instead of merging (+229.1M, find_next is a serial dependence chain). Layout 1x1 against the query-wire record, 10 reps interleaved, 10/10 p=.00: build_graph[hubbard] 0.87x, propagate[hubbard] 0.90x, build_graph[pauli] 0.79x, propagate[pauli] 0.76x. energy/gradient 1.00x as the null control, peak RSS flat. Callgrind: 9,720.99M -> 9,265.72M instructions, floor 0.052%. Co-Authored-By: Claude Opus 5 --- cpp/monoprop/algebra/AlgebraCommon.h | 23 ++- cpp/monoprop/core/CMakeLists.txt | 1 + .../evolution/layer_build/CMakeLists.txt | 3 + .../detail/evolution/layer_build/Engine.h | 62 +++---- .../evolution/layer_build/PartnerMerge.h | 172 ++++++++++++++++++ .../detail/evolution/layer_build/QueryCodec.h | 7 + .../detail/evolution/layer_build/Scan.h | 135 ++++++++++---- cpp/monoprop/detail/operator/OperatorIndex.h | 3 + cpp/tests/mpi_utils_tests.cpp | 33 +++- cpp/tests/partner_merge_tests.cpp | 168 +++++++++++++++++ 10 files changed, 524 insertions(+), 83 deletions(-) create mode 100644 cpp/monoprop/detail/evolution/layer_build/PartnerMerge.h create mode 100644 cpp/tests/partner_merge_tests.cpp diff --git a/cpp/monoprop/algebra/AlgebraCommon.h b/cpp/monoprop/algebra/AlgebraCommon.h index 87da0c16..9b6f8eb5 100644 --- a/cpp/monoprop/algebra/AlgebraCommon.h +++ b/cpp/monoprop/algebra/AlgebraCommon.h @@ -278,20 +278,29 @@ class CutoffEvaluator { return cutoff_fn_(mono); } - // Same decision without cutoff_sums: the caller already knows k, so only d is computed. No - // popcount early-out, deliberately -- the digest is cheaper than the branch. nullopt if opaque. - auto passes_from_dense(const Monomial &mono, size_t k) const -> std::optional { - // paired_mode_count has no active_mask, so it agrees with cutoff_sums(mono, L) only above it. - assert(mono.find_first() >= active_bit_offset_() && "monomial has a set bit below its active offset"); + // The decision from the (k, d) digest alone. The emit site gets both out of the partner merge, so + // nothing here reads a bitset. Same precondition as cutoff_sums(k, d): d must have been folded + // without an active mask, which holds for a well-formed monomial. nullopt if opaque. + auto passes_from_digest(size_t k, size_t d) const -> std::optional { if (length_cutoff_ != nullptr) { - return length_keeps(k, paired_mode_count(mono), length_cutoff_->cutoff); + return length_keeps(k, d, length_cutoff_->cutoff); } if (support_cutoff_ != nullptr) { - return support_keeps(k, paired_mode_count(mono), support_cutoff_->cutoff); + return support_keeps(k, d, support_cutoff_->cutoff); } return std::nullopt; } + // The same decision when only the dense form is at hand, so d must be folded out of it. + auto passes_from_dense(const Monomial &mono, size_t k) const -> std::optional { + // paired_mode_count has no active_mask, so it agrees with cutoff_sums(mono, L) only above it. + assert(mono.find_first() >= active_bit_offset_() && "monomial has a set bit below its active offset"); + if (length_cutoff_ == nullptr && support_cutoff_ == nullptr) { + return std::nullopt; + } + return passes_from_digest(k, paired_mode_count(mono)); + } + // Upper bound on the set bits (physical slots) a surviving term can carry, so the store can size // its packed inline rows. A length cutoff counts set bits directly; a support cutoff counts // modes/qubits, each spanning two slots, hence the x2. diff --git a/cpp/monoprop/core/CMakeLists.txt b/cpp/monoprop/core/CMakeLists.txt index d9faa730..6fb5707c 100644 --- a/cpp/monoprop/core/CMakeLists.txt +++ b/cpp/monoprop/core/CMakeLists.txt @@ -5,4 +5,5 @@ target_sources( TYPE HEADERS FILES "Monomial.h" + "SparseMonomial.h" ) diff --git a/cpp/monoprop/detail/evolution/layer_build/CMakeLists.txt b/cpp/monoprop/detail/evolution/layer_build/CMakeLists.txt index eb4c9853..f0c81013 100644 --- a/cpp/monoprop/detail/evolution/layer_build/CMakeLists.txt +++ b/cpp/monoprop/detail/evolution/layer_build/CMakeLists.txt @@ -7,6 +7,9 @@ target_sources( "Common.h" "Engine.h" "FusedApply.h" + "PartnerMerge.h" + "QueryCodec.h" "Resolve.h" "Scan.h" + "SparseQuery.h" ) diff --git a/cpp/monoprop/detail/evolution/layer_build/Engine.h b/cpp/monoprop/detail/evolution/layer_build/Engine.h index 88983362..0f03f4af 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Engine.h +++ b/cpp/monoprop/detail/evolution/layer_build/Engine.h @@ -29,6 +29,7 @@ #include "monoprop/algebra/Algebra.h" #include "monoprop/detail/evolution/CutoffContext.h" #include "monoprop/detail/evolution/layer_build/Common.h" +#include "monoprop/detail/evolution/layer_build/PartnerMerge.h" #include "monoprop/detail/evolution/layer_build/QueryCodec.h" #include "monoprop/detail/evolution/layer_build/Resolve.h" #include "monoprop/detail/evolution/layer_build/Scan.h" @@ -342,8 +343,9 @@ struct LayerBuildEngine { std::vector deferred_self_misses; // Deferred-miss positions, concatenated in miss order; parallel to deferred_self_misses. std::vector deferred_pos_flat_; - // Per-batch decode scratch for resolve_range_; a member so one allocation serves every batch. - std::vector self_pos_flat_; + // This pass's self-owned queries as positions, straight from the scan: never encoded, so the resolve + // below has nothing to decode. Parallel to src_idx_r[my_rank]. + SelfQueryStage self_stage_; // Scan-captured v_src per query (ContractSink only via Sink::wants_values; empty for GraphSink). std::vector> src_val_r; // Fused query+value send scratch (ContractSink, R>1): shared by a gate's two exchange passes. @@ -371,17 +373,16 @@ struct LayerBuildEngine { // Resolve this rank's own query stream inline, then clear it so the alltoallv never sends to self. auto resolve_self_queries(bool is_leader_pass) -> void { - VecZ &lq = queries_r[my_rank]; std::vector &ls = src_idx_r[my_rank]; std::vector *lv = nullptr; if constexpr (Sink::wants_values) { lv = &src_val_r[my_rank]; } - // One source per query, pushed by the scan, so the count needs neither a walk nor a division. - assert(ls.size() == QueryCodec::count_queries(lq, sink.querier_layout()) - && "the self query buffer does not hold exactly one query per source"); - resolve_range_(lq, ls, lv, ls.size(), is_leader_pass); - lq.clear(); + // The scan routes a self-owned partner to the stage, never to the wire buffer. + assert(queries_r[my_rank].empty() && "a self-owned query was encoded instead of staged"); + assert(ls.size() == self_stage_.size() && "the self stage does not hold exactly one query per source"); + resolve_range_(ls, lv, is_leader_pass); + self_stage_.clear(); ls.clear(); if constexpr (Sink::wants_values) { src_val_r[my_rank].clear(); @@ -397,9 +398,11 @@ struct LayerBuildEngine { auto run_exchange(bool is_leader_pass, std::vector &&queries, std::vector> &&src_idx, - std::vector> &&src_val) -> void { + std::vector> &&src_val, + SelfQueryStage &&self_stage) -> void { queries_r = std::move(queries); src_idx_r = std::move(src_idx); + self_stage_ = std::move(self_stage); // src_val is empty unless Sink::wants_values, so the move is a no-op under GraphSink. src_val_r = std::move(src_val); if (!is_leader_pass && R > 1) { @@ -506,13 +509,11 @@ struct LayerBuildEngine { // Batched self-resolve over the index's group-prefetch find_batch; hits/misses are emitted to the sink // in query order. `lv` is the per-query v_src array parallel to `ls` (read only when Sink::wants_values). static constexpr size_t kResolveBatch = 64; - auto resolve_range_(VecZ &lq, - std::vector &ls, - [[maybe_unused]] std::vector *lv, - size_t hi, - bool is_leader_pass) -> void { + auto resolve_range_(std::vector &ls, [[maybe_unused]] std::vector *lv, bool is_leader_pass) + -> void { const size_t op_size = local_op.store->size(); - // pos_off/k_of index self_pos_flat_, rebuilt per batch but keeping its capacity; no dense keys. + // Gathered per batch because a matched follower is skipped; the offsets stay ABSOLUTE into the + // stage's pos_flat, so find_batch_positions reads it in place and nothing is copied. std::array pos_off; std::array k_of; std::array hashes; @@ -520,27 +521,18 @@ struct LayerBuildEngine { std::array srcs; std::array vals; std::array found; - using QC = QueryCodec; - const QueryLayout layout = sink.querier_layout(); - // A cursor, not an ordinal, and it must advance even when the query is skipped. - size_t off = 0; + const size_t hi = self_stage_.size(); size_t q = 0; while (q < hi) { size_t m = 0; - self_pos_flat_.clear(); for (; q < hi && m < kResolveBatch; ++q) { const size_t src = ls[q]; - const size_t this_off = off; if (!is_leader_pass && matched.is_marked(src)) { - off = QC::next_off(lq, layout, this_off); continue; // follower already matched by a leader → not an independent rotation } - const size_t k = QC::k_at(lq, this_off); - const size_t at = self_pos_flat_.size(); - self_pos_flat_.resize(at + k); // default-init grow: read_positions writes every element - off = QC::read_positions(lq, layout, this_off, self_pos_flat_.data() + at, phases[m]); - pos_off[m] = at; - k_of[m] = static_cast(k); + pos_off[m] = self_stage_.pos_off[q]; + k_of[m] = self_stage_.k_of[q]; + phases[m] = self_stage_.phase_of[q]; srcs[m] = src; if constexpr (Sink::wants_values) { vals[m] = (*lv)[q]; @@ -551,7 +543,7 @@ struct LayerBuildEngine { break; } // The hashes come back because a miss needs one at insert, folded from these same positions. - local_op.store->find_batch_positions(self_pos_flat_.data(), + local_op.store->find_batch_positions(self_stage_.pos_flat.data(), pos_off.data(), k_of.data(), m, @@ -571,9 +563,9 @@ struct LayerBuildEngine { sink.self_hit(srcs[j], found[j], phases[j], v_src); } else { - // self_pos_flat_ is cleared by the next batch, so copy now, into one gate-long buffer. + // The stage dies with this pass and the misses are flushed after both, so copy now. const size_t at = deferred_pos_flat_.size(); - const auto *const first = self_pos_flat_.data() + pos_off[j]; + const auto *const first = self_stage_.pos_flat.data() + pos_off[j]; deferred_pos_flat_.insert(deferred_pos_flat_.end(), first, first + k_of[j]); deferred_self_misses.push_back({at, k_of[j], hashes[j], srcs[j], phases[j], v_src}); } @@ -627,7 +619,7 @@ auto build_layer(MPOperator &local_op, } assert(fused_scale_coeffs == nullptr || (local_coeffs && &local_coeffs->get() == fused_scale_coeffs)); - FusedScanResult fused = [&] { + FusedScanResult fused = [&] { double *const sweep_ptr = fused_scale ? fused_scale_coeffs->data() : nullptr; return with_algebra(basis, [&]() { return fused_find_and_collect(local_op, @@ -669,11 +661,13 @@ auto build_layer(MPOperator &local_op, eng.run_exchange(/*is_leader_pass=*/true, std::move(fused.leader_queries), std::move(fused.leader_src), - std::move(fused.leader_val)); + std::move(fused.leader_val), + std::move(fused.leader_self)); eng.run_exchange(/*is_leader_pass=*/false, std::move(fused.follower_queries), std::move(fused.follower_src), - std::move(fused.follower_val)); + std::move(fused.follower_val), + std::move(fused.follower_self)); return eng.finish(std::move(cos_all), out_cos); }; diff --git a/cpp/monoprop/detail/evolution/layer_build/PartnerMerge.h b/cpp/monoprop/detail/evolution/layer_build/PartnerMerge.h new file mode 100644 index 00000000..712e4cc0 --- /dev/null +++ b/cpp/monoprop/detail/evolution/layer_build/PartnerMerge.h @@ -0,0 +1,172 @@ +// 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 + +// M⊕G as ascending positions. A slot in both M and G cancels (m_p m_p = 1), so the partner is the +// symmetric difference of two ascending position lists, and one merge yields its positions, `overlap` +// and `d` (modes carrying BOTH Majoranas) together -- the (k, d) digest the structural cutoff wants, +// with no second sweep over the dense form and no walk back out of it. + +#include +#include +#include +#include +#include + +#include "monoprop/TypeAliases.h" +#include "monoprop/detail/operator/OperatorIndex.h" + +namespace monoprop::detail { + +// Both inputs must be strictly ascending. The output is their symmetric difference, so it is bounded +// by the universe the positions are drawn from -- 2*NumModes here -- and ka + kb is only the bound that +// ignores cancellation. Returns the merged count. GenT is separate from PosT because the generator's +// positions are the wire's width, not the store's. +template +[[gnu::always_inline]] inline auto merge_partner_positions(const PosT *a, + size_t ka, + const GenT *b, + size_t kb, + PosT *out, + size_t &overlap_out, + size_t &d_out) noexcept -> size_t { + size_t i = 0; + size_t j = 0; + size_t n = 0; + size_t overlap = 0; + size_t d = 0; + // Seeded ODD, so the (prev % 2 == 0) test cannot fire on the first emit and the loops need no + // n != 0 guard; 1 is not a reachable `prev + 1` either, since prev would have to be 0 and even. + size_t prev = 1; + // Ascending output, so a doubly-occupied mode is an even position immediately followed by its + // successor -- the same count paired_mode_count folds out of the bitset. Written out three times + // rather than through a lambda: callgrind measured 15,279,191 CALLS to that lambda at 18 + // instructions each (275.1M, a third of this port's whole delta), because GCC declined to inline a + // closure capturing five locals by reference into three call sites. + while (i < ka && j < kb) { + const size_t pa = static_cast(a[i]); + const size_t pb = static_cast(b[j]); + if (pa == pb) { + ++overlap; + ++i; + ++j; + continue; + } + const size_t p = pa < pb ? pa : pb; + i += static_cast(pa < pb); + j += static_cast(pb < pa); + d += static_cast((prev % 2 == 0) && p == prev + 1); + out[n++] = static_cast(p); + prev = p; + } + for (; i < ka; ++i) { + const size_t p = static_cast(a[i]); + d += static_cast((prev % 2 == 0) && p == prev + 1); + out[n++] = static_cast(p); + prev = p; + } + for (; j < kb; ++j) { + const size_t p = static_cast(b[j]); + d += static_cast((prev % 2 == 0) && p == prev + 1); + out[n++] = static_cast(p); + prev = p; + } + overlap_out = overlap; + d_out = d; + return n; +} + +// Self-owned queries never reach a wire, so they are staged as positions rather than encoded records: +// OperatorIndex's find_batch_positions and set_positions both take exactly this shape, so the resolve +// path consumes the stage with no transformation and the codec is not on the self leg at all. +template +struct SelfQueryStage { + using PosT = typename OperatorIndex::PosT; + + // SIZED, not filled: the vectors carry the capacity and n_/pos_n_ carry the logical length, so a + // push writes rather than appends. Read them through size() and the data pointers only. + // + // DefaultInitVector, the allocator Resolve.h already uses for exactly this: a plain vector's resize + // VALUE-initialises, so sizing pos_flat ahead would memset every byte a push is about to overwrite + // -- trading the append cost for a per-gate zero-fill instead of removing it. + DefaultInitVector pos_flat; // ascending positions, concatenated in push order + DefaultInitVector pos_off; // query -> absolute offset into pos_flat + DefaultInitVector k_of; + DefaultInitVector phase_of; // emit_phase is ternary, so a byte is the whole range + + [[nodiscard]] auto size() const -> size_t { return n_; } + [[nodiscard]] auto positions() const -> size_t { return pos_n_; } + + auto clear() -> void { + n_ = 0; + pos_n_ = 0; + } + + auto reserve(size_t n_queries, size_t positions_per_query) -> void { + if (pos_off.size() < n_queries) { + pos_off.resize(n_queries); + k_of.resize(n_queries); + phase_of.resize(n_queries); + } + if (pos_flat.size() < n_queries * positions_per_query) { + pos_flat.resize(n_queries * positions_per_query); + } + } + + // Four preallocated writes, not four container appends. Callgrind on the pauli cell put + // vector::_M_range_insert at 110.6M instructions and vector::emplace_back at + // 52.2M -- 18% of this port's whole instruction delta -- for a push whose capacity is already + // reserved. `insert` cannot know that, so it re-derives the grow path per query; grow_() is the + // one place that checks, and it runs once per capacity doubling instead of once per push. + auto push(const PosT *pos, size_t k, int phase) -> void { + assert(phase >= -1 && phase <= 1 && "emit_phase is ternary: rotation_sign, or REAL_PARTS entry"); + const size_t n = n_; + const size_t at = pos_n_; + if (n == pos_off.size() || at + k > pos_flat.size()) { + grow_(k); + } + pos_off[n] = at; + // An explicit loop, not std::copy_n: k averages ~5 bytes here and copy_n compiles to a memcpy + // CALL, which callgrind counted 2.54M extra times for a copy smaller than its own prologue. + PosT *dst = pos_flat.data() + at; + for (size_t j = 0; j < k; ++j) { + dst[j] = pos[j]; + } + k_of[n] = static_cast(k); + phase_of[n] = static_cast(phase); + n_ = n + 1; + pos_n_ = at + k; + } + +private: + size_t n_ = 0; // queries pushed + size_t pos_n_ = 0; // positions pushed + + // Amortised doubling, and pos_flat grows by the larger of a double and what this push needs, so a + // single wide term cannot leave it short. + [[gnu::noinline]] auto grow_(size_t k) -> void { + if (n_ == pos_off.size()) { + const size_t want = (pos_off.size() * 2) + 64; + pos_off.resize(want); + k_of.resize(want); + phase_of.resize(want); + } + if (pos_n_ + k > pos_flat.size()) { + pos_flat.resize(std::max((pos_flat.size() * 2) + 256, pos_n_ + k)); + } + } +}; + +} // namespace monoprop::detail diff --git a/cpp/monoprop/detail/evolution/layer_build/QueryCodec.h b/cpp/monoprop/detail/evolution/layer_build/QueryCodec.h index a7439a30..a6da4977 100644 --- a/cpp/monoprop/detail/evolution/layer_build/QueryCodec.h +++ b/cpp/monoprop/detail/evolution/layer_build/QueryCodec.h @@ -66,6 +66,13 @@ struct QueryCodec { return CQ::push_mono(buf, mono, phase); } + // From ascending positions, which is what the partner merge hands the emit site; the dense overload + // above is for callers that hold only a bitset. + template + static auto push_positions(VecZ &buf, const PosU *pos, size_t k, int phase) -> size_t { + return CQ::push(buf, pos, k, phase); + } + // Identical in both formats: the value is one bit_cast word after the query's words. static auto push_value(VecZ &buf, double v) -> void { buf.push_back(encode_value(v)); } diff --git a/cpp/monoprop/detail/evolution/layer_build/Scan.h b/cpp/monoprop/detail/evolution/layer_build/Scan.h index 59c1a030..cce972b0 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Scan.h +++ b/cpp/monoprop/detail/evolution/layer_build/Scan.h @@ -30,6 +30,7 @@ #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/PartnerMerge.h" #include "monoprop/detail/evolution/layer_build/QueryCodec.h" #include "monoprop/detail/graph_encoding/MPGraphEncodingTypes.h" #include "monoprop/detail/mpi/MPIUtils.h" @@ -159,23 +160,65 @@ inline auto rotation_dynamic_gate(std::optional only_rotate_len_k, return true; } +// M⊕G as the emit site needs it. The dense form is unavoidable -- the owner hash folds every word and +// the basis sign reads the source bitset -- so the merge below runs ALONGSIDE it, not instead of it, +// and supplies k, d and the positions without a second sweep (paired_mode_count) or a third +// (push_mono's walk). +template +struct PartnerProduct { + Monomial new_mono; + size_t k = 0; // popcount(M⊕G) + size_t d = 0; // modes of M⊕G carrying BOTH Majoranas + size_t overlap = 0; // slots in both M and G, which cancel + int phase_factor = 0; +}; + // 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. -template +// hermitian_phase at emit; Pauli pauli_rotation_sign, already rotation-ready. `out_pos` receives the +// partner's ascending positions and needs capacity 2*NumModes; a spilled source row has no position +// array, so that case alone walks the dense partner back out. +// +// BOTH SHAPES WERE MEASURED, on the pauli cell over 2,455,950 emit calls (callgrind, jobs cg-sym4 and +// cg-sym5). Walking the dense partner instead -- find_first/find_next for the positions and the same +// running pairing test -- costs 229.1M MORE instructions than this merge, because the walk is a serial +// dependence chain through find_next where the merge streams two ascending arrays. The intuition that +// the merge is redundant work on top of a bitset that exists anyway is wrong: it is cheaper than +// reading that bitset back out. Do not replace it with the walk again. +template [[gnu::always_inline]] inline auto emit_term_products(const OperatorIndex &ham, size_t i, const typename A::GenContext &ctx, - Monomial &new_mono, - size_t &overlap, - int &phase_factor) -> void { - Monomial mono; - ham.for_each_position(i, [&](size_t pos) { mono.set(pos); }); + const GenT *gen_pos, + size_t gen_pop, + PosT *out_pos) -> PartnerProduct { const Monomial &gen = A::generator(ctx); - new_mono = mono ^ gen; - overlap = mono.count_and(gen); - phase_factor = A::rotation_sign(ctx, mono, new_mono); + PartnerProduct out; + Monomial mono; + if (const auto src = ham.row_positions(i); src.inlined()) { + out.k = merge_partner_positions(src.pos, src.count, gen_pos, gen_pop, out_pos, out.overlap, out.d); + for (size_t j = 0; j < src.count; ++j) { + mono.set(static_cast(src.pos[j])); + } + out.new_mono = mono ^ gen; + } + else { + ham.for_each_position(i, [&](size_t pos) { mono.set(pos); }); + out.new_mono = mono ^ gen; + out.overlap = mono.count_and(gen); + size_t prev = 0; + for (size_t b = out.new_mono.find_first(); b < out.new_mono.size(); b = out.new_mono.find_next(b)) { + if (out.k != 0 && (prev % 2 == 0) && b == prev + 1) { + ++out.d; + } + out_pos[out.k++] = static_cast(b); + prev = b; + } + } + out.phase_factor = A::rotation_sign(ctx, mono, out.new_mono); + return out; } +template struct FusedScanResult { std::vector cos_blocks; // ascending, disjoint, chunk order std::vector leader_queries; // size R: serialized leader queries per owner rank @@ -186,6 +229,11 @@ struct FusedScanResult { // leader_src / follower_src. Empty when capture_values is false. std::vector> leader_val; std::vector> follower_val; + // Self-owned queries, staged as positions instead of encoded into leader_queries[my_rank]: they are + // resolved inline and never reach a wire, so the codec is not on this leg. Order matches + // leader_src[my_rank] / follower_src[my_rank], which is the accumulation order. + SelfQueryStage leader_self; + SelfQueryStage follower_self; }; // Classify, cut off and emit in one pass over the anticommuting terms. Queries go to the owner of @@ -204,12 +252,12 @@ auto fused_find_and_collect(const MPOperator &op, size_t my_rank, bool capture_values = false, double *fused_scale_coeffs = nullptr, - double fused_scale_cos = 1.0) -> FusedScanResult { + double fused_scale_cos = 1.0) -> FusedScanResult { validate_only_rotate_len_k_(only_rotate_len_k, 2 * NumModes); const size_t gen_pop = gen.count(); const auto ectx = A::make_gen_context(gen); - FusedScanResult res; + FusedScanResult res; res.leader_queries.assign(rank_count, VecZ{}); res.leader_src.assign(rank_count, std::vector{}); res.follower_queries.assign(rank_count, VecZ{}); @@ -271,16 +319,29 @@ auto fused_find_and_collect(const MPOperator &op, auto &fv = res.follower_val; const OperatorIndex &ham = *op.store; + using RowPosT = typename OperatorIndex::PosT; + + // The generator's positions, once per gate: the merge's second input. + std::vector gen_pos; + gen_pos.reserve(gen_pop); + for (size_t b = gen.find_first(); b < gen.size(); b = gen.find_next(b)) { + gen_pos.push_back(static_cast(b)); + } + // 2*NumModes is the true bound: the partner's positions are distinct and below it. NOT + // thread_local: every access to one from a shared library goes through __tls_get_addr, which + // callgrind measured at 54.7M instructions on the pauli cell -- 6% of this port's delta -- to + // save one allocation per gate. + std::vector pbuf(2 * NumModes); - // Everything after a term survives the structural cutoff, from the dense partner emit built. + // Everything after a term survives the structural cutoff. Self-owned partners are staged as + // positions; only a remote owner's partner is encoded. auto push = [&](const Monomial &dense, - size_t mono_pop, - size_t overlap, - int phase_factor, + const RowPosT *pos, + size_t k, + int phase, size_t i, double v_src, bool is_follower) { - 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. // Must be the SAME function find_rank computes (MPIUtils.h) or a term is placed and queried // on different ranks, which duplicates a row silently; mpi_utils_tests.cpp asserts it. @@ -288,19 +349,15 @@ auto fused_find_and_collect(const MPOperator &op, if (rank_count != 1) { r_prime = monomial_hash(dense) % rank_count; } - if (is_follower) { - QueryCodec::push(fq[r_prime], dense, phase); - fs[r_prime].push_back(i); - if (capture_values) { - fv[r_prime].push_back(v_src); - } + if (r_prime == my_rank) { + (is_follower ? res.follower_self : res.leader_self).push(pos, k, phase); } else { - QueryCodec::push(lq[r_prime], dense, phase); - ls[r_prime].push_back(i); - if (capture_values) { - lv[r_prime].push_back(v_src); - } + QueryCodec::push_positions(is_follower ? fq[r_prime] : lq[r_prime], pos, k, phase); + } + (is_follower ? fs[r_prime] : ls[r_prime]).push_back(i); + if (capture_values) { + (is_follower ? fv[r_prime] : lv[r_prime]).push_back(v_src); } }; @@ -310,19 +367,17 @@ auto fused_find_and_collect(const MPOperator &op, if (!rotation_dynamic_gate(only_rotate_len_k, mono_pop, cut_st, abs_c)) { return; } - Monomial new_mono; - size_t overlap = 0; - int phase_factor = 0; - emit_term_products(ham, i, ectx, new_mono, overlap, phase_factor); + const auto p = emit_term_products(ham, i, ectx, gen_pos.data(), gen_pop, pbuf.data()); + assert(p.k == mono_pop + gen_pop - 2 * p.overlap && "the merge disagrees with the popcount identity"); // 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; // nullopt only for an opaque cutoff_fn_, which has no (k, d) form and must be invoked. - const auto keep = cutoff_eval.passes_from_dense(new_mono, new_pop); - const bool struct_pass = keep.value_or(false) || (!keep.has_value() && cutoff_eval(new_mono)); + const auto keep = cutoff_eval.passes_from_digest(p.k, p.d); + const bool struct_pass = keep.value_or(false) || (!keep.has_value() && cutoff_eval(p.new_mono)); if (!struct_pass && !cut_st.is_above_upper(abs_c)) { return; } - push(new_mono, mono_pop, overlap, phase_factor, i, v_src, is_follower); + const int phase = A::emit_phase(p.phase_factor, mono_pop, gen_pop, p.overlap); + push(p.new_mono, pbuf.data(), p.k, phase, i, v_src, is_follower); }; // Pass 1 and pass 2 stay fused over `nz`: splitting them regressed measurably, as `nz` spills L1 @@ -348,11 +403,11 @@ auto fused_find_and_collect(const MPOperator &op, n_foll); } if (rank_count == 1) { - // A hint only: 2 words covers k <= 6; wider terms grow the buffer rather than be reserved for. - const size_t qw = QueryCodec::kReserveWordsPerQuery; - lq[my_rank].reserve((n_anti - n_foll) * qw); + // A hint only, off the measured mean of 5.33 positions; wider terms grow the buffer. + const size_t pq = QueryCodec::kReservePositionsPerQuery; + res.leader_self.reserve(n_anti - n_foll, pq); ls[my_rank].reserve(n_anti - n_foll); - fq[my_rank].reserve(n_foll * qw); + res.follower_self.reserve(n_foll, pq); fs[my_rank].reserve(n_foll); } auto derive_coeff = [&](size_t i) -> std::pair { diff --git a/cpp/monoprop/detail/operator/OperatorIndex.h b/cpp/monoprop/detail/operator/OperatorIndex.h index 51323e22..be728232 100644 --- a/cpp/monoprop/detail/operator/OperatorIndex.h +++ b/cpp/monoprop/detail/operator/OperatorIndex.h @@ -527,6 +527,9 @@ class OperatorIndex { if (qk != static_cast(c)) { return false; } + // std::equal, i.e. a memcmp CALL, and MEASURED to be the right choice: replacing it with the + // obvious scalar loop cost 54.6M instructions on the pauli cell, because glibc's AVX2 memcmp + // beats a byte loop even at the ~5 PosT this compares. Do not "optimise" the call away again. return std::equal(q, q + qk, &rows_[(i * stride_) + 1]); } diff --git a/cpp/tests/mpi_utils_tests.cpp b/cpp/tests/mpi_utils_tests.cpp index 0f296507..b59f0426 100644 --- a/cpp/tests/mpi_utils_tests.cpp +++ b/cpp/tests/mpi_utils_tests.cpp @@ -130,6 +130,22 @@ auto check_bucket_ownership(const std::vector &buckets, size_t ranks, size } } +// The self-owned bucket is staged as positions, not encoded, so it is invisible to the walk above -- +// without this the r == my_rank arm of the routing decision goes unchecked. +auto check_self_ownership(const detail::SelfQueryStage<32> &stage, + size_t ranks, + size_t my_rank, + size_t &checked) -> void { + for (size_t q = 0; q < stage.size(); ++q) { + Monomial<32> mono; + for (size_t j = 0; j < stage.k_of[q]; ++j) { + mono.set(static_cast(stage.pos_flat[stage.pos_off[q] + j])); + } + BOOST_REQUIRE_EQUAL(find_rank<32>(mono, ranks), my_rank); + ++checked; + } +} + } // namespace // The scan hashes the partner it just built; find_rank hashes what the resolve side decoded off the @@ -155,6 +171,7 @@ BOOST_AUTO_TEST_CASE(mpi_utils_scan_routing_agrees_with_find_rank) { std::optional{0.3}); size_t checked = 0; + size_t self_checked = 0; for (const size_t ranks : {2U, 4U, 8U}) { const auto res = detail::fused_find_and_collect>(op, gen, @@ -168,9 +185,21 @@ BOOST_AUTO_TEST_CASE(mpi_utils_scan_routing_agrees_with_find_rank) { nullptr, 1.0); BOOST_REQUIRE_EQUAL(res.leader_queries.size(), ranks); + // The scan routes a self-owned partner to the stage, so bucket 0 must be empty here. + BOOST_REQUIRE(res.leader_queries[0].empty()); + BOOST_REQUIRE(res.follower_queries[0].empty()); check_bucket_ownership(res.leader_queries, ranks, checked); check_bucket_ownership(res.follower_queries, ranks, checked); + check_self_ownership(res.leader_self, ranks, /*my_rank=*/0, self_checked); + check_self_ownership(res.follower_self, ranks, /*my_rank=*/0, self_checked); } - // Without this the loop above passes trivially if the scan emitted nothing at all. - BOOST_TEST(checked > 1000U); + // Without this the loop above passes trivially if the scan emitted nothing. The floor is on the SUM + // because that is what is invariant across the split: the encoded counter alone fell to 797 of 1161 + // when the self-owned partners moved into the stage, with nothing going unchecked. Each arm still + // carries its own floor -- a routing bug sending everything one way leaves the sum intact -- and the + // message prints the measured 797/364 so those can be re-grounded rather than guessed. + BOOST_TEST_MESSAGE("encoded=" << checked << " staged=" << self_checked); + BOOST_TEST(checked + self_checked > 1000U); + BOOST_TEST(checked > 500U); + BOOST_TEST(self_checked > 200U); } diff --git a/cpp/tests/partner_merge_tests.cpp b/cpp/tests/partner_merge_tests.cpp new file mode 100644 index 00000000..c12a4cf4 --- /dev/null +++ b/cpp/tests/partner_merge_tests.cpp @@ -0,0 +1,168 @@ +// 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 partner merge against the dense XOR it replaces. Random pairs alone do not reach the cases that +// decide it: the merge's whole surface is how many of G's slots the source already holds, so every +// overlap in [0, gen_pop] is drawn on purpose, and the paired case is drawn separately because a +// fully-paired monomial is 94 in 20.9M on production models. + +#include + +#include +#include +#include +#include +#include + +#include "monoprop/algebra/AlgebraCommon.h" +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/evolution/layer_build/PartnerMerge.h" + +using namespace monoprop; + +namespace { + +constexpr size_t kN = 64; +constexpr size_t kBits = 2 * kN; + +// The dense reference: positions of M^G, ascending, straight off the bitset. +auto dense_partner(const Monomial &mono, const Monomial &gen) -> std::vector { + const Monomial nm = mono ^ gen; + std::vector out; + for (size_t b = nm.find_first(); b < nm.size(); b = nm.find_next(b)) { + out.push_back(b); + } + return out; +} + +auto positions_of(const Monomial &m) -> std::vector { + std::vector out; + for (size_t b = m.find_first(); b < m.size(); b = m.find_next(b)) { + out.push_back(static_cast(b)); + } + return out; +} + +// Checks the merge against the dense form on one pair, and returns the merged count so a caller can +// assert it saw work. Every field the emit site consumes is compared, not just the positions. +auto check_pair(const Monomial &mono, const Monomial &gen) -> size_t { + const auto src = positions_of(mono); + const auto gpos = positions_of(gen); + std::vector out(kBits); + size_t overlap = 0; + size_t d = 0; + const size_t k = + detail::merge_partner_positions(src.data(), src.size(), gpos.data(), gpos.size(), out.data(), overlap, d); + + const auto expect = dense_partner(mono, gen); + BOOST_REQUIRE_EQUAL(k, expect.size()); + for (size_t j = 0; j < k; ++j) { + BOOST_REQUIRE_EQUAL(static_cast(out[j]), expect[j]); + } + BOOST_REQUIRE_EQUAL(overlap, mono.count_and(gen)); + // d must be exactly what the cutoff digest would have folded out of the dense partner. + BOOST_REQUIRE_EQUAL(d, paired_mode_count(mono ^ gen)); + // The popcount identity the emit site asserts on. + BOOST_REQUIRE_EQUAL(k, mono.count() + gen.count() - (2 * overlap)); + return k; +} + +// A generator of `gen_pop` slots, and a source holding exactly `overlap` of them plus `extra` others. +auto build_case(std::mt19937_64 &rng, size_t gen_pop, size_t overlap, size_t extra) + -> std::pair, Monomial> { + std::vector all(kBits); + for (size_t i = 0; i < kBits; ++i) { + all[i] = i; + } + std::shuffle(all.begin(), all.end(), rng); + Monomial gen; + for (size_t i = 0; i < gen_pop; ++i) { + gen.set(all[i]); + } + Monomial mono; + for (size_t i = 0; i < overlap; ++i) { + mono.set(all[i]); // a slot G also holds: it cancels + } + for (size_t i = 0; i < extra; ++i) { + mono.set(all[gen_pop + i]); // disjoint from G: it survives + } + return {mono, gen}; +} + +} // namespace + +// Every overlap between the source and the generator, which is the branch the merge exists to take. +BOOST_AUTO_TEST_CASE(partner_merge_matches_dense_at_every_overlap) { + std::mt19937_64 rng(0xA11CEU); + size_t cases = 0; + size_t nonempty = 0; + for (size_t gen_pop = 1; gen_pop <= 6; ++gen_pop) { + for (size_t overlap = 0; overlap <= gen_pop; ++overlap) { + for (const size_t extra : {size_t{0}, size_t{1}, size_t{5}, size_t{20}}) { + for (size_t rep = 0; rep < 8; ++rep) { + const auto [mono, gen] = build_case(rng, gen_pop, overlap, extra); + BOOST_REQUIRE_EQUAL(mono.count_and(gen), overlap); // the case is the one intended + nonempty += (check_pair(mono, gen) != 0) ? 1 : 0; + ++cases; + } + } + } + } + BOOST_REQUIRE_EQUAL(cases, 6U * 4U * 8U + (1U + 2U + 3U + 4U + 5U + 6U) * 4U * 8U); + // Total cancellation (overlap == gen_pop, extra == 0) is the only empty partner, so most must not be. + BOOST_TEST(nonempty > 600U); +} + +// The paired population separately: d is the field the length cutoff's escape hatch rests on, and a +// uniform draw almost never produces a fully paired monomial. +BOOST_AUTO_TEST_CASE(partner_merge_d_matches_on_paired_monomials) { + std::mt19937_64 rng(0xBEEFU); + size_t paired_seen = 0; + for (size_t rep = 0; rep < 400; ++rep) { + // Whole modes only, so both slots of each are set and the result is fully paired. + std::vector modes(kN); + for (size_t i = 0; i < kN; ++i) { + modes[i] = i; + } + std::shuffle(modes.begin(), modes.end(), rng); + Monomial mono; + const size_t n_modes = 1 + (rng() % 8); + for (size_t i = 0; i < n_modes; ++i) { + mono.set(2 * modes[i]); + mono.set((2 * modes[i]) + 1); + } + Monomial gen; + const size_t g_modes = 1 + (rng() % 3); + for (size_t i = 0; i < g_modes; ++i) { + gen.set(2 * modes[kN - 1 - i]); + gen.set((2 * modes[kN - 1 - i]) + 1); + } + check_pair(mono, gen); + paired_seen += is_paired(mono ^ gen) ? 1 : 0; + } + // The draw is meant to land on the paired branch every time; a zero here means it stopped doing so. + BOOST_REQUIRE_EQUAL(paired_seen, 400U); +} + +// An empty generator and an empty source are both reachable (a truncated gate, a fresh row). +BOOST_AUTO_TEST_CASE(partner_merge_handles_empty_inputs) { + Monomial mono; + mono.set(4); + mono.set(5); + mono.set(70); + const Monomial empty; + BOOST_REQUIRE_EQUAL(check_pair(mono, empty), 3U); + BOOST_REQUIRE_EQUAL(check_pair(empty, mono), 3U); + BOOST_REQUIRE_EQUAL(check_pair(empty, empty), 0U); +} From 7e9f5188e7d58495839f753f6fa148a64b071e66 Mon Sep 17 00:00:00 2001 From: Aaron Miller Date: Mon, 24 Aug 2026 22:27:40 +0100 Subject: [PATCH 3/5] =?UTF-8?q?perf(evolution):=20=E2=9A=A1=20collapse=20t?= =?UTF-8?q?he=20query=20record=20to=20gap=20coding=20alone?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The record chose per-term between three closed forms -- raw lanes, gap coding and a raw bitmap -- behind a 2-bit mode field, and sized a stack Writer array from the worst case of the three. Replayed over the 106,368 captured query records, BITMAP was chosen 0 times and FIXED only on word-count ties, so the argmin was paying for two branches that never won. Gap coding alone, with the mode field's two bits and one of k's six spent on a narrower header: 11 bits against 14. gw = bit_width(max gap) <= kPosBits, so kPosBits + (k-1)*gw <= k*kPosBits and gap is never wider than raw lanes at any k. A raw mask does win once k*kPosBits > kBits, from k = 34 at 128 modes, where this record costs one word more -- above the k <= 32 that Pauli cutoff 16 implies and 11 above the widest k ever captured. That is the price of one code path, and sparse_record_documents_what_deleting_the_argmin_cost asserts it against the deleted encoder's own three formulas rather than leaving it to be rediscovered. Dropping the argmin also dropped the array it sized: gap coding emits bits monotonically, so the encoder streams into the output buffer through one accumulator and no longer zeroes 48 B per push. Bytes per record, predicted from the captured streams before the code existed and confirmed by it: hubbard c10 8.289 -> 8.038, pauli c12 15.519 -> 15.180, pauli c6 8.750 -> 8.438. Both formats decode the same positions on all 106,368 records (identical checksums). Callgrind over a standalone replay, where the binary is bit-exactly repeatable: push 345.1 -> 292.5 instructions per record (-15.2%), decode -3.0, every libc counter unchanged. The byte win is ~3% of a payload whose byte-proportional share of exchange is ~24%, so it is worth <= 0.2% of wall clock and no timing claim is made for it. Co-Authored-By: Claude Opus 5 --- .../detail/evolution/layer_build/QueryCodec.h | 2 +- .../evolution/layer_build/SparseQuery.h | 234 ++++++------------ cpp/tests/sparse_query_tests.cpp | 179 +++++++++++--- 3 files changed, 225 insertions(+), 190 deletions(-) diff --git a/cpp/monoprop/detail/evolution/layer_build/QueryCodec.h b/cpp/monoprop/detail/evolution/layer_build/QueryCodec.h index a6da4977..4a6e7411 100644 --- a/cpp/monoprop/detail/evolution/layer_build/QueryCodec.h +++ b/cpp/monoprop/detail/evolution/layer_build/QueryCodec.h @@ -43,7 +43,7 @@ struct QueryCodec { using PosT = typename CQ::PosT; // Words the QUERY at `off` occupies, NOT counting a trailing fused value word. Asked of the record - // rather than derived from k, which does not determine the width once the mode and gw can vary. + // rather than derived from k, which does not determine the width once gw can vary. [[nodiscard]] static auto query_words(const VecZ &buf, size_t off) -> size_t { return CQ::words_at(buf, off); } // Complete mode pairs among ascending positions, exposed so no caller names a concrete record type. diff --git a/cpp/monoprop/detail/evolution/layer_build/SparseQuery.h b/cpp/monoprop/detail/evolution/layer_build/SparseQuery.h index 330a182c..dab6e4e8 100644 --- a/cpp/monoprop/detail/evolution/layer_build/SparseQuery.h +++ b/cpp/monoprop/detail/evolution/layer_build/SparseQuery.h @@ -25,16 +25,20 @@ namespace monoprop::detail { // A variable-width query record: one word-aligned record per term holding the term's ascending set-bit -// positions. push() picks the argmin of three closed forms -- FIXED (k lanes of kPosBits), GAP (first -// position raw, then k-1 gaps of gw bits) and BITMAP (a raw kBits mask) -- so the record is never larger -// than any of the three, which is also what bounds Writer's kMaxWords. Record order is preserved -// everywhere, because it is the floating-point accumulation order (Resolve.h mints misses in it). +// positions, gap-coded. Record order is preserved everywhere, because it is the floating-point +// accumulation order (Resolve.h mints misses in it). // // Header in the low bits of word 0, then the payload LSB-first, both by explicit shift, never punning: -// [0..1] mode (FIXED/GAP/BITMAP) [2..3] phase+1 (emit_phase is TERNARY) [4..9] k, 63 escaping to a -// following 16-bit k then [4 bits] gw, in GAP mode only then payload. -// Header width is a design variable, not overhead: it decides which side of a 64-bit line a term falls -// on. +// [0..1] phase+1 (emit_phase is TERNARY) [2..6] k, 31 escaping to a following kLongKBits-wide k +// then [kGwBits] gw then pos[0] raw at kPosBits, then k-1 gaps of gw bits. +// +// ONE form, no mode field and no per-record argmin. Gap coding is never wider than raw lanes, because +// gw = bit_width(max gap) <= bit_width(kBits - 1) = kPosBits, hence kPosBits + (k-1)*gw <= k*kPosBits. A +// raw kBits mask is narrower only once k*kPosBits > kBits, i.e. from k = 34 at 128 modes, where this +// record costs one word more: measured on 0 of 106,368 captured records (max k = 23, at pauli cutoff +// 12), and that is the price of one code path. Dropping the argmin also dropped the stack array it +// sized -- gap coding emits bits monotonically, so the encoder streams straight into `buf` and no longer +// zeroes 48 B per push. template struct SparseQuery { using PosT = uint16_t; @@ -45,32 +49,30 @@ struct SparseQuery { //: Bits for one raw position in [0, 2*NumModes); compile-time, so the lane width is free. static constexpr size_t kPosBits = static_cast(std::bit_width(kBits - 1)); - static constexpr uint64_t kModeFixed = 0; - static constexpr uint64_t kModeGap = 1; - static constexpr uint64_t kModeBitmap = 2; - - static constexpr size_t kModeBits = 2; static constexpr size_t kPhaseBits = 2; - static constexpr size_t kKBits = 6; - static constexpr size_t kLongKBits = 16; + static constexpr size_t kKBits = 5; + //: k is a popcount of a kBits bitset, so the escape can never need more than this. + static constexpr size_t kLongKBits = static_cast(std::bit_width(kBits)); static constexpr size_t kGwBits = 4; static constexpr size_t kKEscape = (1U << kKBits) - 1U; - static constexpr size_t kBaseHeaderBits = kModeBits + kPhaseBits + kKBits; - static_assert(kPosBits <= (1U << kGwBits) - 1U, "gw <= kPosBits must fit the header's 4-bit gap-width field"); + static constexpr size_t kHeaderBits = kPhaseBits + kKBits + kGwBits; + static_assert(kPosBits <= (1U << kGwBits) - 1U, "gw <= kPosBits must fit the header's gap-width field"); + static_assert(kHeaderBits + kLongKBits <= 64, "the widest header must be readable from word 0 alone"); - static constexpr size_t kMaxPositions = 65535; - //: Words in a full occupancy mask. ceil, NOT kBits/64: a C++ caller may pick any width (LiH: 24). - static constexpr size_t kMaskWords = (kBits + 63U) / 64U; - //: BITMAP bounds every mode from above, because the encoder takes the minimum of the three. - static constexpr size_t kMaxWords = (kBits + kBaseHeaderBits + kLongKBits + kGwBits + 63U) / 64U + 1U; + //: A popcount cannot exceed the width, which the old 65535 never said. + static constexpr size_t kMaxPositions = kBits; // ---- bit stream ------------------------------------------------------------------------------- + // Streams into `buf`: one accumulator, flushed when a word fills, in place of an array sized by the + // worst case of three encodings. struct Writer { - uint64_t w[kMaxWords] = {}; - size_t nbits = 0; + VecZ &buf; + uint64_t cur = 0; + size_t nbits = 0; // bits held in cur, always < 64 + size_t words = 0; - constexpr auto put(uint64_t v, size_t width) noexcept -> void { + [[gnu::always_inline]] auto put(uint64_t v, size_t width) noexcept -> void { if (width == 0) { assert(v == 0 && "a zero-width field cannot carry a value"); return; @@ -80,15 +82,25 @@ struct SparseQuery { if (width < 64) { v &= (uint64_t{1} << width) - 1U; } - const size_t word = nbits >> 6U; - const size_t off = nbits & 63U; - assert(word < kMaxWords && "record overran its worst-case word bound"); - w[word] |= v << off; - // off > 0 is implied here (width <= 64), so the shift is in [1, 63]; `v >> 64` would be UB. - if (off + width > 64) { - w[word + 1] |= v >> (64U - off); + cur |= v << nbits; + if (nbits + width < 64) { + nbits += width; + return; + } + buf.push_back(static_cast(cur)); + ++words; + // nbits == 0 only at width == 64, where every bit is already in cur; `v >> 64` would be UB. + cur = (nbits == 0) ? 0 : (v >> (64U - nbits)); + nbits = nbits + width - 64U; + } + + auto flush() noexcept -> void { + if (nbits != 0) { + buf.push_back(static_cast(cur)); + ++words; + cur = 0; + nbits = 0; } - nbits += width; } }; @@ -115,7 +127,6 @@ struct SparseQuery { // ---- header ----------------------------------------------------------------------------------- struct Header { - uint64_t mode = 0; int phase = 0; size_t k = 0; size_t gw = 0; @@ -126,46 +137,30 @@ struct SparseQuery { [[nodiscard]] static auto header_at(const VecZ &buf, size_t off) noexcept -> Header { const auto w0 = static_cast(buf[off]); Header h; - h.mode = w0 & 0x3U; - h.phase = static_cast((w0 >> kModeBits) & 0x3U) - 1; - h.k = static_cast((w0 >> (kModeBits + kPhaseBits)) & kKEscape); - h.bits = kBaseHeaderBits; + h.phase = static_cast(w0 & 0x3U) - 1; + h.k = static_cast((w0 >> kPhaseBits) & kKEscape); + h.bits = kPhaseBits + kKBits; if (h.k == kKEscape) { - h.k = static_cast((w0 >> h.bits) & 0xFFFFU); + h.k = static_cast((w0 >> h.bits) & ((uint64_t{1} << kLongKBits) - 1U)); h.bits += kLongKBits; } - if (h.mode == kModeGap) { - h.gw = static_cast((w0 >> h.bits) & 0xFU); - h.bits += kGwBits; - } + h.gw = static_cast((w0 >> h.bits) & ((uint64_t{1} << kGwBits) - 1U)); + h.bits += kGwBits; return h; } - [[nodiscard]] static constexpr auto header_bits_for(size_t k, uint64_t mode) noexcept -> size_t { - return kBaseHeaderBits + ((k >= kKEscape) ? kLongKBits : 0U) + ((mode == kModeGap) ? kGwBits : 0U); + [[nodiscard]] static constexpr auto header_bits_for(size_t k) noexcept -> size_t { + return kHeaderBits + ((k >= kKEscape) ? kLongKBits : 0U); } - [[nodiscard]] static constexpr auto fixed_bits(size_t k) noexcept -> size_t { - return header_bits_for(k, kModeFixed) + k * kPosBits; - } [[nodiscard]] static constexpr auto gap_bits(size_t k, size_t gw) noexcept -> size_t { - return header_bits_for(k, kModeGap) + ((k == 0) ? 0U : kPosBits + (k - 1U) * gw); - } - [[nodiscard]] static constexpr auto bitmap_bits(size_t k) noexcept -> size_t { - return header_bits_for(k, kModeBitmap) + kBits; + return header_bits_for(k) + ((k == 0) ? 0U : kPosBits + (k - 1U) * gw); } [[nodiscard]] static constexpr auto words_of(size_t bits) noexcept -> size_t { return (bits + 63U) / 64U; } - //: The record's word count, from the header alone: k does not determine it, mode and gw do too. + //: The record's word count, from the header alone: k does not determine it, gw does too. [[nodiscard]] static constexpr auto words_of_header(const Header &h) noexcept -> size_t { - switch (h.mode) { - case kModeGap: - return words_of(gap_bits(h.k, h.gw)); - case kModeBitmap: - return words_of(bitmap_bits(h.k)); - default: - return words_of(fixed_bits(h.k)); - } + return words_of(gap_bits(h.k, h.gw)); } [[nodiscard]] static auto words_at(const VecZ &buf, size_t off) noexcept -> size_t { @@ -197,30 +192,14 @@ struct SparseQuery { // type is narrower than the wire's below 129 modes, and the encoding does not depend on it. template static auto push(VecZ &buf, const PosU *pos, size_t k, int phase) -> size_t { - assert(k <= kMaxPositions && "term has more positions than the record's k field can hold"); + assert(k <= kMaxPositions && "term has more positions than the record's width admits"); assert(phase >= -1 && phase <= 1 && "emit_phase is ternary: rotation_sign, or REAL_PARTS entry"); for (size_t j = 1; j < k; ++j) { assert(pos[j] > pos[j - 1] && "positions must be strictly ascending"); } const size_t gw = gap_width(pos, k); - const size_t wf = words_of(fixed_bits(k)); - const size_t wg = words_of(gap_bits(k, gw)); - const size_t wb = words_of(bitmap_bits(k)); - - uint64_t mode = kModeFixed; - size_t want = wf; - if (wg < want) { - mode = kModeGap; - want = wg; - } - if (wb < want) { - mode = kModeBitmap; - want = wb; - } - - Writer w; - w.put(mode, kModeBits); + Writer w{buf}; w.put(static_cast(phase + 1), kPhaseBits); if (k >= kKEscape) { w.put(kKEscape, kKBits); @@ -229,41 +208,16 @@ struct SparseQuery { else { w.put(static_cast(k), kKBits); } - - if (mode == kModeGap) { - w.put(static_cast(gw), kGwBits); - if (k != 0) { - w.put(static_cast(pos[0]), kPosBits); - for (size_t j = 1; j < k; ++j) { - w.put(static_cast(pos[j] - pos[j - 1] - 1U), gw); - } - } - } - else if (mode == kModeBitmap) { - // The trailing partial word is NOT optional: at kBits 24 a `kBits / 64` loop writes nothing - // while bitmap_bits() still charges 24 bits, and the decoder then recovers no positions. - uint64_t mask[kMaskWords] = {}; - for (size_t j = 0; j < k; ++j) { - const auto p = static_cast(pos[j]); - mask[p >> 6U] |= uint64_t{1} << (p & 63U); - } - for (size_t done = 0, i = 0; done < kBits; ++i) { - const size_t chunk = (kBits - done < 64U) ? (kBits - done) : 64U; - w.put(mask[i], chunk); - done += chunk; - } - } - else { - for (size_t j = 0; j < k; ++j) { - w.put(static_cast(pos[j]), kPosBits); + w.put(static_cast(gw), kGwBits); + if (k != 0) { + w.put(static_cast(pos[0]), kPosBits); + for (size_t j = 1; j < k; ++j) { + w.put(static_cast(pos[j] - pos[j - 1] - 1U), gw); } } - - assert(words_of(w.nbits) == want && "encoder wrote a different width than it costed"); - for (size_t i = 0; i < want; ++i) { - buf.push_back(static_cast(w.w[i])); - } - return want; + w.flush(); + assert(w.words == words_of(gap_bits(k, gw)) && "encoder wrote a different width than it costed"); + return w.words; } // ---- decode ----------------------------------------------------------------------------------- @@ -273,34 +227,12 @@ struct SparseQuery { static auto read_positions(const VecZ &buf, size_t off, OutT *out) -> size_t { const Header h = header_at(buf, off); Reader r{buf, off, h.bits}; - if (h.mode == kModeGap) { - if (h.k != 0) { - auto prev = static_cast(r.get(kPosBits)); - out[0] = static_cast(prev); - for (size_t j = 1; j < h.k; ++j) { - prev += static_cast(r.get(h.gw)) + 1U; - out[j] = static_cast(prev); - } - } - } - else if (h.mode == kModeBitmap) { - // Symmetric with the encoder, trailing partial word included -- see the note there. - size_t n = 0; - for (size_t done = 0, i = 0; done < kBits; ++i) { - const size_t chunk = (kBits - done < 64U) ? (kBits - done) : 64U; - uint64_t word = r.get(chunk); - while (word != 0) { - const auto b = static_cast(std::countr_zero(word)); - out[n++] = static_cast(done + b); - word &= word - 1U; - } - done += chunk; - } - assert(n == h.k && "bitmap popcount disagrees with the record's k"); - } - else { - for (size_t j = 0; j < h.k; ++j) { - out[j] = static_cast(r.get(kPosBits)); + if (h.k != 0) { + auto prev = static_cast(r.get(kPosBits)); + out[0] = static_cast(prev); + for (size_t j = 1; j < h.k; ++j) { + prev += static_cast(r.get(h.gw)) + 1U; + out[j] = static_cast(prev); } } const size_t next = off + words_of_header(h); @@ -312,9 +244,6 @@ struct SparseQuery { template [[nodiscard]] static auto check_header(const VecZ &buf, size_t off, const OutT *pos) -> bool { const Header h = header_at(buf, off); - if (h.mode > kModeBitmap) { - return false; // an unknown mode is a format the reader does not understand - } if (h.phase < -1 || h.phase > 1) { return false; } @@ -326,18 +255,13 @@ struct SparseQuery { if (h.k != 0 && static_cast(pos[h.k - 1]) >= kBits) { return false; } - if (h.mode == kModeGap && h.k > 1) { - // gw is the MAXIMUM gap width: too small truncates a gap silently, too large wastes bits. - size_t g = 0; - for (size_t j = 1; j < h.k; ++j) { - const auto b = static_cast(std::bit_width(static_cast(pos[j] - pos[j - 1] - 1))); - g = (b > g) ? b : g; - } - if (g != h.gw) { - return false; - } + // gw is the MAXIMUM gap width: too small truncates a gap silently, too large wastes bits. + size_t g = 0; + for (size_t j = 1; j < h.k; ++j) { + const auto b = static_cast(std::bit_width(static_cast(pos[j] - pos[j - 1] - 1))); + g = (b > g) ? b : g; } - return true; + return g == h.gw; } // d, recomputed rather than carried: ascending order makes a pair an even position then its successor. @@ -376,7 +300,7 @@ struct SparseQuery { return next; } - // Encode from a dense monomial, which is what the scan holds: the partner is built densely anyway. + // Encode from a dense monomial, for callers that hold only a bitset; the emit path merges positions. static auto push_mono(VecZ &buf, const Monomial &mono, int phase) -> size_t { const size_t k = mono.count(); if (k <= kStackPositions) { diff --git a/cpp/tests/sparse_query_tests.cpp b/cpp/tests/sparse_query_tests.cpp index a0a8f544..00e6b7ab 100644 --- a/cpp/tests/sparse_query_tests.cpp +++ b/cpp/tests/sparse_query_tests.cpp @@ -92,8 +92,8 @@ auto strided(size_t k, size_t start, size_t step, size_t universe) -> std::vecto return v; } -// Uniform draws are what can reach bitmap mode -- gap coding wins every regular pattern -- but see -// sparse_record_actually_exercises_its_bitmap_mode: at narrow widths they cannot reach it either. +// Uniform draws are the widest gap widths, which is the case a regular pattern never reaches; see +// sparse_record_reaches_the_widest_gap_width for why a narrow universe needs a run plus one outlier. auto scattered(size_t k, size_t universe, std::mt19937_64 &rng) -> std::vector { std::vector pool(universe); for (size_t j = 0; j < universe; ++j) { @@ -168,10 +168,11 @@ BOOST_AUTO_TEST_CASE(sparse_record_handles_widths_that_are_not_whole_words) { } } -BOOST_AUTO_TEST_CASE(sparse_record_actually_exercises_its_bitmap_mode) { - // Two generators because uniform draws CANNOT reach bitmap at kBits=24: bitmap is one word there, so - // gap must need two (19 + (k-1)*gw > 64), which only a dense run plus one far outlier forces. - const auto count_bitmap = [](auto tag, size_t universe) { +BOOST_AUTO_TEST_CASE(sparse_record_reaches_the_widest_gap_width) { + // gw == kPosBits is the record's worst case and uniform draws CANNOT reach it at kBits=24: it needs + // one gap of at least half the universe, which only a dense run plus one far outlier forces. The + // shape of the input is a selection rule, so this generator is kept even though bitmap mode is gone. + const auto count_widest = [](auto tag, size_t universe) { using SQ = SparseQuery; size_t used = 0; size_t bad = 0; @@ -180,14 +181,16 @@ BOOST_AUTO_TEST_CASE(sparse_record_actually_exercises_its_bitmap_mode) { return; } VecZ buf; - (void)SQ::push(buf, pos.data(), pos.size(), 1); - if (SQ::header_at(buf, 0).mode != SQ::kModeBitmap) { + const size_t w = SQ::push(buf, pos.data(), pos.size(), 1); + const size_t gw = SQ::gap_width(pos.data(), pos.size()); + if (gw != SQ::kPosBits) { return; } ++used; std::vector back(pos.size()); SQ::read_positions(buf, 0, back.data()); - bad += static_cast(back != pos || SQ::k_at(buf, 0) != pos.size()); + bad += static_cast(back != pos || SQ::k_at(buf, 0) != pos.size() + || w != SQ::words_of(SQ::gap_bits(pos.size(), gw))); }; std::mt19937_64 rng(0xB1747U ^ universe); for (size_t trial = 0; trial < 600; ++trial) { @@ -206,9 +209,9 @@ BOOST_AUTO_TEST_CASE(sparse_record_actually_exercises_its_bitmap_mode) { } return std::pair{used, bad}; }; - const auto narrow = count_bitmap(std::integral_constant{}, 24); - const auto partial = count_bitmap(std::integral_constant{}, 500); - const auto bucket = count_bitmap(std::integral_constant{}, 256); + const auto narrow = count_widest(std::integral_constant{}, 24); + const auto partial = count_widest(std::integral_constant{}, 500); + const auto bucket = count_widest(std::integral_constant{}, 256); BOOST_TEST(narrow.first > 0U); BOOST_TEST(partial.first > 0U); BOOST_TEST(bucket.first > 0U); @@ -217,8 +220,9 @@ BOOST_AUTO_TEST_CASE(sparse_record_actually_exercises_its_bitmap_mode) { BOOST_TEST(bucket.second == 0U); } -BOOST_AUTO_TEST_CASE(sparse_record_survives_the_six_bit_k_escape) { - for (const size_t k : {size_t{62}, size_t{63}, size_t{64}, size_t{200}}) { +BOOST_AUTO_TEST_CASE(sparse_record_survives_the_five_bit_k_escape) { + // The escape is at k = 31 now, not 63: both boundaries are here so a field-width change is caught. + for (const size_t k : {size_t{30}, size_t{31}, size_t{32}, size_t{62}, size_t{63}, size_t{64}, size_t{200}}) { const auto pos = strided(k, 0, 3, 2048); BOOST_REQUIRE_EQUAL(pos.size(), k); differential<1024>(pos, 1); @@ -226,7 +230,8 @@ BOOST_AUTO_TEST_CASE(sparse_record_survives_the_six_bit_k_escape) { } BOOST_AUTO_TEST_CASE(sparse_record_bounds_the_fully_paired_term) { - // Every bit set: 514 words in FIXED at NumModes=1024, one word in BITMAP. Hence the argmin. + // Every bit set means every gap is 0, so gw is 0 and the payload is one raw position: 23 header + // bits + 11 = one word, where raw lanes would take 514. This is what pays for deleting the argmin. std::vector all(2048); for (size_t j = 0; j < all.size(); ++j) { all[j] = static_cast(j); @@ -235,25 +240,70 @@ BOOST_AUTO_TEST_CASE(sparse_record_bounds_the_fully_paired_term) { BOOST_TEST(sw == 1U); } -BOOST_AUTO_TEST_CASE(sparse_record_never_exceeds_the_dense_words_it_replaces) { - const auto check = [](auto tag, size_t universe, size_t dense_words) { +BOOST_AUTO_TEST_CASE(sparse_record_never_exceeds_its_own_raw_lanes) { + // The one width guarantee that survives deleting the argmin, and it is exhaustive rather than + // sampled: gw = bit_width(max gap) <= kPosBits, so kPosBits + (k-1)*gw <= k*kPosBits at every k. + const auto check = [](auto tag) { using SQ = SparseQuery; - std::mt19937_64 rng(0xC0FFEE ^ universe); - for (size_t k = 0; k <= universe; k += std::max(1, universe / 37)) { - const auto pos = scattered(k, universe, rng); - VecZ buf; - const size_t w = SQ::push(buf, pos.data(), pos.size(), 1); - BOOST_TEST(w <= dense_words + 1, - "k=" << k << " at U=" << universe << " took " << w << " words vs dense " << dense_words); + size_t cells = 0; + size_t bad = 0; + for (size_t k = 0; k <= SQ::kMaxPositions; ++k) { + for (size_t gw = 0; gw <= SQ::kPosBits; ++gw) { + const size_t lanes = SQ::words_of(SQ::header_bits_for(k) + (k * SQ::kPosBits)); + bad += static_cast(SQ::words_of(SQ::gap_bits(k, gw)) > lanes); + ++cells; + } } + return std::pair{cells, bad}; }; - check(std::integral_constant{}, 64, 2); - check(std::integral_constant{}, 256, 5); - check(std::integral_constant{}, 500, 9); + const auto narrow = check(std::integral_constant{}); + const auto bucket = check(std::integral_constant{}); + const auto wide = check(std::integral_constant{}); + BOOST_TEST(narrow.second == 0U); + BOOST_TEST(bucket.second == 0U); + BOOST_TEST(wide.second == 0U); + // A guarded loop that asserted nothing would pass the three above; these are the cell counts. + BOOST_TEST(narrow.first == 25U * 6U); + BOOST_TEST(bucket.first == 257U * 9U); + BOOST_TEST(wide.first == 501U * 10U); +} + +BOOST_AUTO_TEST_CASE(sparse_record_documents_what_deleting_the_argmin_cost) { + // Deleting FIXED and BITMAP has a price, and this pins it to a number rather than leaving it to be + // rediscovered. The three formulas below are the DELETED encoder's, with its own 10-bit header and + // 16-bit k escape, so the comparison is against what actually shipped in #263. + using SQ = SparseQuery<128>; + const auto old_words = [](size_t k, size_t gw) { + const size_t h = 10U + ((k >= 63U) ? 16U : 0U); + return std::min({SQ::words_of(h + (k * SQ::kPosBits)), + SQ::words_of(h + 4U + (k == 0 ? 0U : SQ::kPosBits + ((k - 1U) * gw))), + SQ::words_of(h + SQ::kBits)}); + }; + // Nothing in the supported envelope loses: Pauli cutoff 16 bounds k at 32, and the widest k ever + // captured is 23 (pauli c12, 45,296 records). + size_t crossings = 0; + for (size_t k = 0; k <= 33U; ++k) { + for (size_t gw = 0; gw <= SQ::kPosBits; ++gw) { + crossings += static_cast(SQ::words_of(SQ::gap_bits(k, gw)) > old_words(k, gw)); + } + } + BOOST_TEST(crossings == 0U); + + // Above it, a raw mask wins, and 34 is where. If a field width changes, this number moves and says so. + size_t first = SQ::kMaxPositions + 1U; + for (size_t k = 0; k <= SQ::kMaxPositions && first > SQ::kMaxPositions; ++k) { + for (size_t gw = 0; gw <= SQ::kPosBits; ++gw) { + if (SQ::words_of(SQ::gap_bits(k, gw)) > old_words(k, gw)) { + first = k; + break; + } + } + } + BOOST_TEST(first == 34U); } BOOST_AUTO_TEST_CASE(sparse_record_walks_a_multi_query_buffer_exactly) { - // Mixed width and mixed mode, which is the case a hardcoded stride gets wrong. + // Mixed width, which is the case a hardcoded stride gets wrong. using SQ = SparseQuery<128>; using QC = QueryCodec<128>; const QueryLayout layout{/*fused=*/false}; @@ -287,8 +337,9 @@ BOOST_AUTO_TEST_CASE(sparse_record_walks_a_multi_query_buffer_exactly) { BOOST_TEST(off == buf.size()); } -BOOST_AUTO_TEST_CASE(sparse_record_picks_the_smallest_of_its_three_modes) { - // Load-bearing: gap ALONE is 8.56 B/term against fixed lanes' 8.00 on uniform draws at 250 modes. +BOOST_AUTO_TEST_CASE(sparse_record_is_exactly_the_gap_code_it_costed) { + // What the argmin test became: there is one closed form now, so the encoder's width and the costing + // function must agree on every draw, including the uniform ones that used to select other modes. using SQ = SparseQuery<128>; std::mt19937_64 rng(12345); for (size_t trial = 0; trial < 400; ++trial) { @@ -297,10 +348,10 @@ BOOST_AUTO_TEST_CASE(sparse_record_picks_the_smallest_of_its_three_modes) { VecZ buf; const size_t w = SQ::push(buf, pos.data(), pos.size(), 1); const size_t gwid = SQ::gap_width(pos.data(), pos.size()); - const size_t best = std::min({SQ::words_of(SQ::fixed_bits(pos.size())), - SQ::words_of(SQ::gap_bits(pos.size(), gwid)), - SQ::words_of(SQ::bitmap_bits(pos.size()))}); - BOOST_TEST(w == best, "k=" << k << " wrote " << w << " words, best was " << best); + BOOST_TEST(w == SQ::words_of(SQ::gap_bits(pos.size(), gwid)), + "k=" << k << " wrote " << w << " words, costed " + << SQ::words_of(SQ::gap_bits(pos.size(), gwid))); + BOOST_TEST(w <= SQ::words_of(SQ::header_bits_for(pos.size()) + (pos.size() * SQ::kPosBits))); } } @@ -480,3 +531,63 @@ BOOST_AUTO_TEST_CASE(sparse_record_fused_value_channel_is_bit_exact_and_reusable QC::build_fused(empty, {}, dirty); BOOST_TEST(dirty.empty()); } + +// Positions with exactly k entries whose widest gap is exactly `gw`: one gap of 2^(gw-1) -- the smallest +// value of that bit width -- then a contiguous run. Empty if the shape does not fit the universe. +namespace { +auto gap_shaped(size_t k, size_t gw, size_t universe) -> std::vector { + std::vector v; + if (k == 0) { + return v; + } + if (gw == 0 || k == 1) { + if (gw != 0 || k > universe) { + return v; + } + for (size_t j = 0; j < k; ++j) { + v.push_back(static_cast(j)); + } + return v; + } + const size_t second = (size_t{1} << (gw - 1)) + 1U; // gap value 2^(gw-1), i.e. bit_width == gw + if (second + (k - 2U) >= universe) { + return v; + } + v.push_back(0); + for (size_t j = 0; j + 1 < k; ++j) { + v.push_back(static_cast(second + j)); + } + return v; +} +} // namespace + +BOOST_AUTO_TEST_CASE(sparse_record_round_trips_every_reachable_k_and_gap_width) { + // The whole (k, gw) surface of the one remaining form, constructed rather than drawn: a random draw + // reaches neither gw = kPosBits nor the escape boundary. differential() carries ten assertions per + // cell, and `cells` is here so a shape that stops fitting cannot silently empty the loop. + using SQ = SparseQuery<128>; + size_t cells = 0; + for (size_t k = 0; k <= 40U; ++k) { + for (size_t gw = 0; gw <= SQ::kPosBits; ++gw) { + if (k < 2U && gw > 0U) { + continue; // one gap width is reachable below k = 2, so the other rows are the same cell + } + const auto pos = gap_shaped(k, gw, SQ::kBits); + if (pos.size() != k) { + continue; + } + BOOST_REQUIRE_EQUAL(SQ::gap_width(pos.data(), k), k < 2 ? 0U : gw); + const size_t w = differential<128>(pos, (k % 3U) == 0U ? 0 : ((k % 3U) == 1U ? 1 : -1)); + BOOST_TEST(w == SQ::words_of(SQ::gap_bits(k, k < 2 ? 0U : gw))); + ++cells; + } + } + BOOST_TEST(cells == 353U); + + // The width boundary: k = kBits is a fully paired term, one word because every gap is 0. + for (const size_t k : {size_t{254}, size_t{255}, size_t{256}}) { + const auto pos = gap_shaped(k, 0, SQ::kBits); + BOOST_REQUIRE_EQUAL(pos.size(), k); + BOOST_TEST(differential<128>(pos, 1) == 1U); + } +} From 96e2ec0caaa300f1652b791c1eb51da259999c26 Mon Sep 17 00:00:00 2001 From: Aaron Miller Date: Mon, 24 Aug 2026 23:34:54 +0100 Subject: [PATCH 4/5] test(evolution): :bug: stage the self-resolve fixture as positions #267 landed on main while #263 was open, and its self_resolve_mark_bounded_by_combined_size fixture feeds the engine through detail::query_push -- the dense record #263 retires. The rebase onto that main left the call in place, so origin/pr/query-wire-v3 does not compile: query_push now lives in cpp/tests/dense_query_reference.h as test_ref::query_push, an oracle for the codec differential and not an engine entry point. Ported rather than renamed. Under the positions-staged self leg the fixture would have tripped resolve_self_queries' own assertion -- "a self-owned query was encoded instead of staged" -- because at my_rank == 0 the engine reads self_stage_ and requires queries_r[my_rank] empty. It now pushes each term's ascending positions onto self_stage_, which is what the scan does, and every assertion the case made about marking past combined_size is unchanged. Co-Authored-By: Claude Opus 5 --- cpp/tests/evolution_detail_tests.cpp | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/cpp/tests/evolution_detail_tests.cpp b/cpp/tests/evolution_detail_tests.cpp index 91f2662b..bade37f4 100644 --- a/cpp/tests/evolution_detail_tests.cpp +++ b/cpp/tests/evolution_detail_tests.cpp @@ -154,8 +154,17 @@ BOOST_AUTO_TEST_CASE(self_resolve_mark_bounded_by_combined_size) { matched, combined_size, RecordingSink{}); - detail::query_push<8>(eng.queries_r[0], terms[1], 1); - detail::query_push<8>(eng.queries_r[0], terms[5], -1); + // The self leg is staged as positions, never encoded, so this feeds the stage the scan would fill. + using Eng = detail::LayerBuildEngine<8, RecordingSink>; + const auto stage_self = [&eng](const Monomial<8> &m, int phase) { + std::vector pos; + for (size_t b = m.find_first(); b < m.size(); b = m.find_next(b)) { + pos.push_back(static_cast(b)); + } + eng.self_stage_.push(pos.data(), pos.size(), phase); + }; + stage_self(terms[1], 1); + stage_self(terms[5], -1); eng.src_idx_r[0] = {0, 2}; eng.resolve_self_queries(/*is_leader_pass=*/true); From 9ff9d73238f6e5a1a1eb96ef9fecad8861671e44 Mon Sep 17 00:00:00 2001 From: Aaron Miller Date: Tue, 25 Aug 2026 09:01:54 +0100 Subject: [PATCH 5/5] style: :art: apply clang-format to the lines this stack added prek pins clang-format v21.1.0 and three of the lines PR A and PR B wrote were not in its shape: two comment columns in PartnerMerge.h, a four-line parameter list in mpi_utils_tests.cpp that fits on one, and a BOOST_TEST message split a token earlier than needed. `git diff -w` reports PartnerMerge.h as empty and the other two as line joins with every token conserved, so nothing here can reach codegen -- it is the lint leg only, which was the one red left on the PR. Co-Authored-By: Claude Opus 5 --- cpp/monoprop/detail/evolution/layer_build/PartnerMerge.h | 4 ++-- cpp/tests/mpi_utils_tests.cpp | 6 ++---- cpp/tests/sparse_query_tests.cpp | 3 +-- 3 files changed, 5 insertions(+), 8 deletions(-) diff --git a/cpp/monoprop/detail/evolution/layer_build/PartnerMerge.h b/cpp/monoprop/detail/evolution/layer_build/PartnerMerge.h index 712e4cc0..8b5bc668 100644 --- a/cpp/monoprop/detail/evolution/layer_build/PartnerMerge.h +++ b/cpp/monoprop/detail/evolution/layer_build/PartnerMerge.h @@ -101,8 +101,8 @@ struct SelfQueryStage { // DefaultInitVector, the allocator Resolve.h already uses for exactly this: a plain vector's resize // VALUE-initialises, so sizing pos_flat ahead would memset every byte a push is about to overwrite // -- trading the append cost for a per-gate zero-fill instead of removing it. - DefaultInitVector pos_flat; // ascending positions, concatenated in push order - DefaultInitVector pos_off; // query -> absolute offset into pos_flat + DefaultInitVector pos_flat; // ascending positions, concatenated in push order + DefaultInitVector pos_off; // query -> absolute offset into pos_flat DefaultInitVector k_of; DefaultInitVector phase_of; // emit_phase is ternary, so a byte is the whole range diff --git a/cpp/tests/mpi_utils_tests.cpp b/cpp/tests/mpi_utils_tests.cpp index b59f0426..45fdaf9e 100644 --- a/cpp/tests/mpi_utils_tests.cpp +++ b/cpp/tests/mpi_utils_tests.cpp @@ -132,10 +132,8 @@ auto check_bucket_ownership(const std::vector &buckets, size_t ranks, size // The self-owned bucket is staged as positions, not encoded, so it is invisible to the walk above -- // without this the r == my_rank arm of the routing decision goes unchecked. -auto check_self_ownership(const detail::SelfQueryStage<32> &stage, - size_t ranks, - size_t my_rank, - size_t &checked) -> void { +auto check_self_ownership(const detail::SelfQueryStage<32> &stage, size_t ranks, size_t my_rank, size_t &checked) + -> void { for (size_t q = 0; q < stage.size(); ++q) { Monomial<32> mono; for (size_t j = 0; j < stage.k_of[q]; ++j) { diff --git a/cpp/tests/sparse_query_tests.cpp b/cpp/tests/sparse_query_tests.cpp index 00e6b7ab..b67752e2 100644 --- a/cpp/tests/sparse_query_tests.cpp +++ b/cpp/tests/sparse_query_tests.cpp @@ -349,8 +349,7 @@ BOOST_AUTO_TEST_CASE(sparse_record_is_exactly_the_gap_code_it_costed) { const size_t w = SQ::push(buf, pos.data(), pos.size(), 1); const size_t gwid = SQ::gap_width(pos.data(), pos.size()); BOOST_TEST(w == SQ::words_of(SQ::gap_bits(pos.size(), gwid)), - "k=" << k << " wrote " << w << " words, costed " - << SQ::words_of(SQ::gap_bits(pos.size(), gwid))); + "k=" << k << " wrote " << w << " words, costed " << SQ::words_of(SQ::gap_bits(pos.size(), gwid))); BOOST_TEST(w <= SQ::words_of(SQ::header_bits_for(pos.size()) + (pos.size() * SQ::kPosBits))); } }