diff --git a/cpp/monoprop/algebra/AlgebraCommon.h b/cpp/monoprop/algebra/AlgebraCommon.h index 80ef32b9..9b6f8eb5 100644 --- a/cpp/monoprop/algebra/AlgebraCommon.h +++ b/cpp/monoprop/algebra/AlgebraCommon.h @@ -16,6 +16,7 @@ #include #include +#include #include #include #include @@ -24,6 +25,7 @@ #include "monoprop/TypeAliases.h" #include "monoprop/Utilities.h" +#include "monoprop/core/SparseMonomial.h" #include "monoprop/detail/operator/RowAccess.h" namespace monoprop { @@ -96,6 +98,8 @@ auto is_paired(const VecZ &mono) -> bool { return is_paired(indices_to_bitset(mono)); } +// The (k, d) digest form of the predicate above is is_paired(size_t, size_t) in SparseMonomial.h. + template auto is_fully_paired(const VecZ &inds, const Rows &op) -> VecZ { VecZ result; @@ -155,6 +159,27 @@ template return {(first_pair ^ second_pair).count(), active_mono.count(), (first_pair | second_pair).count()}; } +// The same sums from a (k, d) digest; no logical_num_modes because the masking above is inert for a +// well-formed monomial (every set bit at physical position >= 2 * (NumModes - logical_num_modes)). +[[nodiscard]] inline constexpr auto cutoff_sums(size_t k, size_t d) noexcept -> CutoffSums { + return {k - (2 * d), k, k - d}; +} + +// d alone: mode m owns bits (2m, 2m+1) LSb0, so `w & (w >> 1)` masked to even bits counts each +// doubly-occupied mode once. The shift is word-local because a carry would land on odd bit 63. +// Same well-formedness precondition as cutoff_sums(k, d); a monomial built below the active offset by +// hand (majorana_cutoff_tests.cpp:79,101) must keep the bitset overload, which stays the oracle. +template +[[gnu::always_inline]] [[nodiscard]] inline auto paired_mode_count(const Monomial &mono) noexcept -> size_t { + constexpr auto even = even_bits<2 * NumModes, LSb0>(); + size_t d = 0; + for (size_t w = 0; w < Monomial::num_words(); ++w) { + const uint64_t word = mono.word(w); + d += static_cast(std::popcount(word & (word >> 1) & even.word(w))); + } + return d; +} + // Both cutoffs below keep a fully paired monomial (xor_sum == 0) unconditionally: those are the only // terms contributing to an expectation value against a product reference state, so bounding them by // length or support would discard signal. @@ -170,6 +195,11 @@ auto length_cutoff(const Monomial &mono, unsigned int cutoff) -> bool return length_cutoff(mono, cutoff, NumModes); } +// Digest form, on the same precondition as cutoff_sums(k, d). Width-independent, hence not templated. +[[nodiscard]] inline constexpr auto length_cutoff(size_t k, size_t d, unsigned int cutoff) noexcept -> bool { + return length_keeps(k, d, cutoff); +} + template auto support_cutoff(const Monomial &mono, unsigned int cutoff, size_t logical_num_modes) -> bool { const auto sums = cutoff_sums(mono, logical_num_modes); @@ -181,6 +211,11 @@ auto support_cutoff(const Monomial &mono, unsigned int cutoff) -> bool return support_cutoff(mono, cutoff, NumModes); } +// Digest form, on the same precondition as cutoff_sums(k, d). Width-independent, hence not templated. +[[nodiscard]] inline constexpr auto support_cutoff(size_t k, size_t d, unsigned int cutoff) noexcept -> bool { + return support_keeps(k, d, cutoff); +} + namespace detail { template @@ -243,6 +278,29 @@ class CutoffEvaluator { return cutoff_fn_(mono); } + // The decision from the (k, d) digest alone. The emit site gets both out of the partner merge, so + // nothing here reads a bitset. Same precondition as cutoff_sums(k, d): d must have been folded + // without an active mask, which holds for a well-formed monomial. nullopt if opaque. + auto passes_from_digest(size_t k, size_t d) const -> std::optional { + if (length_cutoff_ != nullptr) { + return length_keeps(k, d, length_cutoff_->cutoff); + } + if (support_cutoff_ != nullptr) { + return support_keeps(k, d, support_cutoff_->cutoff); + } + return std::nullopt; + } + + // The same decision when only the dense form is at hand, so d must be folded out of it. + auto passes_from_dense(const Monomial &mono, size_t k) const -> std::optional { + // paired_mode_count has no active_mask, so it agrees with cutoff_sums(mono, L) only above it. + assert(mono.find_first() >= active_bit_offset_() && "monomial has a set bit below its active offset"); + if (length_cutoff_ == nullptr && support_cutoff_ == nullptr) { + return std::nullopt; + } + return passes_from_digest(k, paired_mode_count(mono)); + } + // Upper bound on the set bits (physical slots) a surviving term can carry, so the store can size // its packed inline rows. A length cutoff counts set bits directly; a support cutoff counts // modes/qubits, each spanning two slots, hence the x2. @@ -257,6 +315,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/CMakeLists.txt b/cpp/monoprop/core/CMakeLists.txt index d9faa730..6fb5707c 100644 --- a/cpp/monoprop/core/CMakeLists.txt +++ b/cpp/monoprop/core/CMakeLists.txt @@ -5,4 +5,5 @@ target_sources( TYPE HEADERS FILES "Monomial.h" + "SparseMonomial.h" ) diff --git a/cpp/monoprop/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/CMakeLists.txt b/cpp/monoprop/detail/evolution/layer_build/CMakeLists.txt index eb4c9853..f0c81013 100644 --- a/cpp/monoprop/detail/evolution/layer_build/CMakeLists.txt +++ b/cpp/monoprop/detail/evolution/layer_build/CMakeLists.txt @@ -7,6 +7,9 @@ target_sources( "Common.h" "Engine.h" "FusedApply.h" + "PartnerMerge.h" + "QueryCodec.h" "Resolve.h" "Scan.h" + "SparseQuery.h" ) diff --git a/cpp/monoprop/detail/evolution/layer_build/Common.h b/cpp/monoprop/detail/evolution/layer_build/Common.h index e3c203d3..8777ce24 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Common.h +++ b/cpp/monoprop/detail/evolution/layer_build/Common.h @@ -104,24 +104,8 @@ struct FusedContract { std::vector cross_half; // R>1: one half per cross-rank query (resolver +φ, querier −φ) }; -// Queries ride flat VecZ buffers: kQueryWords elements per query (W monomial words + one ±1 phase word). -// The source index is not in the payload — the resolver answers by position; the querier holds src_idx_r[r][q]. -template -inline constexpr size_t kQueryWords = mpi_detail::kWords + 1; - -// Fused query+value record width (R>1): the plain query record plus one trailing word holding the source's -// pre-cos coeff (v_src, bit-cast from double), so query + value ride a single alltoallv instead of two. -template -inline constexpr size_t kQueryWordsFused = kQueryWords + 1; - -// The unsigned-int intermediate normalizes the ±1 sign bit into a fixed 32-bit pattern so the round-trip -// is exact for any VecZ element width. Edit encode/decode as a pair. -inline auto encode_phase(int phase) -> size_t { - return static_cast(static_cast(phase)); -} -inline auto decode_phase(size_t word) -> int { - return static_cast(static_cast(word)); -} +// Queries ride flat VecZ buffers in one VARIABLE-WIDTH format (SparseQuery): no stride exists, so every +// offset comes from QueryCodec's walk. The source index is not on the wire; the querier holds src_idx_r. // bit_cast, not a conversion, so v_src arrives over the wire bit-identical. static_assert(sizeof(size_t) == sizeof(double), "fused query value word assumes 64-bit VecZ element"); @@ -132,45 +116,4 @@ inline auto decode_value(size_t word) -> double { return std::bit_cast(word); } -template -inline auto query_push(VecZ &buf, const Monomial &mono, int phase) -> void { - mpi_detail::append_monomial_words(mono, buf); - buf.push_back(encode_phase(phase)); -} - -// The mono + phase words occupy the same leading offsets in the plain and fused record, so readers differ -// only in the per-record stride QW (defaulted to the plain width). -template > -inline auto query_read(const VecZ &buf, size_t q, Monomial &mono_out, int &phase_out) -> void { - const size_t base = q * QW; - mono_out = mpi_detail::read_monomial_from_words(buf, base); - phase_out = decode_phase(buf[base + mpi_detail::kWords]); -} - -// No monomial reconstruction: process_responses needs only the phase. -template > -inline auto query_phase(const VecZ &buf, size_t q) -> int { - return decode_phase(buf[q * QW + mpi_detail::kWords]); -} - -template -inline auto query_value(const VecZ &buf, size_t q) -> double { - return decode_value(buf[q * kQueryWordsFused + mpi_detail::kWords + 1]); -} - -// Requires v.size() == q.size()/kQueryWords: exactly one value per query record. -template -inline auto build_fused_query_value(const VecZ &q, const std::vector &v, VecZ &out) -> void { - constexpr size_t W = kQueryWords; - const size_t nq = q.empty() ? 0 : q.size() / W; - out.clear(); - out.reserve(nq * kQueryWordsFused); - for (size_t i = 0; i < nq; ++i) { - out.insert(out.end(), - q.begin() + static_cast(i * W), - q.begin() + static_cast((i + 1) * W)); - out.push_back(encode_value(v[i])); - } -} - } // namespace monoprop::detail diff --git a/cpp/monoprop/detail/evolution/layer_build/Engine.h b/cpp/monoprop/detail/evolution/layer_build/Engine.h index d073f00c..0f03f4af 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Engine.h +++ b/cpp/monoprop/detail/evolution/layer_build/Engine.h @@ -29,6 +29,8 @@ #include "monoprop/algebra/Algebra.h" #include "monoprop/detail/evolution/CutoffContext.h" #include "monoprop/detail/evolution/layer_build/Common.h" +#include "monoprop/detail/evolution/layer_build/PartnerMerge.h" +#include "monoprop/detail/evolution/layer_build/QueryCodec.h" #include "monoprop/detail/evolution/layer_build/Resolve.h" #include "monoprop/detail/evolution/layer_build/Scan.h" #include "monoprop/detail/graph_encoding/MPGraphEncodingStorage.h" @@ -70,7 +72,10 @@ inline auto append_inserted_endpoints(CosMask &cos_all, size_t combined_size, co template struct GraphSink { static constexpr bool wants_values = false; - static constexpr size_t kStride = kQueryWords; + // Named apart: incoming_layout is what this rank RECEIVES, querier_layout its OWN send buffer. They + // coincide here only because GraphSink never fuses -- see ContractSink::querier_layout. + [[nodiscard]] auto incoming_layout() const -> QueryLayout { return {/*fused=*/false}; } + [[nodiscard]] auto querier_layout() const -> QueryLayout { return {/*fused=*/false}; } using Response = TermIndex; static auto init_response() -> Response { return std::numeric_limits::max(); } @@ -134,11 +139,16 @@ struct GraphSink { auto &out = acc[r].out_entries; const size_t base = out.size(); const size_t nq = resp.size(); + const QueryLayout layout = querier_layout(); out.resize(base + nq); + // Forward walk, not indexing by q: a compact query's width depends on its own popcount. + size_t off = 0; for (size_t q = 0; q < nq; ++q) { assert(resp[q] != std::numeric_limits::max() && "resolver must insert absent cross-rank terms"); - out[base + q] = {srcs[q], query_phase(qbuf, q)}; + out[base + q] = {srcs[q], QueryCodec::phase_at(qbuf, off)}; + off = QueryCodec::next_off(qbuf, layout, off); } + assert(off == qbuf.size() && "querier buffer does not hold exactly one query per response"); } // Drains the per-rank accumulators into the LayerCore's sin_send/sin_recv lists (layout derivation: @@ -186,7 +196,11 @@ struct GraphSink { template struct ContractSink { static constexpr bool wants_values = true; - static constexpr size_t kStride = kQueryWordsFused; + // This rank RECEIVES fused (query+value) records, but the buffer on_response_block is handed is its + // own queries_r, which is PLAIN (build_fused writes the fused form into combined_qv_). Reading the + // phase with the wrong layout takes a neighbouring record's, which is a silent coefficient sign flip. + [[nodiscard]] auto incoming_layout() const -> QueryLayout { return {/*fused=*/true}; } + [[nodiscard]] auto querier_layout() const -> QueryLayout { return {/*fused=*/false}; } using Response = double; static auto init_response() -> Response { return 0.0; } @@ -226,7 +240,7 @@ struct ContractSink { -> std::vector & { scratch.resize(queries.size()); for (size_t r = 0; r < queries.size(); ++r) { - build_fused_query_value(queries[r], vals[r], scratch[r]); + QueryCodec::build_fused(queries[r], vals[r], scratch[r]); } return scratch; } @@ -240,7 +254,7 @@ struct ContractSink { } auto on_resolved(size_t g, size_t s, - size_t q, + size_t /*q*/, size_t ip, const IncomingProbe &pr, const std::vector &incoming) -> Response { @@ -249,16 +263,19 @@ struct ContractSink { v_tgt = fused_scale ? op_coeffs[ip] * inv_cos : op_coeffs[ip]; } else if (schrodinger) { - v_tgt = - is_paired(pr.mono[g]) ? algebra_state_phase(basis, pr.mono[g], state_mask_) : 0.0; + // Through the probe's accessors: it holds position lists, and mono_at builds a bitset only + // for the fully paired minority that is_paired_at admits. + v_tgt = pr.is_paired_at(g) ? algebra_state_phase(basis, pr.mono_at(g), state_mask_) : 0.0; } else { v_tgt = 0.0; // Heisenberg fresh insert } - fc.cross_half[cross_base_ + g] = HalfRotationRec{ip, - query_value(incoming[s], q), - static_cast(pr.phase_of[g]), - /*is_insert=*/ip >= pr.base}; + // pr.off_of[g], not q: under the compact record a query ordinal does not name a buffer position. + fc.cross_half[cross_base_ + g] = + HalfRotationRec{ip, + QueryCodec::value_at(incoming[s], incoming_layout(), pr.off_of[g]), + static_cast(pr.phase_of[g]), + /*is_insert=*/ip >= pr.base}; return v_tgt; } auto process_reserve(const std::vector> &inc_r, size_t rank_count, size_t my_rank_) -> void { @@ -276,10 +293,14 @@ struct ContractSink { const std::vector &srcs, const VecZ &qbuf) -> void { const size_t nq = rval.size(); + const QueryLayout layout = querier_layout(); + size_t off = 0; for (size_t q = 0; q < nq; ++q) { - const auto nphase = static_cast(-query_phase(qbuf, q)); + const auto nphase = static_cast(-QueryCodec::phase_at(qbuf, off)); fc.cross_half.push_back(HalfRotationRec{srcs[q], rval[q], nphase, /*is_insert=*/false}); + off = QueryCodec::next_off(qbuf, layout, off); } + assert(off == qbuf.size() && "querier buffer does not hold exactly one query per response"); } // No LayerCore in the fused path → nullptr. Two-pass fused (k>0 / cos==0 fallback) appends inserted @@ -297,8 +318,14 @@ struct ContractSink { // Owns build_layer's machinery over a compile-time Sink policy. combined_size = the pre-layer operator size. template struct LayerBuildEngine { + // The store's position type, narrower than the wire's below 129 modes; decoded straight into. + using RowPosT = typename OperatorIndex::PosT; + + // A miss keeps its decoded positions (pos_at indexes deferred_pos_flat_) and the probe's hash. struct DeferredSelfMiss { - Monomial mono; + size_t pos_at; + uint32_t k; + uint32_t hash; size_t src; int phase; double v_src = 0.0; // ContractSink only: op_pre[src] captured at scan emit; 0 for GraphSink @@ -314,6 +341,11 @@ 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_; + // This pass's self-owned queries as positions, straight from the scan: never encoded, so the resolve + // below has nothing to decode. Parallel to src_idx_r[my_rank]. + SelfQueryStage self_stage_; // Scan-captured v_src per query (ContractSink only via Sink::wants_values; empty for GraphSink). std::vector> src_val_r; // Fused query+value send scratch (ContractSink, R>1): shared by a gate's two exchange passes. @@ -341,15 +373,16 @@ struct LayerBuildEngine { // Resolve this rank's own query stream inline, then clear it so the alltoallv never sends to self. auto resolve_self_queries(bool is_leader_pass) -> void { - VecZ &lq = queries_r[my_rank]; std::vector &ls = src_idx_r[my_rank]; std::vector *lv = nullptr; if constexpr (Sink::wants_values) { lv = &src_val_r[my_rank]; } - const size_t nq = lq.empty() ? 0 : lq.size() / kQueryWords; - resolve_range_(lq, ls, lv, 0, nq, is_leader_pass); - lq.clear(); + // The scan routes a self-owned partner to the stage, never to the wire buffer. + assert(queries_r[my_rank].empty() && "a self-owned query was encoded instead of staged"); + assert(ls.size() == self_stage_.size() && "the self stage does not hold exactly one query per source"); + resolve_range_(ls, lv, is_leader_pass); + self_stage_.clear(); ls.clear(); if constexpr (Sink::wants_values) { src_val_r[my_rank].clear(); @@ -365,9 +398,11 @@ struct LayerBuildEngine { auto run_exchange(bool is_leader_pass, std::vector &&queries, std::vector> &&src_idx, - std::vector> &&src_val) -> void { + std::vector> &&src_val, + SelfQueryStage &&self_stage) -> void { queries_r = std::move(queries); src_idx_r = std::move(src_idx); + self_stage_ = std::move(self_stage); // src_val is empty unless Sink::wants_values, so the move is a no-op under GraphSink. src_val_r = std::move(src_val); if (!is_leader_pass && R > 1) { @@ -389,7 +424,8 @@ struct LayerBuildEngine { // Followers a leader already matched must not be re-resolved over the wire, so compact them out. auto drop_matched_cross_rank_followers() -> void { - constexpr size_t W = kQueryWords; + using QC = QueryCodec; + const QueryLayout layout = sink.querier_layout(); for (size_t r = 0; r < R; ++r) { if (r == my_rank) { continue; @@ -403,22 +439,23 @@ struct LayerBuildEngine { } const size_t nq = s.size(); size_t kept = 0; + // Two cursors, since a dropped query has no fixed width; order is the accumulation order. + size_t src_off = 0; + size_t dst_off = 0; for (size_t k = 0; k < nq; ++k) { - if (matched.is_marked(s[k])) { - continue; - } - if (kept != k) { - std::copy(q.begin() + static_cast(k * W), - q.begin() + static_cast((k + 1) * W), - q.begin() + static_cast(kept * W)); - } - s[kept] = s[k]; - if (v != nullptr) { - (*v)[kept] = (*v)[k]; + const size_t next = QC::next_off(q, layout, src_off); + if (!matched.is_marked(s[k])) { + dst_off += QC::move_query(q, layout, src_off, dst_off); + s[kept] = s[k]; + if (v != nullptr) { + (*v)[kept] = (*v)[k]; + } + ++kept; } - ++kept; + src_off = next; } - q.resize(kept * W); + assert(src_off == q.size() && "follower compaction did not consume the whole query buffer"); + q.resize(dst_off); s.resize(kept); if (v != nullptr) { v->resize(kept); @@ -435,13 +472,19 @@ struct LayerBuildEngine { if (n_miss == 0) { return; } - auto key_at = [&](size_t k) -> const Monomial & { return deferred_self_misses[k].mono; }; sink.prepare_deferred(n_miss); - insert_absent_terms(local_op, n_miss, key_at, [&](size_t k, size_t base) { + // insert_absent_terms' three steps without its dense round trips, on the same ordering contract: + // miss k lands at base+k, in leader-then-follower order. + // insert_absent_terms is the dense reference this path is differentially tested against + // (sparse_resolve_tests.cpp), so it must not be deleted for having no library caller. + const size_t base = local_op.store->grow_rows_geometric(n_miss); + for (size_t k = 0; k < n_miss; ++k) { const auto &m = deferred_self_misses[k]; - assign_row(*local_op.store, base + k, m.mono); + local_op.store->set_positions(base + k, deferred_pos_flat_.data() + m.pos_at, m.k); sink.emit_deferred(k, base + k, m.src, m.phase, m.v_src); - }); + } + local_op.store->bulk_insert_hashed(n_miss, base, [&](size_t j) { return deferred_self_misses[j].hash; }); + local_op.reindex_after_growth(base, n_miss); } auto finish(CosMask &&cos_all, CosMask *out_cos = nullptr) -> std::shared_ptr { @@ -455,7 +498,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; } @@ -463,19 +509,20 @@ struct LayerBuildEngine { // Batched self-resolve over the index's group-prefetch find_batch; hits/misses are emitted to the sink // in query order. `lv` is the per-query v_src array parallel to `ls` (read only when Sink::wants_values). static constexpr size_t kResolveBatch = 64; - auto resolve_range_(VecZ &lq, - std::vector &ls, - [[maybe_unused]] std::vector *lv, - size_t lo, - size_t hi, - bool is_leader_pass) -> void { + auto resolve_range_(std::vector &ls, [[maybe_unused]] std::vector *lv, bool is_leader_pass) + -> void { const size_t op_size = local_op.store->size(); - std::array, kResolveBatch> keys; + // Gathered per batch because a matched follower is skipped; the offsets stay ABSOLUTE into the + // stage's pos_flat, so find_batch_positions reads it in place and nothing is copied. + std::array pos_off; + std::array k_of; + std::array hashes; std::array phases; std::array srcs; std::array vals; std::array found; - size_t q = lo; + const size_t hi = self_stage_.size(); + size_t q = 0; while (q < hi) { size_t m = 0; for (; q < hi && m < kResolveBatch; ++q) { @@ -483,7 +530,9 @@ struct LayerBuildEngine { if (!is_leader_pass && matched.is_marked(src)) { continue; // follower already matched by a leader → not an independent rotation } - query_read(lq, q, keys[m], phases[m]); + pos_off[m] = self_stage_.pos_off[q]; + k_of[m] = self_stage_.k_of[q]; + phases[m] = self_stage_.phase_of[q]; srcs[m] = src; if constexpr (Sink::wants_values) { vals[m] = (*lv)[q]; @@ -493,7 +542,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_stage_.pos_flat.data(), + pos_off.data(), + k_of.data(), + m, + found.data(), + hashes.data()); for (size_t j = 0; j < m; ++j) { double v_src = 0.0; if constexpr (Sink::wants_values) { @@ -508,7 +563,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}); + // The stage dies with this pass and the misses are flushed after both, so copy now. + const size_t at = deferred_pos_flat_.size(); + const auto *const first = self_stage_.pos_flat.data() + pos_off[j]; + deferred_pos_flat_.insert(deferred_pos_flat_.end(), first, first + k_of[j]); + deferred_self_misses.push_back({at, k_of[j], hashes[j], srcs[j], phases[j], v_src}); } } } @@ -560,7 +619,7 @@ auto build_layer(MPOperator &local_op, } assert(fused_scale_coeffs == nullptr || (local_coeffs && &local_coeffs->get() == fused_scale_coeffs)); - FusedScanResult fused = [&] { + FusedScanResult fused = [&] { double *const sweep_ptr = fused_scale ? fused_scale_coeffs->data() : nullptr; return with_algebra(basis, [&]() { return fused_find_and_collect(local_op, @@ -602,11 +661,13 @@ auto build_layer(MPOperator &local_op, eng.run_exchange(/*is_leader_pass=*/true, std::move(fused.leader_queries), std::move(fused.leader_src), - std::move(fused.leader_val)); + std::move(fused.leader_val), + std::move(fused.leader_self)); eng.run_exchange(/*is_leader_pass=*/false, std::move(fused.follower_queries), std::move(fused.follower_src), - std::move(fused.follower_val)); + std::move(fused.follower_val), + std::move(fused.follower_self)); return eng.finish(std::move(cos_all), out_cos); }; diff --git a/cpp/monoprop/detail/evolution/layer_build/PartnerMerge.h b/cpp/monoprop/detail/evolution/layer_build/PartnerMerge.h new file mode 100644 index 00000000..8b5bc668 --- /dev/null +++ b/cpp/monoprop/detail/evolution/layer_build/PartnerMerge.h @@ -0,0 +1,172 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +// M⊕G as ascending positions. A slot in both M and G cancels (m_p m_p = 1), so the partner is the +// symmetric difference of two ascending position lists, and one merge yields its positions, `overlap` +// and `d` (modes carrying BOTH Majoranas) together -- the (k, d) digest the structural cutoff wants, +// with no second sweep over the dense form and no walk back out of it. + +#include +#include +#include +#include +#include + +#include "monoprop/TypeAliases.h" +#include "monoprop/detail/operator/OperatorIndex.h" + +namespace monoprop::detail { + +// Both inputs must be strictly ascending. The output is their symmetric difference, so it is bounded +// by the universe the positions are drawn from -- 2*NumModes here -- and ka + kb is only the bound that +// ignores cancellation. Returns the merged count. GenT is separate from PosT because the generator's +// positions are the wire's width, not the store's. +template +[[gnu::always_inline]] inline auto merge_partner_positions(const PosT *a, + size_t ka, + const GenT *b, + size_t kb, + PosT *out, + size_t &overlap_out, + size_t &d_out) noexcept -> size_t { + size_t i = 0; + size_t j = 0; + size_t n = 0; + size_t overlap = 0; + size_t d = 0; + // Seeded ODD, so the (prev % 2 == 0) test cannot fire on the first emit and the loops need no + // n != 0 guard; 1 is not a reachable `prev + 1` either, since prev would have to be 0 and even. + size_t prev = 1; + // Ascending output, so a doubly-occupied mode is an even position immediately followed by its + // successor -- the same count paired_mode_count folds out of the bitset. Written out three times + // rather than through a lambda: callgrind measured 15,279,191 CALLS to that lambda at 18 + // instructions each (275.1M, a third of this port's whole delta), because GCC declined to inline a + // closure capturing five locals by reference into three call sites. + while (i < ka && j < kb) { + const size_t pa = static_cast(a[i]); + const size_t pb = static_cast(b[j]); + if (pa == pb) { + ++overlap; + ++i; + ++j; + continue; + } + const size_t p = pa < pb ? pa : pb; + i += static_cast(pa < pb); + j += static_cast(pb < pa); + d += static_cast((prev % 2 == 0) && p == prev + 1); + out[n++] = static_cast(p); + prev = p; + } + for (; i < ka; ++i) { + const size_t p = static_cast(a[i]); + d += static_cast((prev % 2 == 0) && p == prev + 1); + out[n++] = static_cast(p); + prev = p; + } + for (; j < kb; ++j) { + const size_t p = static_cast(b[j]); + d += static_cast((prev % 2 == 0) && p == prev + 1); + out[n++] = static_cast(p); + prev = p; + } + overlap_out = overlap; + d_out = d; + return n; +} + +// Self-owned queries never reach a wire, so they are staged as positions rather than encoded records: +// OperatorIndex's find_batch_positions and set_positions both take exactly this shape, so the resolve +// path consumes the stage with no transformation and the codec is not on the self leg at all. +template +struct SelfQueryStage { + using PosT = typename OperatorIndex::PosT; + + // SIZED, not filled: the vectors carry the capacity and n_/pos_n_ carry the logical length, so a + // push writes rather than appends. Read them through size() and the data pointers only. + // + // DefaultInitVector, the allocator Resolve.h already uses for exactly this: a plain vector's resize + // VALUE-initialises, so sizing pos_flat ahead would memset every byte a push is about to overwrite + // -- trading the append cost for a per-gate zero-fill instead of removing it. + DefaultInitVector pos_flat; // ascending positions, concatenated in push order + DefaultInitVector pos_off; // query -> absolute offset into pos_flat + DefaultInitVector k_of; + DefaultInitVector phase_of; // emit_phase is ternary, so a byte is the whole range + + [[nodiscard]] auto size() const -> size_t { return n_; } + [[nodiscard]] auto positions() const -> size_t { return pos_n_; } + + auto clear() -> void { + n_ = 0; + pos_n_ = 0; + } + + auto reserve(size_t n_queries, size_t positions_per_query) -> void { + if (pos_off.size() < n_queries) { + pos_off.resize(n_queries); + k_of.resize(n_queries); + phase_of.resize(n_queries); + } + if (pos_flat.size() < n_queries * positions_per_query) { + pos_flat.resize(n_queries * positions_per_query); + } + } + + // Four preallocated writes, not four container appends. Callgrind on the pauli cell put + // vector::_M_range_insert at 110.6M instructions and vector::emplace_back at + // 52.2M -- 18% of this port's whole instruction delta -- for a push whose capacity is already + // reserved. `insert` cannot know that, so it re-derives the grow path per query; grow_() is the + // one place that checks, and it runs once per capacity doubling instead of once per push. + auto push(const PosT *pos, size_t k, int phase) -> void { + assert(phase >= -1 && phase <= 1 && "emit_phase is ternary: rotation_sign, or REAL_PARTS entry"); + const size_t n = n_; + const size_t at = pos_n_; + if (n == pos_off.size() || at + k > pos_flat.size()) { + grow_(k); + } + pos_off[n] = at; + // An explicit loop, not std::copy_n: k averages ~5 bytes here and copy_n compiles to a memcpy + // CALL, which callgrind counted 2.54M extra times for a copy smaller than its own prologue. + PosT *dst = pos_flat.data() + at; + for (size_t j = 0; j < k; ++j) { + dst[j] = pos[j]; + } + k_of[n] = static_cast(k); + phase_of[n] = static_cast(phase); + n_ = n + 1; + pos_n_ = at + k; + } + +private: + size_t n_ = 0; // queries pushed + size_t pos_n_ = 0; // positions pushed + + // Amortised doubling, and pos_flat grows by the larger of a double and what this push needs, so a + // single wide term cannot leave it short. + [[gnu::noinline]] auto grow_(size_t k) -> void { + if (n_ == pos_off.size()) { + const size_t want = (pos_off.size() * 2) + 64; + pos_off.resize(want); + k_of.resize(want); + phase_of.resize(want); + } + if (pos_n_ + k > pos_flat.size()) { + pos_flat.resize(std::max((pos_flat.size() * 2) + 256, pos_n_ + k)); + } + } +}; + +} // namespace monoprop::detail diff --git a/cpp/monoprop/detail/evolution/layer_build/QueryCodec.h b/cpp/monoprop/detail/evolution/layer_build/QueryCodec.h new file mode 100644 index 00000000..4a6e7411 --- /dev/null +++ b/cpp/monoprop/detail/evolution/layer_build/QueryCodec.h @@ -0,0 +1,145 @@ +// 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 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); + } + + // From ascending positions, which is what the partner merge hands the emit site; the dense overload + // above is for callers that hold only a bitset. + template + static auto push_positions(VecZ &buf, const PosU *pos, size_t k, int phase) -> size_t { + return CQ::push(buf, pos, k, phase); + } + + // Identical in both formats: the value is one bit_cast word after the query's words. + static auto push_value(VecZ &buf, double v) -> void { buf.push_back(encode_value(v)); } + + // Inflates the record back into a dense Monomial, for callers that cannot consume positions. + static auto read_mono(const VecZ &buf, size_t off, Monomial &mono_out, int &phase_out) -> void { + (void)CQ::read_mono(buf, off, mono_out, phase_out); + } + + // The query's popcount, straight out of the record's header field. + [[nodiscard]] static auto k_at(const VecZ &buf, size_t off) -> size_t { return CQ::k_at(buf, off); } + + // Positions plus phase, into the CALLER's element type: the store's PosT is narrower below 129 modes. + template + static auto read_positions(const VecZ &buf, QueryLayout layout, size_t off, OutT *out, int &phase_out) -> size_t { + phase_out = CQ::phase_at(buf, off); + return CQ::read_positions(buf, off, out) + (layout.fused ? 1U : 0U); + } + + [[nodiscard]] static auto phase_at(const VecZ &buf, size_t off) -> int { return CQ::phase_at(buf, off); } + + [[nodiscard]] static auto value_at(const VecZ &buf, [[maybe_unused]] QueryLayout layout, size_t off) -> double { + assert(layout.fused && "there is no value word in a plain query buffer"); + return decode_value(buf[off + query_words(buf, off)]); + } + + // The number of QUERIES. Genuinely a walk: records vary in width, so there is no stride to divide by. + [[nodiscard]] static auto count_queries(const VecZ &buf, QueryLayout layout) -> size_t { + size_t off = 0; + size_t n = 0; + while (off < buf.size()) { + off = next_off(buf, layout, off); + ++n; + } + assert(off == buf.size() && "a compact query ran past the end of the buffer"); + return n; + } + + // Interleave a plain query stream with its parallel v_src array; a size mismatch shifts every coeff. + static auto build_fused(const VecZ &queries, const std::vector &vals, VecZ &out) -> void { + out.clear(); + out.reserve(queries.size() + vals.size()); + size_t off = 0; + size_t i = 0; + while (off < queries.size()) { + const size_t n = query_words(queries, off); + out.insert(out.end(), + queries.begin() + static_cast(off), + queries.begin() + static_cast(off + n)); + assert(i < vals.size() && "fused build needs exactly one value per query"); + out.push_back(encode_value(vals[i])); + off += n; + ++i; + } + assert(i == vals.size() && "fused build needs exactly one value per query"); + } + + // Copy the query at `src_off` (with its value word, if fused) to `dst_off`; returns words written. + static auto move_query(VecZ &buf, QueryLayout layout, size_t src_off, size_t dst_off) -> size_t { + const size_t n = query_words(buf, src_off) + (layout.fused ? 1U : 0U); + if (src_off != dst_off) { + assert(dst_off < src_off && "compaction only ever moves a query earlier"); + std::copy(buf.begin() + static_cast(src_off), + buf.begin() + static_cast(src_off + n), + buf.begin() + static_cast(dst_off)); + } + return n; + } +}; + +} // namespace monoprop::detail diff --git a/cpp/monoprop/detail/evolution/layer_build/Resolve.h b/cpp/monoprop/detail/evolution/layer_build/Resolve.h index ef90d91f..0c49706f 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Resolve.h +++ b/cpp/monoprop/detail/evolution/layer_build/Resolve.h @@ -14,6 +14,7 @@ #pragma once +#include #include #include #include @@ -22,6 +23,7 @@ #include "monoprop/algebra/Algebra.h" #include "monoprop/detail/evolution/CutoffContext.h" #include "monoprop/detail/evolution/layer_build/Common.h" +#include "monoprop/detail/evolution/layer_build/QueryCodec.h" #include "monoprop/detail/operator/MPOperator.h" #include "monoprop/detail/operator/RowAccess.h" @@ -33,28 +35,58 @@ namespace monoprop::detail { // pairwise distinct ⇒ misses distinct and absent. template struct IncomingProbe { - std::vector goff; // rank_count+1 flat offsets: g = goff[s] + q - DefaultInitVector sender_of; // g → sender rank - DefaultInitVector> mono; // g → deserialized query monomial - DefaultInitVector phase_of; // g → query phase - DefaultInitVector idx_of; // g → resolved index (hit: < base; miss: base+j) - std::vector miss_g; // j → the g that became miss j (Phase 4 reads mono[miss_g[j]]) - size_t base = 0; // op size before the miss inserts (the miss-index base) + // The STORE's position width, not the wire's: these positions exist to become rows. + using PosT = typename OperatorIndex::PosT; + + std::vector goff; // rank_count+1 flat offsets: g = goff[s] + q + DefaultInitVector sender_of; // g → sender rank + DefaultInitVector phase_of; // g → query phase + // g → WORD offset of that query inside incoming[sender_of[g]]; a query ordinal names no position. + DefaultInitVector off_of; + DefaultInitVector idx_of; // g → resolved index (hit: < base; miss: base+j) + std::vector miss_g; // j → the g that became miss j (Phase 4 reads the key of miss_g[j]) + size_t base = 0; // op size before the miss inserts (the miss-index base) size_t nq_total = 0; + + // The queries as they arrived, flat: query g owns pos_flat[pos_off[g] .. pos_off[g] + k_of[g]). + DefaultInitVector pos_flat; + DefaultInitVector pos_off; + DefaultInitVector k_of; + // g → fold_hash of the query key, folded by the probe and reused by the insert. + DefaultInitVector hash_of; + + // BUILDS a bitset, so cold consumers only -- the fully paired minority, never anything per-term. + [[nodiscard]] auto mono_at(size_t g) const -> Monomial { + Monomial m; + const PosT *p = pos_flat.data() + pos_off[g]; + for (size_t j = 0; j < k_of[g]; ++j) { + m.set(static_cast(p[j])); + } + return m; + } + + // is_paired from the positions' (k, d) digest, no bitset built. + [[nodiscard]] auto is_paired_at(size_t g) const -> bool { + const PosT *p = pos_flat.data() + pos_off[g]; + const size_t k = k_of[g]; + return monoprop::is_paired(k, QueryCodec::pair_count(p, k)); + } }; -// Phases 1-2, read-only w.r.t. operator contents. QW = per-record stride: the plain query width, or -// kQueryWordsFused for the fused resolver. The caller runs Phase 3, then insert_incoming_misses. -template > +// Phases 1-2, read-only w.r.t. operator contents. `layout` describes the records this rank RECEIVES: +// fused for the ContractSink resolver, plain for GraphSink. The caller runs Phase 3, then +// insert_incoming_misses. Counts and offsets come from the decode walk; there is no record stride. +template auto probe_incoming_queries(const std::vector &incoming, // serialized, one VecZ per sender MPOperator &op, - size_t rank_count) -> IncomingProbe { - constexpr size_t W = QW; + size_t rank_count, + QueryLayout layout) -> IncomingProbe { + using QC = QueryCodec; IncomingProbe pr; pr.goff.assign(rank_count + 1, 0); for (size_t s = 0; s < rank_count; ++s) { - const size_t nq = incoming[s].empty() ? 0 : incoming[s].size() / W; + const size_t nq = QC::count_queries(incoming[s], layout); pr.goff[s + 1] = pr.goff[s] + nq; } pr.nq_total = pr.goff[rank_count]; @@ -69,22 +101,39 @@ auto probe_incoming_queries(const std::vector &incoming, // serialized, on static_cast(s)); } - // Phase 1 (read-only): deserialize, then probe with the group-prefetch batch find. - pr.mono.resize(pr.nq_total); + // Phase 1 (read-only): deserialize, then probe with the group-prefetch batch find. One walk per sender. pr.phase_of.resize(pr.nq_total); + pr.off_of.resize(pr.nq_total); pr.idx_of.resize(pr.nq_total); - for (size_t g = 0; g < pr.nq_total; ++g) { - const size_t s = pr.sender_of[g]; - const size_t q = g - pr.goff[s]; - Monomial m; - int ph = 0; - query_read(incoming[s], q, m, ph); - pr.mono[g] = m; - pr.phase_of[g] = ph; + pr.pos_off.resize(pr.nq_total); + pr.k_of.resize(pr.nq_total); + pr.hash_of.resize(pr.nq_total); + pr.pos_flat.clear(); + // A hint only: the measured mean is 5.33 positions, so this is one allocation but for an outlier. + pr.pos_flat.reserve(pr.nq_total * QueryCodec::kReservePositionsPerQuery); + for (size_t s = 0; s < rank_count; ++s) { + size_t off = 0; + for (size_t g = pr.goff[s]; g < pr.goff[s + 1]; ++g) { + int ph = 0; + const size_t k = QC::k_at(incoming[s], off); + const size_t at = pr.pos_flat.size(); + pr.pos_flat.resize(at + k); // default-init grow: read_positions writes every element + pr.pos_off[g] = at; + pr.k_of[g] = static_cast(k); + pr.off_of[g] = off; + off = QC::read_positions(incoming[s], layout, off, pr.pos_flat.data() + at, ph); + pr.phase_of[g] = ph; + } + assert(off == incoming[s].size() && "the query walk did not consume the sender's whole buffer"); } { const size_t op_size = op.store->size(); - op.store->find_batch(pr.mono.data(), pr.nq_total, pr.idx_of.data()); + op.store->find_batch_positions(pr.pos_flat.data(), + pr.pos_off.data(), + pr.k_of.data(), + pr.nq_total, + pr.idx_of.data(), + pr.hash_of.data()); for (size_t g = 0; g < pr.nq_total; ++g) { if (pr.idx_of[g] >= op_size) { // kNotFound is size_t max → also lands here pr.idx_of[g] = kMissingIndex; @@ -111,11 +160,15 @@ auto insert_incoming_misses(MPOperator &op, const IncomingProbe( - op, - n_miss, - [&](size_t j) -> const Monomial & { return pr.mono[pr.miss_g[j]]; }, - [&](size_t j, size_t base) { assign_row(*op.store, base + j, pr.mono[pr.miss_g[j]]); }); + // insert_absent_terms' three steps without its two dense round-trips, and on the same ordering + // contract, which is what matters: slot j lands at base+j, in miss order = (sender, record) order. + const size_t base = op.store->grow_rows_geometric(n_miss); + for (size_t j = 0; j < n_miss; ++j) { + const size_t g = pr.miss_g[j]; + op.store->set_positions(base + j, pr.pos_flat.data() + pr.pos_off[g], pr.k_of[g]); + } + op.store->bulk_insert_hashed(n_miss, base, [&](size_t j) { return pr.hash_of[pr.miss_g[j]]; }); + op.reindex_after_growth(base, n_miss); } // resolve_incoming / process_responses are the picture-independent cross-rank exchange skeletons; what @@ -134,7 +187,8 @@ auto resolve_incoming(const std::vector &incoming, // serialized, one VecZ size_t combined_size, // pre-layer op size: bounds the matched set Sink &sink) -> std::vector> { using Resp = typename Sink::Response; - const IncomingProbe pr = probe_incoming_queries(incoming, op, rank_count); + const IncomingProbe pr = + probe_incoming_queries(incoming, op, rank_count, sink.incoming_layout()); std::vector> responses(rank_count); for (size_t s = 0; s < rank_count; ++s) { responses[s].assign(pr.goff[s + 1] - pr.goff[s], Sink::init_response()); diff --git a/cpp/monoprop/detail/evolution/layer_build/Scan.h b/cpp/monoprop/detail/evolution/layer_build/Scan.h index d5b8a77e..cce972b0 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Scan.h +++ b/cpp/monoprop/detail/evolution/layer_build/Scan.h @@ -21,6 +21,7 @@ #include #include #include +#include #include #include "monoprop/TypeAliases.h" @@ -29,6 +30,8 @@ #include "monoprop/core/Monomial.h" #include "monoprop/detail/evolution/CutoffContext.h" #include "monoprop/detail/evolution/layer_build/Common.h" +#include "monoprop/detail/evolution/layer_build/PartnerMerge.h" +#include "monoprop/detail/evolution/layer_build/QueryCodec.h" #include "monoprop/detail/graph_encoding/MPGraphEncodingTypes.h" #include "monoprop/detail/mpi/MPIUtils.h" #include "monoprop/detail/operator/InvertedIndex.h" @@ -157,23 +160,65 @@ inline auto rotation_dynamic_gate(std::optional only_rotate_len_k, return true; } +// M⊕G as the emit site needs it. The dense form is unavoidable -- the owner hash folds every word and +// the basis sign reads the source bitset -- so the merge below runs ALONGSIDE it, not instead of it, +// and supplies k, d and the positions without a second sweep (paired_mode_count) or a third +// (push_mono's walk). +template +struct PartnerProduct { + Monomial new_mono; + size_t k = 0; // popcount(M⊕G) + size_t d = 0; // modes of M⊕G carrying BOTH Majoranas + size_t overlap = 0; // slots in both M and G, which cancel + int phase_factor = 0; +}; + // phase_factor is the basis-specific sign only: Majorana interleave_phase, still to be folded with -// hermitian_phase at emit; Pauli pauli_rotation_sign, already rotation-ready. -template +// hermitian_phase at emit; Pauli pauli_rotation_sign, already rotation-ready. `out_pos` receives the +// partner's ascending positions and needs capacity 2*NumModes; a spilled source row has no position +// array, so that case alone walks the dense partner back out. +// +// BOTH SHAPES WERE MEASURED, on the pauli cell over 2,455,950 emit calls (callgrind, jobs cg-sym4 and +// cg-sym5). Walking the dense partner instead -- find_first/find_next for the positions and the same +// running pairing test -- costs 229.1M MORE instructions than this merge, because the walk is a serial +// dependence chain through find_next where the merge streams two ascending arrays. The intuition that +// the merge is redundant work on top of a bitset that exists anyway is wrong: it is cheaper than +// reading that bitset back out. Do not replace it with the walk again. +template [[gnu::always_inline]] inline auto emit_term_products(const OperatorIndex &ham, size_t i, const typename A::GenContext &ctx, - Monomial &new_mono, - size_t &overlap, - int &phase_factor) -> void { - Monomial mono; - ham.for_each_position(i, [&](size_t pos) { mono.set(pos); }); + const GenT *gen_pos, + size_t gen_pop, + PosT *out_pos) -> PartnerProduct { const Monomial &gen = A::generator(ctx); - new_mono = mono ^ gen; - overlap = mono.count_and(gen); - phase_factor = A::rotation_sign(ctx, mono, new_mono); + PartnerProduct out; + Monomial mono; + if (const auto src = ham.row_positions(i); src.inlined()) { + out.k = merge_partner_positions(src.pos, src.count, gen_pos, gen_pop, out_pos, out.overlap, out.d); + for (size_t j = 0; j < src.count; ++j) { + mono.set(static_cast(src.pos[j])); + } + out.new_mono = mono ^ gen; + } + else { + ham.for_each_position(i, [&](size_t pos) { mono.set(pos); }); + out.new_mono = mono ^ gen; + out.overlap = mono.count_and(gen); + size_t prev = 0; + for (size_t b = out.new_mono.find_first(); b < out.new_mono.size(); b = out.new_mono.find_next(b)) { + if (out.k != 0 && (prev % 2 == 0) && b == prev + 1) { + ++out.d; + } + out_pos[out.k++] = static_cast(b); + prev = b; + } + } + out.phase_factor = A::rotation_sign(ctx, mono, out.new_mono); + return out; } +template struct FusedScanResult { std::vector cos_blocks; // ascending, disjoint, chunk order std::vector leader_queries; // size R: serialized leader queries per owner rank @@ -184,6 +229,11 @@ struct FusedScanResult { // leader_src / follower_src. Empty when capture_values is false. std::vector> leader_val; std::vector> follower_val; + // Self-owned queries, staged as positions instead of encoded into leader_queries[my_rank]: they are + // resolved inline and never reach a wire, so the codec is not on this leg. Order matches + // leader_src[my_rank] / follower_src[my_rank], which is the accumulation order. + SelfQueryStage leader_self; + SelfQueryStage follower_self; }; // Classify, cut off and emit in one pass over the anticommuting terms. Queries go to the owner of @@ -202,12 +252,12 @@ auto fused_find_and_collect(const MPOperator &op, size_t my_rank, bool capture_values = false, double *fused_scale_coeffs = nullptr, - double fused_scale_cos = 1.0) -> FusedScanResult { + double fused_scale_cos = 1.0) -> FusedScanResult { validate_only_rotate_len_k_(only_rotate_len_k, 2 * NumModes); const size_t gen_pop = gen.count(); const auto ectx = A::make_gen_context(gen); - FusedScanResult res; + FusedScanResult res; res.leader_queries.assign(rank_count, VecZ{}); res.leader_src.assign(rank_count, std::vector{}); res.follower_queries.assign(rank_count, VecZ{}); @@ -268,40 +318,66 @@ auto fused_find_and_collect(const MPOperator &op, auto &fs = res.follower_src; auto &fv = res.follower_val; + const OperatorIndex &ham = *op.store; + using RowPosT = typename OperatorIndex::PosT; + + // The generator's positions, once per gate: the merge's second input. + std::vector gen_pos; + gen_pos.reserve(gen_pop); + for (size_t b = gen.find_first(); b < gen.size(); b = gen.find_next(b)) { + gen_pos.push_back(static_cast(b)); + } + // 2*NumModes is the true bound: the partner's positions are distinct and below it. NOT + // thread_local: every access to one from a shared library goes through __tls_get_addr, which + // callgrind measured at 54.7M instructions on the pauli cell -- 6% of this port's delta -- to + // save one allocation per gate. + std::vector pbuf(2 * NumModes); + + // Everything after a term survives the structural cutoff. Self-owned partners are staged as + // positions; only a remote owner's partner is encoded. + auto push = [&](const Monomial &dense, + const RowPosT *pos, + size_t k, + int phase, + size_t i, + double v_src, + bool is_follower) { + // 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 (r_prime == my_rank) { + (is_follower ? res.follower_self : res.leader_self).push(pos, k, phase); + } + else { + QueryCodec::push_positions(is_follower ? fq[r_prime] : lq[r_prime], pos, k, phase); + } + (is_follower ? fs[r_prime] : ls[r_prime]).push_back(i); + if (capture_values) { + (is_follower ? fv[r_prime] : lv[r_prime]).push_back(v_src); + } + }; + // 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) { if (!rotation_dynamic_gate(only_rotate_len_k, mono_pop, cut_st, abs_c)) { return; } - Monomial new_mono; - size_t overlap = 0; - int phase_factor = 0; - emit_term_products(*op.store, i, ectx, new_mono, overlap, phase_factor); + const auto p = emit_term_products(ham, i, ectx, gen_pos.data(), gen_pop, pbuf.data()); + assert(p.k == mono_pop + gen_pop - 2 * p.overlap && "the merge disagrees with the popcount identity"); // Structural cutoff on the partner M⊕G, unless upper_atol rescues it (CutoffContext::is_above_upper). - const size_t new_pop = mono_pop + gen_pop - 2 * overlap; - 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_digest(p.k, p.d); + const bool struct_pass = keep.value_or(false) || (!keep.has_value() && cutoff_eval(p.new_mono)); if (!struct_pass && !cut_st.is_above_upper(abs_c)) { return; } - 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); - } - } + const int phase = A::emit_phase(p.phase_factor, mono_pop, gen_pop, p.overlap); + push(p.new_mono, pbuf.data(), p.k, phase, i, v_src, is_follower); }; // Pass 1 and pass 2 stay fused over `nz`: splitting them regressed measurably, as `nz` spills L1 @@ -327,9 +403,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, off the measured mean of 5.33 positions; wider terms grow the buffer. + const size_t pq = QueryCodec::kReservePositionsPerQuery; + res.leader_self.reserve(n_anti - n_foll, pq); ls[my_rank].reserve(n_anti - n_foll); - fq[my_rank].reserve(n_foll * kQueryWords); + res.follower_self.reserve(n_foll, pq); fs[my_rank].reserve(n_foll); } auto derive_coeff = [&](size_t i) -> std::pair { diff --git a/cpp/monoprop/detail/evolution/layer_build/SparseQuery.h b/cpp/monoprop/detail/evolution/layer_build/SparseQuery.h new file mode 100644 index 00000000..dab6e4e8 --- /dev/null +++ b/cpp/monoprop/detail/evolution/layer_build/SparseQuery.h @@ -0,0 +1,325 @@ +// 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, gap-coded. Record order is preserved everywhere, because it is the floating-point +// accumulation order (Resolve.h mints misses in it). +// +// Header in the low bits of word 0, then the payload LSB-first, both by explicit shift, never punning: +// [0..1] phase+1 (emit_phase is TERNARY) [2..6] k, 31 escaping to a following kLongKBits-wide k +// then [kGwBits] gw then pos[0] raw at kPosBits, then k-1 gaps of gw bits. +// +// ONE form, no mode field and no per-record argmin. Gap coding is never wider than raw lanes, because +// gw = bit_width(max gap) <= bit_width(kBits - 1) = kPosBits, hence kPosBits + (k-1)*gw <= k*kPosBits. A +// raw kBits mask is narrower only once k*kPosBits > kBits, i.e. from k = 34 at 128 modes, where this +// record costs one word more: measured on 0 of 106,368 captured records (max k = 23, at pauli cutoff +// 12), and that is the price of one code path. Dropping the argmin also dropped the stack array it +// sized -- gap coding emits bits monotonically, so the encoder streams straight into `buf` and no longer +// zeroes 48 B per push. +template +struct SparseQuery { + using PosT = uint16_t; + + 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 size_t kPhaseBits = 2; + static constexpr size_t kKBits = 5; + //: k is a popcount of a kBits bitset, so the escape can never need more than this. + static constexpr size_t kLongKBits = static_cast(std::bit_width(kBits)); + static constexpr size_t kGwBits = 4; + static constexpr size_t kKEscape = (1U << kKBits) - 1U; + static constexpr size_t kHeaderBits = kPhaseBits + kKBits + kGwBits; + static_assert(kPosBits <= (1U << kGwBits) - 1U, "gw <= kPosBits must fit the header's gap-width field"); + static_assert(kHeaderBits + kLongKBits <= 64, "the widest header must be readable from word 0 alone"); + + //: A popcount cannot exceed the width, which the old 65535 never said. + static constexpr size_t kMaxPositions = kBits; + + // ---- bit stream ------------------------------------------------------------------------------- + + // Streams into `buf`: one accumulator, flushed when a word fills, in place of an array sized by the + // worst case of three encodings. + struct Writer { + VecZ &buf; + uint64_t cur = 0; + size_t nbits = 0; // bits held in cur, always < 64 + size_t words = 0; + + [[gnu::always_inline]] auto put(uint64_t v, size_t width) noexcept -> void { + if (width == 0) { + assert(v == 0 && "a zero-width field cannot carry a value"); + return; + } + // 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; + } + cur |= v << nbits; + if (nbits + width < 64) { + nbits += width; + return; + } + buf.push_back(static_cast(cur)); + ++words; + // nbits == 0 only at width == 64, where every bit is already in cur; `v >> 64` would be UB. + cur = (nbits == 0) ? 0 : (v >> (64U - nbits)); + nbits = nbits + width - 64U; + } + + auto flush() noexcept -> void { + if (nbits != 0) { + buf.push_back(static_cast(cur)); + ++words; + cur = 0; + nbits = 0; + } + } + }; + + 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 { + 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.phase = static_cast(w0 & 0x3U) - 1; + h.k = static_cast((w0 >> kPhaseBits) & kKEscape); + h.bits = kPhaseBits + kKBits; + if (h.k == kKEscape) { + h.k = static_cast((w0 >> h.bits) & ((uint64_t{1} << kLongKBits) - 1U)); + h.bits += kLongKBits; + } + h.gw = static_cast((w0 >> h.bits) & ((uint64_t{1} << kGwBits) - 1U)); + h.bits += kGwBits; + return h; + } + + [[nodiscard]] static constexpr auto header_bits_for(size_t k) noexcept -> size_t { + return kHeaderBits + ((k >= kKEscape) ? kLongKBits : 0U); + } + + [[nodiscard]] static constexpr auto gap_bits(size_t k, size_t gw) noexcept -> size_t { + return header_bits_for(k) + ((k == 0) ? 0U : kPosBits + (k - 1U) * gw); + } + [[nodiscard]] static constexpr auto words_of(size_t bits) noexcept -> size_t { return (bits + 63U) / 64U; } + + //: The record's word count, from the header alone: k does not determine it, gw does too. + [[nodiscard]] static constexpr auto words_of_header(const Header &h) noexcept -> size_t { + return words_of(gap_bits(h.k, h.gw)); + } + + [[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 width admits"); + assert(phase >= -1 && phase <= 1 && "emit_phase is ternary: rotation_sign, or REAL_PARTS entry"); + for (size_t j = 1; j < k; ++j) { + assert(pos[j] > pos[j - 1] && "positions must be strictly ascending"); + } + + const size_t gw = gap_width(pos, k); + Writer w{buf}; + 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); + } + 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); + } + } + w.flush(); + assert(w.words == words_of(gap_bits(k, gw)) && "encoder wrote a different width than it costed"); + return w.words; + } + + // ---- 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.k != 0) { + auto prev = static_cast(r.get(kPosBits)); + out[0] = static_cast(prev); + for (size_t j = 1; j < h.k; ++j) { + prev += static_cast(r.get(h.gw)) + 1U; + out[j] = static_cast(prev); + } + } + const size_t next = off + words_of_header(h); + 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.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; + } + // gw is the MAXIMUM gap width: too small truncates a gap silently, too large wastes bits. + size_t g = 0; + for (size_t j = 1; j < h.k; ++j) { + const auto b = static_cast(std::bit_width(static_cast(pos[j] - pos[j - 1] - 1))); + g = (b > g) ? b : g; + } + return g == h.gw; + } + + // 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, for callers that hold only a bitset; the emit path merges positions. + static auto push_mono(VecZ &buf, const Monomial &mono, int phase) -> size_t { + const size_t k = mono.count(); + if (k <= kStackPositions) { + PosT scratch[kStackPositions]; + size_t j = 0; + for (size_t b = mono.find_first(); b < mono.size(); b = mono.find_next(b)) { + scratch[j++] = static_cast(b); + } + assert(j == k && "find_first/find_next walk disagrees with count()"); + return push(buf, scratch, k, phase); + } + std::vector scratch(k); + size_t j = 0; + for (size_t b = mono.find_first(); b < mono.size(); b = mono.find_next(b)) { + scratch[j++] = static_cast(b); + } + assert(j == k && "find_first/find_next walk disagrees with count()"); + return push(buf, scratch.data(), k, phase); + } +}; + +} // namespace monoprop::detail diff --git a/cpp/monoprop/detail/operator/OperatorIndex.h b/cpp/monoprop/detail/operator/OperatorIndex.h index 8a3c29b2..be728232 100644 --- a/cpp/monoprop/detail/operator/OperatorIndex.h +++ b/cpp/monoprop/detail/operator/OperatorIndex.h @@ -17,8 +17,10 @@ #include #include #include +#include #include #include +#include #include #include #include @@ -96,6 +98,9 @@ class OperatorIndex { [[nodiscard]] auto size() const -> size_t { return size_; } + // Rows that exceeded inline_width_ and spilled; observable so a test can compare the two insert paths. + [[nodiscard]] auto overflow_size() const -> size_t { return overflow_.size(); } + auto reserve(size_t n) -> void { reserve_rows(n); reserve_index(n); @@ -137,6 +142,33 @@ class OperatorIndex { } } + // set() from the row's own form: a row IS an ascending position list. Same postcondition as set(), + // including the dropped stale overflow entry. + // + // Precondition: `pos` strictly ascending, every entry < 2*NumModes. A violation is silent in release + // -- an unsorted row simply never matches, and an out-of-range one decodes to a different term. + auto set_positions(size_t i, const PosT *pos, size_t count) -> void { + assert(std::adjacent_find(pos, pos + count, std::greater_equal{}) == pos + count + && "row positions must be strictly ascending"); + assert((count == 0 || static_cast(pos[count - 1]) < 2 * NumModes) && "row position out of range"); + PosT *row = &rows_[i * stride_]; + if (count > inline_width_) { + // The spill path has no position array, so build the dense form -- only here. + row[0] = kOverflowMarker; + value_type mono; + for (size_t j = 0; j < count; ++j) { + mono.set(pos[j]); + } + overflow_[i] = mono; + return; + } + if (!overflow_.empty()) { + overflow_.erase(i); + } + row[0] = static_cast(count); + std::copy_n(pos, count, row + 1); + } + [[nodiscard]] auto row(size_t i) const -> value_type { const PosT c = rows_[i * stride_]; if (c == kOverflowMarker) { @@ -170,6 +202,19 @@ class OperatorIndex { } return overflow_.at(i).count(); } + // The row's stored ascending positions; (nullptr, 0) for a spilled row, and invalidated by any insert. + struct RowPositions { + const PosT *pos; + size_t count; + [[nodiscard]] auto inlined() const -> bool { return pos != nullptr; } + }; + [[nodiscard]] auto row_positions(size_t i) const -> RowPositions { + const PosT c = rows_[i * stride_]; + if (c == kOverflowMarker) { + return {nullptr, 0}; + } + return {&rows_[(i * stride_) + 1], static_cast(c)}; + } [[nodiscard]] auto memory_bytes() const -> size_t { size_t total = rows_.capacity() * sizeof(PosT); total += overflow_.size() * (sizeof(value_type) + sizeof(size_t) + 24); @@ -233,6 +278,66 @@ class OperatorIndex { } } + // find_batch over ascending position lists: query q is pos_flat[pos_off[q] .. pos_off[q] + k_of[q]). + // Identical results to find_batch on the monomials those positions describe. Same three-stage + // prefetch pipeline, so the positions stay the currency without giving up find_batch's shape. + auto find_batch_positions(const PosT *pos_flat, + const size_t *pos_off, + const uint32_t *k_of, + size_t n, + size_t *out, + uint32_t *hash_out = nullptr) const -> void { + static constexpr size_t G = 16; + std::array hh; + std::array sp; + std::array cand; + for (size_t base = 0; base < n; base += G) { + const size_t g = std::min(G, n - base); + for (size_t j = 0; j < g; ++j) { + hh[j] = fold_hash_positions(pos_flat + pos_off[base + j], k_of[base + j]); + sp[j] = spread(hh[j]); + __builtin_prefetch(&table_.slots[sp[j] & table_.mask], 0, 0); + } + if (hash_out != nullptr) { + std::copy_n(hh.begin(), g, hash_out + base); + } + for (size_t j = 0; j < g; ++j) { + cand[j] = kEmptySlot; + if (table_.count == 0) { + continue; + } + cand[j] = probe_hash_match_(hh[j], sp[j] & table_.mask); + if (cand[j] != kEmptySlot) { + __builtin_prefetch(&rows_[static_cast(cand[j]) * stride_], 0, 0); + } + } + for (size_t j = 0; j < g; ++j) { + const size_t q = base + j; + const PosT *qpos = pos_flat + pos_off[q]; + const size_t qk = k_of[q]; + if (cand[j] == kEmptySlot) { + out[q] = kNotFound; + } + else if (row_eq_positions(static_cast(cand[j]), qpos, qk)) { + out[q] = static_cast(cand[j]); + } + else { + // A 32-bit collision: rare enough to walk the chain from the top rather than resume it. + out[q] = find_positions_(hh[j], qpos, qk); + } + } + } + } + + // fold_hash of the monomial `pos` describes, through the same fold, so it is equal by construction. + [[nodiscard]] static auto fold_hash_positions(const PosT *pos, size_t count) noexcept -> uint32_t { + key_type mono; + for (size_t j = 0; j < count; ++j) { + mono.set(pos[j]); + } + return fold_hash(mono); + } + // Insert-or-no-op. Row at `value` must already be written (the confirm reads dense rows). auto emplace(const key_type &key, mapped_type value) -> void { check_index_fits(value); @@ -251,12 +356,35 @@ class OperatorIndex { // Insert n distinct rows with consecutive indices [base, base+n). Rows must already be written. template auto bulk_insert(size_t n, mapped_type base, KeyFn &&key_at) -> void { + if (n == 0) { + return; + } + // Delegating means both entry points share one insert loop, prefetch pipeline included. + bulk_insert_hashed(n, base, [&](size_t k) { return fold_hash(key_at(k)); }); + } + // bulk_insert with the hashes already in hand: same precondition (n distinct rows, already written, + // at consecutive indices) and the same slot assignment. `hashes[k]` MUST be fold_hash of the key of + // row base+k -- a wrong one leaves the row unfindable, which surfaces later as a duplicate insert. + // + // Group-prefetched like find_batch: correctness does not depend on it (a prefetch is a hint and the + // insert re-reads the slot), but hash_at is called exactly ONCE per element and buffered. + template + auto bulk_insert_hashed(size_t n, mapped_type base, HashFn &&hash_at) -> void { if (n == 0) { return; } check_index_fits(base + n - 1); - for (size_t k = 0; k < n; ++k) { - insert_slot_(static_cast(base + k), fold_hash(key_at(k))); + static constexpr size_t G = 16; // same group width as find_batch, for the same reason + std::array hh; + for (size_t b = 0; b < n; b += G) { + const size_t g = std::min(G, n - b); + for (size_t j = 0; j < g; ++j) { + hh[j] = hash_at(b + j); + __builtin_prefetch(&table_.slots[spread(hh[j]) & table_.mask], /*rw=*/1, /*locality=*/0); + } + for (size_t j = 0; j < g; ++j) { + insert_slot_(static_cast(base + b + j), hh[j]); + } } } template @@ -386,6 +514,41 @@ 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; + } + // std::equal, i.e. a memcmp CALL, and MEASURED to be the right choice: replacing it with the + // obvious scalar loop cost 54.6M instructions on the pauli cell, because glibc's AVX2 memcmp + // beats a byte loop even at the ~5 PosT this compares. Do not "optimise" the call away again. + return std::equal(q, q + qk, &rows_[(i * stride_) + 1]); + } + + // find()'s chain walk for a position-list key, hash already folded; only the collision arm reaches it. + [[nodiscard]] auto find_positions_(uint32_t h, const PosT *q, size_t qk) const -> size_t { + if (table_.count == 0) { + return kNotFound; + } + for (size_t s = spread(h) & table_.mask;; s = (s + 1) & table_.mask) { + const Slot &e = table_.slots[s]; + if (e.idx == kEmptySlot) { + return kNotFound; + } + if (e.h == h && row_eq_positions(static_cast(e.idx), q, qk)) { + return static_cast(e.idx); + } + } + } + static auto check_index_fits(size_t value) -> void { if (value >= kIndexCeiling) { throw TermIndexCeilingReached("OperatorIndex: operator index reached the TermIndex ceiling; rebuild with " diff --git a/cpp/tests/README.md b/cpp/tests/README.md index 6fc34490..6cd35461 100644 --- a/cpp/tests/README.md +++ b/cpp/tests/README.md @@ -73,6 +73,9 @@ name and cannot address suite-nested cases, tests use flat - **`ExchangeLayoutOracle.h`**: `build_layer_exchange_layout` — the independent reference for a layer's exchange counts and displacements, which `derive_exchange_layout` in the library is checked against. +- **`dense_query_reference.h`**: the retired dense query record, frozen as the + independent oracle for `sparse_query_tests.cpp`. Test-only, and not kept in + sync with the wire format. - **`TestData.{h,cpp}`**: the `CaseData` struct and msgpack fixture loader. - **`boost-test.cmake` / `boostAddTests.cmake`**: CMake test discovery. @@ -83,17 +86,26 @@ name and cannot address suite-nested cases, tests use flat vs a std::bitset oracle), `mpfunctions.cpp` (MP utilities + bit-flip helpers), `pauli_algebra_tests.cpp`, `majorana_cutoff_tests.cpp` (length/support cutoff, CutoffEvaluator, interleave phase, coeff encode/decode), `validation_tests.cpp` - (parameter validators), `mpi_utils_tests.cpp` (find_rank + word serialization), - `evolution_detail_tests.cpp` (MatchedEpochSet + CutoffContext), - `row_accessor_tests.cpp` (dense vs OperatorIndex row accessors). + (parameter validators), `mpi_utils_tests.cpp` (find_rank, word serialization, + scan routing agreement), `evolution_detail_tests.cpp` (MatchedEpochSet + + CutoffContext), + `row_accessor_tests.cpp` (dense vs OperatorIndex row accessors), + `sparse_monomial_tests.cpp` (the `(k, d)` cutoff predicates vs their bitset + forms). - **Operator store**: `operator_index_tests.cpp`, `inverted_index_tests.cpp`, `mp_operator_tests.cpp` (MPOperator get_state Pauli/Majorana scoring, get_operator init-map drain, update_initial_operator picture branches, - insert_absent_terms, inverted-index sync, memory estimate, deep copy). + insert_absent_terms, inverted-index sync, memory estimate, deep copy), + `bulk_insert_tests.cpp` (the grouped-prefetch insert vs a one-key-at-a-time + reference: table state and enumeration order). - **Layer build / evolution**: `build_graph_tests.cpp`, `pauli_build_layer_tests.cpp`, `fused_cos_sweep_tests.cpp`, - `fused_query_codec_tests.cpp`, `combined_recompute_equivalence.cpp` - (recompute equivalence + snapshot invariance), `exact_upper_atol_rescue.cpp`, + `sparse_query_tests.cpp` (the SparseQuery wire record against the frozen dense + oracle, plus the fused value channel), `sparse_resolve_tests.cpp` (probe and + insert from wire positions vs the dense Monomial-keyed path), + `digest_cutoff_tests.cpp` (paired_mode_count and the digest cutoff predicate + vs cutoff_sums), `combined_recompute_equivalence.cpp` (recompute equivalence + + snapshot invariance), `exact_upper_atol_rescue.cpp`, `large_cosine_storage_tests.cpp`, `gate_boundaries.cpp`. - **Graph encoding / packing**: `graph_encoding_tests.cpp` (CosineWordBuilder coalescer, checked_* overflow guards, packed-phase storage + int8 read, diff --git a/cpp/tests/bulk_insert_tests.cpp b/cpp/tests/bulk_insert_tests.cpp new file mode 100644 index 00000000..2221b340 --- /dev/null +++ b/cpp/tests/bulk_insert_tests.cpp @@ -0,0 +1,188 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// bulk_insert prefetches 16 slot addresses; being a pure hint it must leave the table in EXACTLY the +// state an unpipelined loop leaves it in, slot order included -- for_each makes it Python-visible. + +#include + +#include +#include +#include +#include +#include +#include + +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/operator/OperatorIndex.h" + +using namespace monoprop; + +namespace { + +constexpr size_t kN = 250; +using Index = detail::OperatorIndex; + +auto draw_distinct(std::mt19937_64 &rng, size_t n) -> std::vector> { + std::vector> out; + std::set> seen; + std::uniform_int_distribution bit(0, Monomial::size() - 1); + std::uniform_int_distribution pop(0, 12); + while (out.size() < n) { + Monomial m; + const size_t k = pop(rng); + for (size_t placed = 0; placed < k;) { + const size_t b = bit(rng); + if (!m.test(b)) { + m.set(b); + ++placed; + } + } + std::vector key; + for (size_t w = 0; w < Monomial::num_words(); ++w) { + key.push_back(m.word(w)); + } + if (seen.insert(key).second) { + out.push_back(m); + } + } + return out; +} + +// Deliberately NOT reserved: rehash_if_needed firing mid-group frees the table already-issued +// addresses point into, the interesting case for a prefetch. +auto build(const std::vector> &terms) -> std::unique_ptr { + auto idx = std::make_unique(); + const size_t base = idx->grow_rows_geometric(terms.size()); + for (size_t k = 0; k < terms.size(); ++k) { + idx->set(base + k, terms[k]); + } + idx->bulk_insert(terms.size(), base, [&](size_t k) -> const Monomial & { return terms[k]; }); + return idx; +} + +// The oracle: N one-key calls, so every group is of size one and none of the grouped loop's boundary +// arithmetic runs. It shares insert_slot_ but not the GROUPING, which is the thing under test. +auto build_reference(const std::vector> &terms) -> std::unique_ptr { + auto idx = std::make_unique(); + const size_t base = idx->grow_rows_geometric(terms.size()); + for (size_t k = 0; k < terms.size(); ++k) { + idx->set(base + k, terms[k]); + idx->bulk_insert(1, base + k, [&](size_t) -> const Monomial & { return terms[k]; }); + } + return idx; +} + +} // namespace + +BOOST_AUTO_TEST_CASE(bulk_insert_finds_every_key) { + std::mt19937_64 rng(20260814); + const auto terms = draw_distinct(rng, 4000); + const auto idx = build(terms); + BOOST_REQUIRE_EQUAL(idx->size(), terms.size()); + for (size_t i = 0; i < terms.size(); ++i) { + const auto found = idx->find(terms[i]); + BOOST_REQUIRE(found.has_value()); + BOOST_TEST(*found == i); + } +} + +BOOST_AUTO_TEST_CASE(bulk_insert_batch_find_agrees_with_one_key_at_a_time) { + std::mt19937_64 rng(20260815); + const auto terms = draw_distinct(rng, 4000); + const auto absent = draw_distinct(rng, 500); + const auto grouped = build(terms); + const auto one_by_one = build_reference(terms); + + std::vector a(terms.size(), 0); + std::vector b(terms.size(), 0); + one_by_one->find_batch(terms.data(), terms.size(), a.data()); + grouped->find_batch(terms.data(), terms.size(), b.data()); + for (size_t i = 0; i < terms.size(); ++i) { + BOOST_TEST(a[i] == b[i]); + } + + // A table that answered everything would satisfy the loop above and prove nothing. + std::vector ma(absent.size(), 0); + std::vector mb(absent.size(), 0); + one_by_one->find_batch(absent.data(), absent.size(), ma.data()); + grouped->find_batch(absent.data(), absent.size(), mb.data()); + size_t genuinely_absent = 0; + for (size_t i = 0; i < absent.size(); ++i) { + BOOST_TEST(ma[i] == mb[i]); + if (!one_by_one->find(absent[i]).has_value()) { + ++genuinely_absent; + } + } + BOOST_TEST(genuinely_absent > 0U); +} + +BOOST_AUTO_TEST_CASE(bulk_insert_preserves_the_enumeration_order) { + // Identical, not merely equal as a set: this sequence is what the Python API returns terms in. + std::mt19937_64 rng(20260816); + const auto terms = draw_distinct(rng, 4000); + const auto grouped = build(terms); + const auto one_by_one = build_reference(terms); + + std::vector order_ref; + std::vector order_grouped; + one_by_one->for_each([&](const Monomial &, size_t i) { order_ref.push_back(i); }); + grouped->for_each([&](const Monomial &, size_t i) { order_grouped.push_back(i); }); + + BOOST_REQUIRE_EQUAL(order_grouped.size(), order_ref.size()); + BOOST_REQUIRE_EQUAL(order_ref.size(), terms.size()); + // Both coming out in index order would make the comparison blind to a reshuffle preserving it. + bool is_sorted_by_index = true; + for (size_t i = 1; i < order_ref.size(); ++i) { + if (order_ref[i] < order_ref[i - 1]) { + is_sorted_by_index = false; + break; + } + } + BOOST_TEST(!is_sorted_by_index); + for (size_t i = 0; i < order_ref.size(); ++i) { + BOOST_TEST(order_grouped[i] == order_ref[i]); + } +} + +BOOST_AUTO_TEST_CASE(bulk_insert_handles_a_partial_final_group) { + // The group width is 16, so these sizes straddle the boundary, including fewer than one group. + std::mt19937_64 rng(20260817); + for (const size_t n : {size_t{1}, size_t{15}, size_t{16}, size_t{17}, size_t{31}, size_t{33}}) { + const auto terms = draw_distinct(rng, n); + const auto grouped = build(terms); + const auto one_by_one = build_reference(terms); + BOOST_REQUIRE_EQUAL(grouped->size(), n); + for (size_t i = 0; i < n; ++i) { + const auto fa = one_by_one->find(terms[i]); + const auto fb = grouped->find(terms[i]); + BOOST_REQUIRE(fa.has_value()); + BOOST_REQUIRE(fb.has_value()); + BOOST_TEST(*fb == *fa); + } + } +} + +BOOST_AUTO_TEST_CASE(bulk_insert_of_nothing_is_a_no_op) { + std::mt19937_64 rng(20260818); + const auto terms = draw_distinct(rng, 100); + auto idx = build(terms); + const size_t before = idx->size(); + const auto key_at = [&](size_t k) -> const Monomial & { return terms[k]; }; + idx->bulk_insert(0, 0, key_at); + BOOST_TEST(idx->size() == before); + for (size_t i = 0; i < terms.size(); ++i) { + BOOST_REQUIRE(idx->find(terms[i]).has_value()); + } +} diff --git a/cpp/tests/dense_query_reference.h b/cpp/tests/dense_query_reference.h new file mode 100644 index 00000000..f190ce26 --- /dev/null +++ b/cpp/tests/dense_query_reference.h @@ -0,0 +1,72 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// The retired dense query record, frozen here as a deliberately independent oracle for +// sparse_query_tests.cpp's differential: test-only, and it does NOT follow the wire format. + +#pragma once + +#include +#include + +#include "monoprop/TypeAliases.h" +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/evolution/layer_build/Common.h" +#include "monoprop/detail/mpi/MPIUtils.h" + +namespace monoprop::test_ref { + +using monoprop::VecZ; + +// W monomial words + one ±1 phase word. +template +inline constexpr size_t kQueryWords = mpi_detail::kWords + 1; + +template +inline constexpr size_t kQueryWordsFused = kQueryWords + 1; + +template +inline auto query_push(VecZ &buf, const Monomial &mono, int phase) -> void { + mpi_detail::append_monomial_words(mono, buf); + buf.push_back(static_cast(static_cast(phase))); +} + +template > +inline auto query_read(const VecZ &buf, size_t q, Monomial &mono_out, int &phase_out) -> void { + const size_t base = q * QW; + mono_out = mpi_detail::read_monomial_from_words(buf, base); + phase_out = static_cast(static_cast(buf[base + mpi_detail::kWords])); +} + +template +inline auto query_value(const VecZ &buf, size_t q) -> double { + return detail::decode_value(buf[q * kQueryWordsFused + mpi_detail::kWords + 1]); +} + +// Requires v.size() == q.size()/kQueryWords: exactly one value per query record. +template +inline auto build_fused_query_value(const VecZ &q, const std::vector &v, VecZ &out) -> void { + constexpr size_t W = kQueryWords; + const size_t nq = q.empty() ? 0 : q.size() / W; + out.clear(); + out.reserve(nq * kQueryWordsFused); + for (size_t i = 0; i < nq; ++i) { + out.insert(out.end(), + q.begin() + static_cast(i * W), + q.begin() + static_cast((i + 1) * W)); + out.push_back(detail::encode_value(v[i])); + } +} + +} // namespace monoprop::test_ref diff --git a/cpp/tests/digest_cutoff_tests.cpp b/cpp/tests/digest_cutoff_tests.cpp new file mode 100644 index 00000000..b34a6114 --- /dev/null +++ b/cpp/tests/digest_cutoff_tests.cpp @@ -0,0 +1,143 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// paired_mode_count's d, differentially against the bitset cutoff_sums: the IDENTITY +// (d == popcount_sum - or_sum) and the PREDICATE built on it are separable, so they are separate cases. + +#include + +#include +#include +#include +#include +#include + +#include "monoprop/algebra/Algebra.h" +#include "monoprop/algebra/AlgebraCommon.h" +#include "monoprop/detail/evolution/CutoffContext.h" +#include "monoprop/detail/evolution/layer_build/Scan.h" +#include "monoprop/detail/operator/MPOperator.h" + +using namespace monoprop; + +namespace { + +// indices_to_bitset places every bit at or above the active offset: the precondition inherited here. +template +auto draw_well_formed(std::mt19937_64 &rng, size_t logical, size_t weight) -> Monomial { + VecZ idx; + std::uniform_int_distribution dist(0, (2 * logical) - 1); + while (idx.size() < weight) { + const size_t v = dist(rng); + bool dup = false; + for (const auto x : idx) { + dup = dup || (x == v); + } + if (!dup) { + idx.push_back(v); + } + } + return indices_to_bitset(idx); +} + +// Weight spans both extremes: d == 0 and k == 2d are the cases the cutoffs branch on. +template +auto check_identity(std::mt19937_64 &rng, size_t logical, size_t &checked, size_t &paired_seen) -> void { + for (size_t weight = 1; weight <= 2 * logical && weight <= 24; ++weight) { + for (int rep = 0; rep < 40; ++rep) { + const auto mono = draw_well_formed(rng, logical, weight); + const auto sums = cutoff_sums(mono, logical); + const size_t d = paired_mode_count(mono); + BOOST_REQUIRE_EQUAL(d, sums.popcount_sum - sums.or_sum); + const auto rebuilt = cutoff_sums(sums.popcount_sum, d); + BOOST_REQUIRE_EQUAL(rebuilt.xor_sum, sums.xor_sum); + BOOST_REQUIRE_EQUAL(rebuilt.or_sum, sums.or_sum); + BOOST_REQUIRE_EQUAL(rebuilt.popcount_sum, sums.popcount_sum); + paired_seen += static_cast(sums.xor_sum == 0); + ++checked; + } + } +} + +} // namespace + +BOOST_AUTO_TEST_CASE(paired_mode_count_matches_cutoff_sums_across_widths) { + std::mt19937_64 rng(0xD16E57U); + size_t checked = 0; + size_t paired_seen = 0; + + check_identity<32>(rng, 32, checked, paired_seen); // W = 64, one word, no active offset + check_identity<32>(rng, 30, checked, paired_seen); // W = 64, active_bit_offset = 4 + check_identity<48>(rng, 45, checked, paired_seen); // W = 96 -- not a multiple of 64 + check_identity<64>(rng, 64, checked, paired_seen); // W = 128, exactly two words + check_identity<128>(rng, 120, checked, paired_seen); + check_identity<256>(rng, 250, checked, paired_seen); // the production shape + + BOOST_TEST(checked > 3000U); + // The identity is only interesting on the fully-paired branch, so the draw must reach it. + BOOST_TEST(paired_seen > 0U); +} + +BOOST_AUTO_TEST_CASE(paired_mode_count_exhaustive_at_small_width) { + constexpr size_t kN = 5; // W = 10 + size_t paired = 0; + for (uint64_t bits = 0; bits < (uint64_t{1} << (2 * kN)); ++bits) { + Monomial mono; + for (size_t b = 0; b < 2 * kN; ++b) { + if ((bits >> b) & 1U) { + mono.set(b); + } + } + const auto sums = cutoff_sums(mono, kN); + BOOST_REQUIRE_EQUAL(paired_mode_count(mono), sums.popcount_sum - sums.or_sum); + paired += static_cast(sums.xor_sum == 0); + } + BOOST_TEST(paired == 32U); // 2^5: each mode independently empty or doubly occupied +} + +// At the PREDICATE level, not the scan level: cutoff_sums is the independent form to compare against. +BOOST_AUTO_TEST_CASE(digest_predicate_matches_cutoff_sums_predicate) { + std::mt19937_64 rng(0xC0FFEEU); + size_t checked = 0; + size_t kept = 0; + size_t rejected = 0; + + // The popcount <= cutoff early-out is the asymmetry between the two, so cutoffs straddle it. + for (const unsigned int cutoff : {1U, 2U, 4U, 6U, 10U, 20U}) { + for (const bool support : {false, true}) { + constexpr size_t kN = 32; + constexpr size_t kLogical = 30; + const CutoffFn fn = support ? CutoffFn{detail::SupportCutoff{cutoff, kLogical}} + : CutoffFn{detail::LengthCutoff{cutoff, kLogical}}; + const detail::CutoffEvaluator eval(fn); + for (size_t w = 1; w <= 12; ++w) { + for (int rep = 0; rep < 40; ++rep) { + const auto mono = draw_well_formed(rng, kLogical, w); + const size_t k = mono.count(); + const auto digest = eval.passes_from_dense(mono, k); + BOOST_REQUIRE(digest.has_value()); // a concrete cutoff must always decide + const bool reference = eval.passes_with_popcount(mono, k); + BOOST_REQUIRE_EQUAL(*digest, reference); + ++checked; + kept += static_cast(*digest); + rejected += static_cast(!*digest); + } + } + } + } + BOOST_TEST(checked > 5000U); + // A sweep that only ever kept would agree with any predicate that returns true. + BOOST_TEST(kept > 0U); + BOOST_TEST(rejected > 0U); +} diff --git a/cpp/tests/evolution_detail_tests.cpp b/cpp/tests/evolution_detail_tests.cpp index 91f2662b..bade37f4 100644 --- a/cpp/tests/evolution_detail_tests.cpp +++ b/cpp/tests/evolution_detail_tests.cpp @@ -154,8 +154,17 @@ BOOST_AUTO_TEST_CASE(self_resolve_mark_bounded_by_combined_size) { matched, combined_size, RecordingSink{}); - detail::query_push<8>(eng.queries_r[0], terms[1], 1); - detail::query_push<8>(eng.queries_r[0], terms[5], -1); + // The self leg is staged as positions, never encoded, so this feeds the stage the scan would fill. + using Eng = detail::LayerBuildEngine<8, RecordingSink>; + const auto stage_self = [&eng](const Monomial<8> &m, int phase) { + std::vector pos; + for (size_t b = m.find_first(); b < m.size(); b = m.find_next(b)) { + pos.push_back(static_cast(b)); + } + eng.self_stage_.push(pos.data(), pos.size(), phase); + }; + stage_self(terms[1], 1); + stage_self(terms[5], -1); eng.src_idx_r[0] = {0, 2}; eng.resolve_self_queries(/*is_leader_pass=*/true); diff --git a/cpp/tests/fused_query_codec_tests.cpp b/cpp/tests/fused_query_codec_tests.cpp deleted file mode 100644 index 72f76763..00000000 --- a/cpp/tests/fused_query_codec_tests.cpp +++ /dev/null @@ -1,125 +0,0 @@ -// Copyright 2026 Algorithmiq -// -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at -// -// http://www.apache.org/licenses/LICENSE-2.0 -// -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. - -#include - -#include -#include -#include - -#include "monoprop/TypeAliases.h" -#include "monoprop/core/Monomial.h" -#include "monoprop/detail/evolution/layer_build/Common.h" - -// Query+value fusion codec: the fused R>1 exchange rides the source coefficient (v_src) on each query -// record as a trailing bit-cast word so one alltoallv carries both streams. Round-tripping must be -// byte-for-byte, including the FP corner cases a lossy value channel would mangle. - -namespace { - -using namespace monoprop; -using monoprop::detail::build_fused_query_value; -using monoprop::detail::kQueryWords; -using monoprop::detail::kQueryWordsFused; -using monoprop::detail::query_push; -using monoprop::detail::query_read; -using monoprop::detail::query_value; - -constexpr size_t kModes = 8; // 2*kModes = 16 majorana bits, one 64-bit word - -// A deterministic, distinct majorana bit pattern per record index. -auto make_mono(size_t r) -> Monomial { - Monomial m; - for (size_t b = 0; b < 2 * kModes; ++b) { - if (((r * 2654435761u + b * 40503u) & 3u) == 0u) { - m.set(b); - } - } - return m; -} - -BOOST_AUTO_TEST_CASE(fused_record_roundtrip_exact) { - const std::vector phases = {1, -1, 1, -1, 1, 1, -1}; - const std::vector values = { - 0.0, - -0.0, - 1.0, - -1.0, - 3.141592653589793, - -2.718281828459045e-300, // near-denormal magnitude - std::numeric_limits::min(), // smallest normal - }; - const size_t nq = values.size(); - - VecZ plain; - std::vector> monos(nq); - for (size_t r = 0; r < nq; ++r) { - monos[r] = make_mono(r); - query_push(plain, monos[r], phases[r]); - } - BOOST_REQUIRE_EQUAL(plain.size(), nq * kQueryWords); - - VecZ fused; - build_fused_query_value(plain, values, fused); - BOOST_REQUIRE_EQUAL(fused.size(), nq * kQueryWordsFused); - - for (size_t q = 0; q < nq; ++q) { - Monomial m_out; - int ph_out = 0; - query_read>(fused, q, m_out, ph_out); - BOOST_CHECK(m_out == monos[q]); - BOOST_CHECK_EQUAL(ph_out, phases[q]); - // Compare the raw payload, so -0.0 and denormals stay distinguished from 0.0. - const double v_out = query_value(fused, q); - BOOST_CHECK(std::memcmp(&v_out, &values[q], sizeof(double)) == 0); - } -} - -// Reusing `out` across calls must leak no stale words: capacity is a high-water mark, size is exact. -// This is the reuse pattern LayerBuildEngine::combined_qv_ relies on gate to gate. -BOOST_AUTO_TEST_CASE(fused_buffer_reuse_shrinks_logical_size) { - VecZ plain_big; - std::vector vbig; - for (size_t r = 0; r < 32; ++r) { - query_push(plain_big, make_mono(r), (r % 2 == 0) ? 1 : -1); - vbig.push_back(static_cast(r) * 1.5 - 7.0); - } - VecZ out; - build_fused_query_value(plain_big, vbig, out); - const size_t cap_after_big = out.capacity(); - - VecZ plain_small; - std::vector vsmall = {42.0, -42.0, 0.25}; - for (size_t r = 0; r < vsmall.size(); ++r) { - query_push(plain_small, make_mono(100 + r), 1); - } - build_fused_query_value(plain_small, vsmall, out); - BOOST_CHECK_EQUAL(out.size(), vsmall.size() * kQueryWordsFused); - BOOST_CHECK_GE(out.capacity(), cap_after_big); - for (size_t q = 0; q < vsmall.size(); ++q) { - const double v_out = query_value(out, q); - BOOST_CHECK(std::memcmp(&v_out, &vsmall[q], sizeof(double)) == 0); - } -} - -// Empty input arises for the self slot, which resolve_self_queries clears before the exchange. -BOOST_AUTO_TEST_CASE(fused_empty_input) { - VecZ empty; - std::vector no_values; - VecZ out{1, 2, 3}; // pre-dirtied; build must clear it - build_fused_query_value(empty, no_values, out); - BOOST_CHECK(out.empty()); -} - -} // namespace diff --git a/cpp/tests/mpi_utils_tests.cpp b/cpp/tests/mpi_utils_tests.cpp index 8372e88f..45fdaf9e 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,120 @@ 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()); + } +} + +// The self-owned bucket is staged as positions, not encoded, so it is invisible to the walk above -- +// without this the r == my_rank arm of the routing decision goes unchecked. +auto check_self_ownership(const detail::SelfQueryStage<32> &stage, size_t ranks, size_t my_rank, size_t &checked) + -> void { + for (size_t q = 0; q < stage.size(); ++q) { + Monomial<32> mono; + for (size_t j = 0; j < stage.k_of[q]; ++j) { + mono.set(static_cast(stage.pos_flat[stage.pos_off[q] + j])); + } + BOOST_REQUIRE_EQUAL(find_rank<32>(mono, ranks), my_rank); + ++checked; + } +} + +} // namespace + +// The scan hashes the partner it just built; find_rank hashes what the resolve side decoded off the +// 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; + size_t self_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); + // The scan routes a self-owned partner to the stage, so bucket 0 must be empty here. + BOOST_REQUIRE(res.leader_queries[0].empty()); + BOOST_REQUIRE(res.follower_queries[0].empty()); + check_bucket_ownership(res.leader_queries, ranks, checked); + check_bucket_ownership(res.follower_queries, ranks, checked); + check_self_ownership(res.leader_self, ranks, /*my_rank=*/0, self_checked); + check_self_ownership(res.follower_self, ranks, /*my_rank=*/0, self_checked); + } + // Without this the loop above passes trivially if the scan emitted nothing. The floor is on the SUM + // because that is what is invariant across the split: the encoded counter alone fell to 797 of 1161 + // when the self-owned partners moved into the stage, with nothing going unchecked. Each arm still + // carries its own floor -- a routing bug sending everything one way leaves the sum intact -- and the + // message prints the measured 797/364 so those can be re-grounded rather than guessed. + BOOST_TEST_MESSAGE("encoded=" << checked << " staged=" << self_checked); + BOOST_TEST(checked + self_checked > 1000U); + BOOST_TEST(checked > 500U); + BOOST_TEST(self_checked > 200U); +} diff --git a/cpp/tests/partner_merge_tests.cpp b/cpp/tests/partner_merge_tests.cpp new file mode 100644 index 00000000..c12a4cf4 --- /dev/null +++ b/cpp/tests/partner_merge_tests.cpp @@ -0,0 +1,168 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// The partner merge against the dense XOR it replaces. Random pairs alone do not reach the cases that +// decide it: the merge's whole surface is how many of G's slots the source already holds, so every +// overlap in [0, gen_pop] is drawn on purpose, and the paired case is drawn separately because a +// fully-paired monomial is 94 in 20.9M on production models. + +#include + +#include +#include +#include +#include +#include + +#include "monoprop/algebra/AlgebraCommon.h" +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/evolution/layer_build/PartnerMerge.h" + +using namespace monoprop; + +namespace { + +constexpr size_t kN = 64; +constexpr size_t kBits = 2 * kN; + +// The dense reference: positions of M^G, ascending, straight off the bitset. +auto dense_partner(const Monomial &mono, const Monomial &gen) -> std::vector { + const Monomial nm = mono ^ gen; + std::vector out; + for (size_t b = nm.find_first(); b < nm.size(); b = nm.find_next(b)) { + out.push_back(b); + } + return out; +} + +auto positions_of(const Monomial &m) -> std::vector { + std::vector out; + for (size_t b = m.find_first(); b < m.size(); b = m.find_next(b)) { + out.push_back(static_cast(b)); + } + return out; +} + +// Checks the merge against the dense form on one pair, and returns the merged count so a caller can +// assert it saw work. Every field the emit site consumes is compared, not just the positions. +auto check_pair(const Monomial &mono, const Monomial &gen) -> size_t { + const auto src = positions_of(mono); + const auto gpos = positions_of(gen); + std::vector out(kBits); + size_t overlap = 0; + size_t d = 0; + const size_t k = + detail::merge_partner_positions(src.data(), src.size(), gpos.data(), gpos.size(), out.data(), overlap, d); + + const auto expect = dense_partner(mono, gen); + BOOST_REQUIRE_EQUAL(k, expect.size()); + for (size_t j = 0; j < k; ++j) { + BOOST_REQUIRE_EQUAL(static_cast(out[j]), expect[j]); + } + BOOST_REQUIRE_EQUAL(overlap, mono.count_and(gen)); + // d must be exactly what the cutoff digest would have folded out of the dense partner. + BOOST_REQUIRE_EQUAL(d, paired_mode_count(mono ^ gen)); + // The popcount identity the emit site asserts on. + BOOST_REQUIRE_EQUAL(k, mono.count() + gen.count() - (2 * overlap)); + return k; +} + +// A generator of `gen_pop` slots, and a source holding exactly `overlap` of them plus `extra` others. +auto build_case(std::mt19937_64 &rng, size_t gen_pop, size_t overlap, size_t extra) + -> std::pair, Monomial> { + std::vector all(kBits); + for (size_t i = 0; i < kBits; ++i) { + all[i] = i; + } + std::shuffle(all.begin(), all.end(), rng); + Monomial gen; + for (size_t i = 0; i < gen_pop; ++i) { + gen.set(all[i]); + } + Monomial mono; + for (size_t i = 0; i < overlap; ++i) { + mono.set(all[i]); // a slot G also holds: it cancels + } + for (size_t i = 0; i < extra; ++i) { + mono.set(all[gen_pop + i]); // disjoint from G: it survives + } + return {mono, gen}; +} + +} // namespace + +// Every overlap between the source and the generator, which is the branch the merge exists to take. +BOOST_AUTO_TEST_CASE(partner_merge_matches_dense_at_every_overlap) { + std::mt19937_64 rng(0xA11CEU); + size_t cases = 0; + size_t nonempty = 0; + for (size_t gen_pop = 1; gen_pop <= 6; ++gen_pop) { + for (size_t overlap = 0; overlap <= gen_pop; ++overlap) { + for (const size_t extra : {size_t{0}, size_t{1}, size_t{5}, size_t{20}}) { + for (size_t rep = 0; rep < 8; ++rep) { + const auto [mono, gen] = build_case(rng, gen_pop, overlap, extra); + BOOST_REQUIRE_EQUAL(mono.count_and(gen), overlap); // the case is the one intended + nonempty += (check_pair(mono, gen) != 0) ? 1 : 0; + ++cases; + } + } + } + } + BOOST_REQUIRE_EQUAL(cases, 6U * 4U * 8U + (1U + 2U + 3U + 4U + 5U + 6U) * 4U * 8U); + // Total cancellation (overlap == gen_pop, extra == 0) is the only empty partner, so most must not be. + BOOST_TEST(nonempty > 600U); +} + +// The paired population separately: d is the field the length cutoff's escape hatch rests on, and a +// uniform draw almost never produces a fully paired monomial. +BOOST_AUTO_TEST_CASE(partner_merge_d_matches_on_paired_monomials) { + std::mt19937_64 rng(0xBEEFU); + size_t paired_seen = 0; + for (size_t rep = 0; rep < 400; ++rep) { + // Whole modes only, so both slots of each are set and the result is fully paired. + std::vector modes(kN); + for (size_t i = 0; i < kN; ++i) { + modes[i] = i; + } + std::shuffle(modes.begin(), modes.end(), rng); + Monomial mono; + const size_t n_modes = 1 + (rng() % 8); + for (size_t i = 0; i < n_modes; ++i) { + mono.set(2 * modes[i]); + mono.set((2 * modes[i]) + 1); + } + Monomial gen; + const size_t g_modes = 1 + (rng() % 3); + for (size_t i = 0; i < g_modes; ++i) { + gen.set(2 * modes[kN - 1 - i]); + gen.set((2 * modes[kN - 1 - i]) + 1); + } + check_pair(mono, gen); + paired_seen += is_paired(mono ^ gen) ? 1 : 0; + } + // The draw is meant to land on the paired branch every time; a zero here means it stopped doing so. + BOOST_REQUIRE_EQUAL(paired_seen, 400U); +} + +// An empty generator and an empty source are both reachable (a truncated gate, a fresh row). +BOOST_AUTO_TEST_CASE(partner_merge_handles_empty_inputs) { + Monomial mono; + mono.set(4); + mono.set(5); + mono.set(70); + const Monomial empty; + BOOST_REQUIRE_EQUAL(check_pair(mono, empty), 3U); + BOOST_REQUIRE_EQUAL(check_pair(empty, mono), 3U); + BOOST_REQUIRE_EQUAL(check_pair(empty, empty), 0U); +} 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..b67752e2 --- /dev/null +++ b/cpp/tests/sparse_query_tests.cpp @@ -0,0 +1,592 @@ +// 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 the widest gap widths, which is the case a regular pattern never reaches; see +// sparse_record_reaches_the_widest_gap_width for why a narrow universe needs a run plus one outlier. +auto scattered(size_t k, size_t universe, std::mt19937_64 &rng) -> std::vector { + std::vector pool(universe); + for (size_t j = 0; j < universe; ++j) { + 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_reaches_the_widest_gap_width) { + // gw == kPosBits is the record's worst case and uniform draws CANNOT reach it at kBits=24: it needs + // one gap of at least half the universe, which only a dense run plus one far outlier forces. The + // shape of the input is a selection rule, so this generator is kept even though bitmap mode is gone. + const auto count_widest = [](auto tag, size_t universe) { + using SQ = SparseQuery; + size_t used = 0; + size_t bad = 0; + const auto tally = [&](const std::vector &pos) { + if (pos.size() < 2 || pos.size() > SQ::kMaxPositions) { + return; + } + VecZ buf; + const size_t w = SQ::push(buf, pos.data(), pos.size(), 1); + const size_t gw = SQ::gap_width(pos.data(), pos.size()); + if (gw != SQ::kPosBits) { + return; + } + ++used; + std::vector back(pos.size()); + SQ::read_positions(buf, 0, back.data()); + bad += static_cast(back != pos || SQ::k_at(buf, 0) != pos.size() + || w != SQ::words_of(SQ::gap_bits(pos.size(), gw))); + }; + 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_widest(std::integral_constant{}, 24); + const auto partial = count_widest(std::integral_constant{}, 500); + const auto bucket = count_widest(std::integral_constant{}, 256); + BOOST_TEST(narrow.first > 0U); + BOOST_TEST(partial.first > 0U); + BOOST_TEST(bucket.first > 0U); + BOOST_TEST(narrow.second == 0U); + BOOST_TEST(partial.second == 0U); + BOOST_TEST(bucket.second == 0U); +} + +BOOST_AUTO_TEST_CASE(sparse_record_survives_the_five_bit_k_escape) { + // The escape is at k = 31 now, not 63: both boundaries are here so a field-width change is caught. + for (const size_t k : {size_t{30}, size_t{31}, size_t{32}, size_t{62}, size_t{63}, size_t{64}, size_t{200}}) { + const auto pos = strided(k, 0, 3, 2048); + BOOST_REQUIRE_EQUAL(pos.size(), k); + differential<1024>(pos, 1); + } +} + +BOOST_AUTO_TEST_CASE(sparse_record_bounds_the_fully_paired_term) { + // Every bit set means every gap is 0, so gw is 0 and the payload is one raw position: 23 header + // bits + 11 = one word, where raw lanes would take 514. This is what pays for deleting the argmin. + std::vector all(2048); + for (size_t j = 0; j < all.size(); ++j) { + all[j] = static_cast(j); + } + const size_t sw = differential<1024>(all, 1); + BOOST_TEST(sw == 1U); +} + +BOOST_AUTO_TEST_CASE(sparse_record_never_exceeds_its_own_raw_lanes) { + // The one width guarantee that survives deleting the argmin, and it is exhaustive rather than + // sampled: gw = bit_width(max gap) <= kPosBits, so kPosBits + (k-1)*gw <= k*kPosBits at every k. + const auto check = [](auto tag) { + using SQ = SparseQuery; + size_t cells = 0; + size_t bad = 0; + for (size_t k = 0; k <= SQ::kMaxPositions; ++k) { + for (size_t gw = 0; gw <= SQ::kPosBits; ++gw) { + const size_t lanes = SQ::words_of(SQ::header_bits_for(k) + (k * SQ::kPosBits)); + bad += static_cast(SQ::words_of(SQ::gap_bits(k, gw)) > lanes); + ++cells; + } + } + return std::pair{cells, bad}; + }; + const auto narrow = check(std::integral_constant{}); + const auto bucket = check(std::integral_constant{}); + const auto wide = check(std::integral_constant{}); + BOOST_TEST(narrow.second == 0U); + BOOST_TEST(bucket.second == 0U); + BOOST_TEST(wide.second == 0U); + // A guarded loop that asserted nothing would pass the three above; these are the cell counts. + BOOST_TEST(narrow.first == 25U * 6U); + BOOST_TEST(bucket.first == 257U * 9U); + BOOST_TEST(wide.first == 501U * 10U); +} + +BOOST_AUTO_TEST_CASE(sparse_record_documents_what_deleting_the_argmin_cost) { + // Deleting FIXED and BITMAP has a price, and this pins it to a number rather than leaving it to be + // rediscovered. The three formulas below are the DELETED encoder's, with its own 10-bit header and + // 16-bit k escape, so the comparison is against what actually shipped in #263. + using SQ = SparseQuery<128>; + const auto old_words = [](size_t k, size_t gw) { + const size_t h = 10U + ((k >= 63U) ? 16U : 0U); + return std::min({SQ::words_of(h + (k * SQ::kPosBits)), + SQ::words_of(h + 4U + (k == 0 ? 0U : SQ::kPosBits + ((k - 1U) * gw))), + SQ::words_of(h + SQ::kBits)}); + }; + // Nothing in the supported envelope loses: Pauli cutoff 16 bounds k at 32, and the widest k ever + // captured is 23 (pauli c12, 45,296 records). + size_t crossings = 0; + for (size_t k = 0; k <= 33U; ++k) { + for (size_t gw = 0; gw <= SQ::kPosBits; ++gw) { + crossings += static_cast(SQ::words_of(SQ::gap_bits(k, gw)) > old_words(k, gw)); + } + } + BOOST_TEST(crossings == 0U); + + // Above it, a raw mask wins, and 34 is where. If a field width changes, this number moves and says so. + size_t first = SQ::kMaxPositions + 1U; + for (size_t k = 0; k <= SQ::kMaxPositions && first > SQ::kMaxPositions; ++k) { + for (size_t gw = 0; gw <= SQ::kPosBits; ++gw) { + if (SQ::words_of(SQ::gap_bits(k, gw)) > old_words(k, gw)) { + first = k; + break; + } + } + } + BOOST_TEST(first == 34U); +} + +BOOST_AUTO_TEST_CASE(sparse_record_walks_a_multi_query_buffer_exactly) { + // Mixed width, 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_is_exactly_the_gap_code_it_costed) { + // What the argmin test became: there is one closed form now, so the encoder's width and the costing + // function must agree on every draw, including the uniform ones that used to select other modes. + using SQ = SparseQuery<128>; + std::mt19937_64 rng(12345); + for (size_t trial = 0; trial < 400; ++trial) { + 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()); + BOOST_TEST(w == SQ::words_of(SQ::gap_bits(pos.size(), gwid)), + "k=" << k << " wrote " << w << " words, costed " << SQ::words_of(SQ::gap_bits(pos.size(), gwid))); + BOOST_TEST(w <= SQ::words_of(SQ::header_bits_for(pos.size()) + (pos.size() * SQ::kPosBits))); + } +} + +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()); +} + +// Positions with exactly k entries whose widest gap is exactly `gw`: one gap of 2^(gw-1) -- the smallest +// value of that bit width -- then a contiguous run. Empty if the shape does not fit the universe. +namespace { +auto gap_shaped(size_t k, size_t gw, size_t universe) -> std::vector { + std::vector v; + if (k == 0) { + return v; + } + if (gw == 0 || k == 1) { + if (gw != 0 || k > universe) { + return v; + } + for (size_t j = 0; j < k; ++j) { + v.push_back(static_cast(j)); + } + return v; + } + const size_t second = (size_t{1} << (gw - 1)) + 1U; // gap value 2^(gw-1), i.e. bit_width == gw + if (second + (k - 2U) >= universe) { + return v; + } + v.push_back(0); + for (size_t j = 0; j + 1 < k; ++j) { + v.push_back(static_cast(second + j)); + } + return v; +} +} // namespace + +BOOST_AUTO_TEST_CASE(sparse_record_round_trips_every_reachable_k_and_gap_width) { + // The whole (k, gw) surface of the one remaining form, constructed rather than drawn: a random draw + // reaches neither gw = kPosBits nor the escape boundary. differential() carries ten assertions per + // cell, and `cells` is here so a shape that stops fitting cannot silently empty the loop. + using SQ = SparseQuery<128>; + size_t cells = 0; + for (size_t k = 0; k <= 40U; ++k) { + for (size_t gw = 0; gw <= SQ::kPosBits; ++gw) { + if (k < 2U && gw > 0U) { + continue; // one gap width is reachable below k = 2, so the other rows are the same cell + } + const auto pos = gap_shaped(k, gw, SQ::kBits); + if (pos.size() != k) { + continue; + } + BOOST_REQUIRE_EQUAL(SQ::gap_width(pos.data(), k), k < 2 ? 0U : gw); + const size_t w = differential<128>(pos, (k % 3U) == 0U ? 0 : ((k % 3U) == 1U ? 1 : -1)); + BOOST_TEST(w == SQ::words_of(SQ::gap_bits(k, k < 2 ? 0U : gw))); + ++cells; + } + } + BOOST_TEST(cells == 353U); + + // The width boundary: k = kBits is a fully paired term, one word because every gap is 0. + for (const size_t k : {size_t{254}, size_t{255}, size_t{256}}) { + const auto pos = gap_shaped(k, 0, SQ::kBits); + BOOST_REQUIRE_EQUAL(pos.size(), k); + BOOST_TEST(differential<128>(pos, 1) == 1U); + } +} 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); +}