From 249aee85cc31e24cd0ed55a76a2e3f84efb8bddc Mon Sep 17 00:00:00 2001 From: Aaron Miller Date: Thu, 20 Aug 2026 20:59:11 +0100 Subject: [PATCH] =?UTF-8?q?perf(evolution):=20=E2=9A=A1=20carry=20cross-ra?= =?UTF-8?q?nk=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 | 124 ----- 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(+), 285 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 278cc9ff..018d2750 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" namespace monoprop { @@ -95,6 +97,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; @@ -154,6 +158,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. @@ -169,6 +194,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); @@ -180,6 +210,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 @@ -242,6 +277,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. @@ -256,6 +305,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 216eb259..7256c882 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Common.h +++ b/cpp/monoprop/detail/evolution/layer_build/Common.h @@ -98,24 +98,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"); @@ -126,45 +110,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 e507e96d..3ef7ef1c 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" @@ -69,7 +70,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(); } @@ -133,11 +137,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: @@ -185,7 +194,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; } @@ -225,7 +238,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; } @@ -239,7 +252,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 { @@ -248,16 +261,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 { @@ -275,10 +291,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 @@ -296,8 +316,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 @@ -313,6 +339,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. @@ -346,8 +376,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) { @@ -388,7 +420,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; @@ -402,22 +435,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); @@ -434,13 +468,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 { @@ -454,7 +494,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; } @@ -465,24 +508,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]; @@ -492,7 +549,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) { @@ -506,7 +569,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 c5f1e651..1272ee94 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" namespace monoprop::detail { @@ -32,28 +34,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]; @@ -68,22 +100,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; @@ -110,11 +159,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 @@ -133,7 +186,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 e46a7890..88043b33 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" @@ -28,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/graph_encoding/MPGraphEncodingTypes.h" #include "monoprop/detail/mpi/MPIUtils.h" #include "monoprop/detail/operator/InvertedIndex.h" @@ -267,6 +269,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) { @@ -276,31 +312,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 @@ -326,9 +347,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 c094c157..8f86211b 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 @@ -95,6 +97,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); @@ -136,6 +141,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) { @@ -169,6 +201,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); @@ -232,6 +277,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); @@ -250,12 +355,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 @@ -385,6 +513,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 f167fe1f..e19dd33c 100644 --- a/cpp/tests/README.md +++ b/cpp/tests/README.md @@ -70,6 +70,9 @@ name and cannot address suite-nested cases, tests use flat - **`GraphBuildHarness.h`**: direct Layer/MPGraph construction helpers (`core_with_gate`, `layer_with_gate`, `graph_with_gates`) for white-box MPGraph transform tests. +- **`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. @@ -80,17 +83,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 2eb69ab0..00000000 --- a/cpp/tests/fused_query_codec_tests.cpp +++ /dev/null @@ -1,124 +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/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); +}