refactor: unify runtime SIMD kernel dispatch

AGENTS.md

@@ -18,7 +18,7 @@ factors estimate L2 distance or inner product.
 | `include/rabitqlib/index/{ivf,hnsw,symqg}/` | Index construction, persistence, and search |
 | `include/rabitqlib/index/{query,estimator}.hpp` | Query state and distance estimation |
 | `include/rabitqlib/simd/`, `src/simd/` | Kernel declarations, implementations, and dispatch |
-| `src/index/` | Compiled HNSW search kernels |
+| `src/index/` | Compiled SymphonyQG implementation, HNSW search kernels, and IVF candidate insertion |
 | `src/utils/cpu_features.cpp` | Generic x86 feature detection |
 | `include/rabitqlib/utils/` | Rotation, allocation, buffers, I/O, and helpers |
 | `python_bindings/` | pybind11 extension and index wrappers |
@@ -110,6 +110,9 @@ Recommended:
 - Public/generic code calls centralized dispatch entry points. Keep ISA-specific translation units
   and their flags in `CMakeLists.txt` synchronized with feature predicates in
   `src/simd/dispatch.cpp` and detection in `src/utils/cpu_features.cpp`, including HNSW source groups.
+- Use the shared resolver in `src/simd/dispatch.cpp` for cached selection, including HNSW.
+  Keep calculations in backend source files; see the dispatch coverage table in
+  [CONTRIBUTING.md](CONTRIBUTING.md#dispatch-conventions-and-coverage).
 - Dispatch resolves function pointers during static initialization. Detection must use safe generic
   code; never execute a high-ISA kernel to find out whether the CPU supports it.
 - Semantic kernel changes must cover every implementation and a backend-independent reference

CMakeLists.txt

@@ -45,9 +45,15 @@ function(rabitq_enable_warnings target)
 endfunction()
 
 set(RABITQ_COMMON_SOURCES
+    src/index/ivf_search.cpp
+    src/index/qg.cpp
     src/utils/cpu_features.cpp
     src/simd/dispatch.cpp
-    src/simd/space_float.cpp
+    src/simd/space_generic.cpp
+    src/simd/fastscan_generic.cpp
+    src/simd/quantization_generic.cpp
+    src/simd/estimator_generic.cpp
+    src/simd/matrix_generic.cpp
 )
 
 set(RABITQ_AVX2_SOURCES
@@ -56,6 +62,8 @@ set(RABITQ_AVX2_SOURCES
     src/simd/space_excode_avx2.cpp
     src/simd/space_avx2.cpp
     src/simd/fastscan_avx2.cpp
+    src/simd/estimator_avx2.cpp
+    src/simd/matrix_avx2.cpp
     src/simd/warmup_avx2.cpp
     src/simd/rotator_avx2.cpp
 )
@@ -66,6 +74,8 @@ set(RABITQ_AVX512_SOURCES
     src/simd/space_excode_avx512.cpp
     src/simd/space_avx512.cpp
     src/simd/fastscan_avx512.cpp
+    src/simd/estimator_avx512.cpp
+    src/simd/matrix_avx512.cpp
     src/simd/rotator_avx512.cpp
 )
 

CONTRIBUTING.md

@@ -231,7 +231,8 @@ cmake -S . -B build-includes -G Ninja \
 
 The script checks library sources using their compilation database and checks
 headers as main files, since this clang-tidy check does not report findings in
-included headers. Private headers are checked with AVX2 and AVX-512 flags.
+included headers. Private headers are checked with AVX2 and AVX-512 flags. New
+untracked files are included, and tracked files deleted from the working tree are skipped.
 Vendored files are excluded. The script also ignores suggestions to include
 Eigen and hnswlib implementation headers behind their existing public headers;
 these vendor snapshots lack the export annotations needed by include-cleaner.
@@ -399,6 +400,51 @@ requests labeled `duplicate` or `invalid` are omitted from release notes.
 Never execute a high-ISA implementation merely to test whether that ISA is supported; detection
 must happen in generic code first.
 
+#### Dispatch conventions and coverage
+
+- Keep backend-neutral declarations in `include/rabitqlib/simd/*_dispatch.hpp` and
+  implementations in `src/simd/*_{generic,avx2,avx512}.cpp`. HNSW keeps its existing
+  `src/index/` implementations and compatibility namespaces.
+- Select a cached function pointer through `resolve_kernel` in `src/simd/dispatch.cpp`.
+  The order is AVX-512, AVX2/FMA, then the existing generic implementation or a descriptive
+  unsupported-operation exception. Public wrappers do not repeat feature checks.
+- Preserve stricter predicates: population-count kernels need AVX512_VPOPCNTDQ, and HNSW's
+  AVX-512 core variant also needs AVX2 for its warmup implementation. Do not infer support
+  from a backend name or from `__AVX*__` macros in a public header.
+- Put calculations and scratch-storage helpers outside the dispatcher. Shared implementation
+  headers use internal linkage so independently compiled backends retain their own bodies.
+  Include the FFHT implementation inside the private kernel namespace for the same reason.
+- Pass ordinary pointers, sizes, and library-owned query state across ISA boundaries. Never
+  pass Eigen matrix/packet objects between backends. The private matrix implementation
+  header includes Eigen under a namespace selected by each backend translation unit.
+  Otherwise, Eigen emits identically named out-of-line template helpers, which the linker
+  can merge across incompatible ISA builds.
+  Keep this isolation when adding matrix kernels; do not modify the vendor snapshot.
+- Preserve existing public names and namespaces, including legacy functions with `_avx` in
+  their name that now dispatch at runtime. The explicit `_avx2` and `_avx512` entry points
+  are for selected kernels and capability-guarded backend tests.
+
+The current first-party kernel audit is summarized below. Dispatching an index operation
+covers its arithmetic kernels, not every scalar loop in construction and search.
+
+| Area | Dispatch coverage / deliberate boundary |
+| --- | --- |
+| Raw float distances, norms, packed-code products | Runtime AVX2/AVX-512 selection; generic raw-float fallback |
+| Quantizer rescale search | SIMD bounded search; the certified scalar event sweep remains the fallback |
+| Integer scalar quantization, extra-code packing, transpose, sign masks | Existing runtime selection; byte layouts and rounding rules are unchanged |
+| Standard and high-accuracy FastScan | Runtime selection; generic LUT construction stays separate from selection |
+| IVF and float SymphonyQG batch correction | Complete estimator runs in the selected backend; non-float template paths remain generic |
+| FHT/Kac rotation | Complete rotation and scaling run in selected ISA translation units; imported FFHT AVX butterflies are unchanged |
+| Float matrix rotation and PiPNN construction | Matrix products, row norms, and lower-triangle pairwise distances use isolated matrix backends |
+| HNSW search and IVF centroid routing | Cached HNSW search selection; centroid routing uses the common raw-distance dispatcher |
+| Quantization orchestration, reconstruction, non-float utilities | Template/control code remains generic; no blanket native tuning or reduction-order rewrite |
+| Graph scheduling, candidate queues, I/O, allocation, random initialization | Generic control code; IVF one-bit candidate insertion stays in a small compiled function to avoid inlining-induced register spills; thread scheduling and seeds remain caller-owned |
+| Example KMeans training | Uses external FAISS, whose build and dispatch are independent of this package |
+
+Portable wheels continue to disable `RABITQ_ENABLE_NATIVE_OPTIMIZATION`. Adding an optimized
+backend does not add support for generic-CPU quantized search or AArch64; those require
+complete implementations and separate compatibility validation.
+
 ### Change quantization or packing
 
 Check all of these together:

docs/docs/index/qg.md

@@ -39,17 +39,19 @@ window; `ef` controls the query search window. Python defaults to one thread.
 
 ### C++
 
-The C++ API uses `rabitqlib::symqg::QuantizedGraph` and `QGBuilder`:
+The C++ API uses the float-only `rabitqlib::symqg::QuantizedGraph` and `QGBuilder`.
+`QuantizedGraph` is a non-template class; use `QuantizedGraph` instead of
+`QuantizedGraph<float>`. Its implementation is compiled in `src/index/qg.cpp`:
 
 ```cpp
-QuantizedGraph<float>(
+QuantizedGraph(
     size_t num, size_t dim, size_t max_deg,
     MetricType metric_type = METRIC_L2,
     RotatorType rotator_type = RotatorType::FhtKacRotator,
     size_t quantization_bits = 0
 );
 QGBuilder(
-    QuantizedGraph<float>& index, uint32_t ef_build, const float* data,
+    QuantizedGraph& index, uint32_t ef_build, const float* data,
     size_t num_threads = std::numeric_limits<size_t>::max(),
     QGInitialization init = QGInitialization::PiPNN
 );
@@ -64,7 +66,7 @@ The builder handles initialization internally.
 using namespace rabitqlib::symqg;
 
 // data contains rows * cols floats.
-QuantizedGraph<float> qg(rows, cols, 32);
+QuantizedGraph qg(rows, cols, 32);
 {
     QGBuilder builder(qg, 200, data.data(), 32, QGInitialization::PiPNN);
     builder.build();
@@ -97,7 +99,7 @@ C++ search accepts one vector in the original input dimension and writes `k` IDs
 and distances:
 
 ```cpp
-QuantizedGraph<float> qg;
+QuantizedGraph qg;
 qg.load("qg_example.index");
 qg.set_ef(100);
 

include/rabitqlib/index/estimator.hpp

@@ -3,12 +3,13 @@
 #include <array>
 #include <cstddef>
 #include <cstdint>
+#include <type_traits>
 
 #include "rabitqlib/defines.hpp"
 #include "rabitqlib/fastscan/fastscan.hpp"
-#include "rabitqlib/fastscan/highacc_fastscan.hpp"
 #include "rabitqlib/index/query.hpp"
 #include "rabitqlib/quantization/data_layout.hpp"
+#include "rabitqlib/simd/estimator_dispatch.hpp"
 #include "rabitqlib/utils/space.hpp"
 #include "rabitqlib/utils/warmup_space.hpp"
 
@@ -33,74 +34,9 @@ inline void split_batch_estdist(
     float* ip_x0_qr,
     bool use_hacc
 ) {
-    constexpr size_t kSafeChunkDim = 1024;
-    ConstBatchDataMap<float> cur_batch(batch_data, padded_dim);
-    std::array<int32_t, fastscan::kBatchSize> accu_values{};
-    RowMajorArrayMap<int32_t> accu_arr(accu_values.data(), 1, fastscan::kBatchSize);
-    const auto* codes_ptr = cur_batch.bin_code();
-    const auto* lut_ptr = q_obj.lut();
-    if (use_hacc) {
-        std::array<int32_t, fastscan::kBatchSize> accu_res;
-        size_t remaining_dim = padded_dim;
-
-        while (remaining_dim > kSafeChunkDim) {
-            fastscan::accumulate_hacc(codes_ptr, lut_ptr, accu_res.data(), kSafeChunkDim);
-            codes_ptr += kSafeChunkDim << 2;
-            lut_ptr += kSafeChunkDim << 3;
-            for (size_t i = 0; i < fastscan::kBatchSize; ++i) {
-                accu_arr.data()[i] += accu_res[i];
-            }
-            remaining_dim -= kSafeChunkDim;
-        }
-
-        fastscan::accumulate_hacc(codes_ptr, lut_ptr, accu_res.data(), remaining_dim);
-        for (size_t i = 0; i < fastscan::kBatchSize; ++i) {
-            accu_arr.data()[i] += accu_res[i];
-        }
-    } else {
-        std::array<uint16_t, fastscan::kBatchSize> accu_res;
-        size_t remaining_dim = padded_dim;
-
-        while (remaining_dim > kSafeChunkDim) {
-            fastscan::accumulate(codes_ptr, lut_ptr, accu_res.data(), kSafeChunkDim);
-            codes_ptr += kSafeChunkDim << 2;
-            lut_ptr += kSafeChunkDim << 2;
-            for (size_t i = 0; i < fastscan::kBatchSize; ++i) {
-                accu_arr.data()[i] += accu_res[i];
-            }
-            remaining_dim -= kSafeChunkDim;
-        }
-
-        fastscan::accumulate(codes_ptr, lut_ptr, accu_res.data(), remaining_dim);
-        for (size_t i = 0; i < fastscan::kBatchSize; ++i) {
-            accu_arr.data()[i] += accu_res[i];
-        }
-    }
-
-    std::array<float, fastscan::kBatchSize> f_add_values;
-    std::array<float, fastscan::kBatchSize> f_rescale_values;
-    std::array<float, fastscan::kBatchSize> f_error_values;
-    cur_batch.f_add().copy_to(f_add_values.data(), f_add_values.size());
-    cur_batch.f_rescale().copy_to(f_rescale_values.data(), f_rescale_values.size());
-    cur_batch.f_error().copy_to(f_error_values.data(), f_error_values.size());
-    ConstRowMajorArrayMap<float> f_add_arr(f_add_values.data(), 1, fastscan::kBatchSize);
-    ConstRowMajorArrayMap<float> f_rescale_arr(
-        f_rescale_values.data(), 1, fastscan::kBatchSize
-    );
-    ConstRowMajorArrayMap<float> f_error_arr(
-        f_error_values.data(), 1, fastscan::kBatchSize
+    simd::split_batch_estdist(
+        batch_data, q_obj, padded_dim, est_distance, low_distance, ip_x0_qr, use_hacc
     );
-
-    RowMajorArrayMap<float> est_dist_arr(est_distance, 1, fastscan::kBatchSize);
-    RowMajorArrayMap<float> ip_x0_qr_arr(ip_x0_qr, 1, fastscan::kBatchSize);
-    RowMajorArrayMap<float> low_dist_arr(low_distance, 1, fastscan::kBatchSize);
-
-    ip_x0_qr_arr = q_obj.delta() * (accu_arr.template cast<float>()) + q_obj.sum_vl_lut();
-
-    est_dist_arr =
-        f_add_arr + q_obj.g_add() + f_rescale_arr * (ip_x0_qr_arr + q_obj.k1xsumq());
-
-    low_dist_arr = est_dist_arr - f_error_arr * q_obj.g_error();
 }
 
 /**
@@ -178,6 +114,11 @@ template <typename T, typename TA = uint16_t>
 inline void qg_batch_estdist(
     const char* batch_data, const BatchQuery<T>& q_obj, size_t padded_dim, T* est_distance
 ) {
+    if constexpr (std::is_same_v<T, float> && std::is_same_v<TA, uint16_t>) {
+        simd::qg_batch_estdist(batch_data, q_obj, padded_dim, est_distance);
+        return;
+    }
+
     // Each 4-dimensional codebook can contribute at most 255, so 1024 dimensions
     // produce at most 255 * (1024 / 4) = 65280 in the uint16_t FastScan result.
     constexpr size_t kSafeChunkDim = 1024;

include/rabitqlib/index/hnsw/hnsw.hpp

@@ -27,8 +27,8 @@
 #include "rabitqlib/index/query.hpp"
 #include "rabitqlib/quantization/data_layout.hpp"
 #include "rabitqlib/quantization/rabitq.hpp"
+#include "rabitqlib/simd/hnsw_dispatch.hpp"
 #include "rabitqlib/utils/buffer.hpp"
-#include "rabitqlib/utils/cpu_features.hpp"
 #include "rabitqlib/utils/memory.hpp"
 #include "rabitqlib/utils/rotator.hpp"
 #include "rabitqlib/utils/space.hpp"
@@ -44,22 +44,6 @@ using maxheap = std::priority_queue<T>;
 template <typename T>
 using minheap = std::priority_queue<T, std::vector<T>, std::greater<T>>;
 
-class HierarchicalNSW;
-
-namespace detail {
-
-maxheap<std::pair<float, PID>> search_knn_avx2(HierarchicalNSW&, const float*, size_t);
-
-maxheap<std::pair<float, PID>> search_knn_avx512_core(
-    HierarchicalNSW&, const float*, size_t
-);
-
-maxheap<std::pair<float, PID>> search_knn_avx512_popcnt(
-    HierarchicalNSW&, const float*, size_t
-);
-
-}  // namespace detail
-
 class HierarchicalNSW {
    public:
     explicit HierarchicalNSW(){};
@@ -1108,17 +1092,7 @@ inline std::vector<std::vector<std::pair<float, PID>>> HierarchicalNSW::search(
 inline maxheap<std::pair<float, PID>> HierarchicalNSW::search_knn(
     const float* rotated_query, size_t TOPK
 ) {
-    if (rabitqlib::cpu::has_avx512_popcnt()) {
-        return detail::search_knn_avx512_popcnt(*this, rotated_query, TOPK);
-    }
-    if (rabitqlib::cpu::has_avx512_core() && rabitqlib::cpu::has_avx2()) {
-        return detail::search_knn_avx512_core(*this, rotated_query, TOPK);
-    }
-    if (rabitqlib::cpu::has_avx2()) {
-        return detail::search_knn_avx2(*this, rotated_query, TOPK);
-    }
-
-    throw std::runtime_error("HNSW search requires AVX2/FMA or AVX512 support");
+    return detail::search_knn(*this, rotated_query, TOPK);
 }
 
 template <class Kernel>

include/rabitqlib/index/ivf/initializer.hpp

@@ -172,12 +172,33 @@ class FlatInitializer : public Initializer {
     }
 };
 
+// Keep centroid routing on the same runtime-selected distance kernels as flat IVF.
+class CentroidL2Space : public hnswlib::SpaceInterface<float> {
+   private:
+    size_t dim_;
+
+    static float distance(const void* a, const void* b, const void* dim) {
+        return euclidean_sqr(
+            static_cast<const float*>(a),
+            static_cast<const float*>(b),
+            *static_cast<const size_t*>(dim)
+        );
+    }
+
+   public:
+    explicit CentroidL2Space(size_t dim) : dim_(dim) {}
+
+    size_t get_data_size() override { return dim_ * sizeof(float); }
+    hnswlib::DISTFUNC<float> get_dist_func() override { return distance; }
+    void* get_dist_func_param() override { return &dim_; }
+};
+
 class HNSWInitializer : public Initializer {
    private:
     int M_ = 16;
     int ef_construction_ = 400;
     hnswlib::HierarchicalNSW<float>* alg_hnsw_ = nullptr;
-    hnswlib::L2Space space_;
+    CentroidL2Space space_;
 
    public:
     explicit HNSWInitializer(size_t d, size_t k) : Initializer(d, k), space_(d) {

include/rabitqlib/index/ivf/ivf.hpp

@@ -3,6 +3,7 @@
 #include <algorithm>
 #include <array>
 #include <cassert>
+#include <cmath>
 #include <cstddef>
 #include <cstdint>
 #include <cstring>
@@ -29,6 +30,10 @@
 #include "rabitqlib/utils/tools.hpp"
 
 namespace rabitqlib::ivf {
+namespace detail {
+void insert_candidates(buffer::SearchBuffer<float>&, const PID*, const float*, size_t);
+}  // namespace detail
+
 class IVF {
    private:
     using ByteStorage =
@@ -717,11 +722,7 @@ inline void IVF::scan_one_batch(
 
     // Without reranking data, return the one-bit estimates directly.
     if (ex_bits_ == 0 && !raw_reranking_) {
-        for (size_t i = 0; i < num_points; ++i) {
-            PID id = ids[i];
-            float ex_dist = est_distance[i];
-            knns.insert(id, ex_dist);
-        }
+        detail::insert_candidates(knns, ids, est_distance.data(), num_points);
         return;
     }