diff --git a/include/rabitqlib/clustering/detail/kmeans.hpp b/include/rabitqlib/clustering/detail/kmeans.hpp
index 2c7f835a..77ff701f 100644
--- a/include/rabitqlib/clustering/detail/kmeans.hpp
+++ b/include/rabitqlib/clustering/detail/kmeans.hpp
@@ -9,7 +9,6 @@
 #include <cstddef>
 #include <cstdint>
 #include <cstring>
-#include <exception>
 #include <functional>
 #include <iostream>
 #include <limits>
@@ -22,6 +21,7 @@
 
 #include "rabitqlib/core/detail/simd/matrix_dispatch.hpp"
 #include "rabitqlib/core/rotator.hpp"
+#include "rabitqlib/index/detail/parallel_exception.hpp"
 #include "rabitqlib/index/types.hpp"
 #include "rabitqlib/utils/cpu_features.hpp"
 #include "rabitqlib/utils/space.hpp"
@@ -86,7 +86,7 @@ inline void exact_assign(
     std::vector<std::vector<float>> products(
         num_workers, std::vector<float>(point_block * std::min(k, kCentroidBlock))
     );
-    std::vector<std::exception_ptr> errors(num_workers);
+    index_detail::ParallelException failure;
 
     const auto exact_distance = [d, spherical](const float* vector, const float* centroid) {
         double distance = 0.0;
@@ -136,7 +136,7 @@ inline void exact_assign(
 #pragma omp for schedule(static)
         for (std::ptrdiff_t block = 0; block < static_cast<std::ptrdiff_t>(blocks);
              ++block) {
-            if (errors[thread]) {
+            if (failure.failed()) {
                 continue;
             }
             try {
@@ -201,14 +201,10 @@ inline void exact_assign(
                 for (size_t row = 0; row < rows; ++row) {
                     distances[first_point + row] = static_cast<float>(best_distances[row]);
                 }
-            } catch (...) { errors[thread] = std::current_exception(); }
-        }
-    }
-    for (const auto& error : errors) {
-        if (error) {
-            std::rethrow_exception(error);
+            } catch (...) { failure.capture(); }
         }
     }
+    failure.rethrow();
 }
 
 template <typename Parameters, typename Assignment>
diff --git a/include/rabitqlib/clustering/qgkmeans.hpp b/include/rabitqlib/clustering/qgkmeans.hpp
index c25c21e8..58b6400f 100644
--- a/include/rabitqlib/clustering/qgkmeans.hpp
+++ b/include/rabitqlib/clustering/qgkmeans.hpp
@@ -5,7 +5,6 @@
 #include <algorithm>
 #include <cstddef>
 #include <cstdint>
-#include <exception>
 #include <limits>
 #include <memory>
 #include <numeric>
@@ -17,6 +16,7 @@
 #include "rabitqlib/core/rotator.hpp"
 #include "rabitqlib/fastscan/pack_codes.hpp"
 #include "rabitqlib/fastscan/query.hpp"
+#include "rabitqlib/index/detail/parallel_exception.hpp"
 #include "rabitqlib/index/detail/visited_set.hpp"
 #include "rabitqlib/index/symqg/detail/search_scratch.hpp"
 #include "rabitqlib/index/symqg/qg.hpp"
@@ -256,7 +256,7 @@ class QGAssigner {
                     static_cast<uint32_t>(point);
             }
         }
-        std::vector<std::exception_ptr> errors(num_workers_);
+        index_detail::ParallelException failure;
 #pragma omp parallel for num_threads(num_workers_) schedule(dynamic, 64)
         for (std::ptrdiff_t loop_index = 0; loop_index < static_cast<std::ptrdiff_t>(n);
              ++loop_index) {
@@ -264,7 +264,7 @@ class QGAssigner {
                                  ? static_cast<size_t>(loop_index)
                                  : query_order[static_cast<size_t>(loop_index)];
             const int thread_id = omp_get_thread_num();
-            if (errors[static_cast<size_t>(thread_id)]) {
+            if (failure.failed()) {
                 continue;
             }
             try {
@@ -300,16 +300,10 @@ class QGAssigner {
                     distances[i] =
                         distance_func(x + i * dim_, centroids + labels[i] * dim_, dim_);
                 }
-            } catch (...) {
-                errors[static_cast<size_t>(thread_id)] = std::current_exception();
-            }
+            } catch (...) { failure.capture(); }
         }
 
-        for (const auto& error : errors) {
-            if (error) {
-                std::rethrow_exception(error);
-            }
-        }
+        failure.rethrow();
 
         previous_labels_.assign(labels, labels + n);
         query_sums.complete();
diff --git a/include/rabitqlib/clustering/rabitqkmeans.hpp b/include/rabitqlib/clustering/rabitqkmeans.hpp
index 5f4d4296..40bdcb49 100644
--- a/include/rabitqlib/clustering/rabitqkmeans.hpp
+++ b/include/rabitqlib/clustering/rabitqkmeans.hpp
@@ -6,7 +6,6 @@
 #include <array>
 #include <cstddef>
 #include <cstdint>
-#include <exception>
 #include <limits>
 #include <stdexcept>
 #include <vector>
@@ -19,6 +18,7 @@
 #include "rabitqlib/fastscan/estimator.hpp"
 #include "rabitqlib/fastscan/pack_codes.hpp"
 #include "rabitqlib/fastscan/query.hpp"
+#include "rabitqlib/index/detail/parallel_exception.hpp"
 #include "rabitqlib/index/types.hpp"
 #include "rabitqlib/utils/space.hpp"
 
@@ -126,12 +126,12 @@ class RaBitQAssigner {
         std::vector<std::vector<float>> rotated_blocks(
             num_workers_, std::vector<float>(fastscan::kBatchSize * padded_dim_)
         );
-        std::vector<std::exception_ptr> errors(num_workers_);
+        index_detail::ParallelException failure;
 #pragma omp parallel for num_threads(num_workers_) schedule(static)
         for (std::ptrdiff_t batch = 0; batch < static_cast<std::ptrdiff_t>(num_batches_);
              ++batch) {
             const auto thread = static_cast<size_t>(omp_get_thread_num());
-            if (errors[thread]) {
+            if (failure.failed()) {
                 continue;
             }
             try {
@@ -153,13 +153,9 @@ class RaBitQAssigner {
                         block * QGBatchDataMap<float>::data_bytes(padded_dim_),
                     metric_type_
                 );
-            } catch (...) { errors[thread] = std::current_exception(); }
-        }
-        for (const auto& error : errors) {
-            if (error) {
-                std::rethrow_exception(error);
-            }
+            } catch (...) { failure.capture(); }
         }
+        failure.rethrow();
     }
 
     void assign(const float* centroids, PID* labels, float* distances, bool = false) {
@@ -175,12 +171,11 @@ class RaBitQAssigner {
             query.set_g_add(distance_func(centroid, center_.data(), dim_));
         }
 
-        std::vector<std::exception_ptr> errors(num_workers_);
+        index_detail::ParallelException failure;
 #pragma omp parallel for num_threads(num_workers_) schedule(static)
         for (std::ptrdiff_t batch = 0; batch < static_cast<std::ptrdiff_t>(num_batches_);
              ++batch) {
-            const auto thread = static_cast<size_t>(omp_get_thread_num());
-            if (errors[thread]) {
+            if (failure.failed()) {
                 continue;
             }
             try {
@@ -218,13 +213,9 @@ class RaBitQAssigner {
                         data_ + id * dim_, centroids + best_labels[point] * dim_, dim_
                     );
                 }
-            } catch (...) { errors[thread] = std::current_exception(); }
-        }
-        for (const auto& error : errors) {
-            if (error) {
-                std::rethrow_exception(error);
-            }
+            } catch (...) { failure.capture(); }
         }
+        failure.rethrow();
     }
 };
 
diff --git a/include/rabitqlib/core/detail/rabitq_impl.hpp b/include/rabitqlib/core/detail/rabitq_impl.hpp
index 359f2aee..0fa4f4ee 100644
--- a/include/rabitqlib/core/detail/rabitq_impl.hpp
+++ b/include/rabitqlib/core/detail/rabitq_impl.hpp
@@ -19,6 +19,7 @@
 #include "rabitqlib/core/metric.hpp"
 #include "rabitqlib/core/pack_excode.hpp"
 #include "rabitqlib/core/types.hpp"
+#include "rabitqlib/utils/bitops.hpp"
 #include "rabitqlib/utils/openmp.hpp"  // IWYU pragma: keep (configuration macros)
 #include "rabitqlib/utils/space.hpp"
 
@@ -90,6 +91,29 @@ inline CodeFactorTerms<T> code_factor_terms(
     throw std::invalid_argument(invalid_metric_message);
 }
 
+// Rescale factor for a code centered at its own width, given code_factor_terms() of
+// that centered code. Normalized ex-bit codes use their own rescale expression.
+template <typename T>
+inline T direct_rescale(T l2_sqr, T residual_code_ip, MetricType metric_type) {
+    return l2_sqr == 0                ? T{0}
+           : metric_type == METRIC_L2 ? -2 * l2_sqr / residual_code_ip
+                                      : -l2_sqr / residual_code_ip;
+}
+
+/**
+ * Inverse of pack_binary_to_bytes<uint64_t>: write each coordinate's sign bit (0 or
+ * 1) to signs. Full 64-coordinate words come first, followed by a final 32-bit word.
+ */
+inline void unpack_sign_bits(const uint8_t* packed, size_t dim, uint16_t* signs) {
+    for (size_t start = 0; start < dim; start += 64) {
+        const size_t width = std::min(size_t{64}, dim - start);
+        const uint64_t word = bitops::load_word(packed + start / 8, width / 8);
+        for (size_t i = 0; i < width; ++i) {
+            signs[start + i] = static_cast<uint16_t>((word >> (width - 1 - i)) & 1U);
+        }
+    }
+}
+
 namespace one_bit {
 
 /**
@@ -147,9 +171,7 @@ inline void one_bit_code_with_factor_into(
         "Unsupported metric type in quantization"
     );
     f_add = factors.add;
-    f_rescale = l2_sqr == 0                ? T{0}
-                : metric_type == METRIC_L2 ? -2 * l2_sqr / factors.residual_code_ip
-                                           : -l2_sqr / factors.residual_code_ip;
+    f_rescale = direct_rescale(l2_sqr, factors.residual_code_ip, metric_type);
     f_error = factors.error;
 }
 
@@ -976,9 +998,7 @@ inline void code_factors(
         "Unsupported metric type in code_factors()"
     );
     f_add = factors.add;
-    f_rescale = l2_sqr == 0                ? T{0}
-                : metric_type == METRIC_L2 ? -2 * l2_sqr / factors.residual_code_ip
-                                           : -l2_sqr / factors.residual_code_ip;
+    f_rescale = direct_rescale(l2_sqr, factors.residual_code_ip, metric_type);
     f_error = factors.error;
 }
 
diff --git a/include/rabitqlib/index/detail/exact_nn.hpp b/include/rabitqlib/index/detail/exact_nn.hpp
index 00362a3b..1a6f29c8 100644
--- a/include/rabitqlib/index/detail/exact_nn.hpp
+++ b/include/rabitqlib/index/detail/exact_nn.hpp
@@ -20,14 +20,7 @@ inline PID exact_nn(
     T (*dist_func)(const T*, const T*, size_t)
 ) {
 #if RABITQ_USE_OPENMP
-    const auto num_workers = static_cast<int>(std::max<size_t>(
-        1,
-        std::min(
-            {resolve_num_threads(num_threads),
-             std::max<size_t>(num_points, 1),
-             static_cast<size_t>(std::numeric_limits<int>::max())}
-        )
-    ));
+    const int num_workers = detail::openmp_team_size(num_threads, num_points);
 #else
     (void)num_threads;
     constexpr int num_workers = 1;
diff --git a/include/rabitqlib/index/detail/parallel_for.hpp b/include/rabitqlib/index/detail/parallel_for.hpp
index 4d5d6af1..576326e2 100644
--- a/include/rabitqlib/index/detail/parallel_for.hpp
+++ b/include/rabitqlib/index/detail/parallel_for.hpp
@@ -8,6 +8,7 @@
 #include <thread>
 #include <vector>
 
+#include "rabitqlib/index/detail/parallel_exception.hpp"
 #include "rabitqlib/utils/openmp.hpp"
 #include "rabitqlib/utils/tools.hpp"
 
@@ -37,21 +38,17 @@ inline void parallel_for(
             // its existing scheduling. Binding requires OpenMP placement because
             // native threads inherit only the caller's place. Keep unbound nested
             // calls on the existing native-thread path.
-            std::vector<std::exception_ptr> errors(num_threads);
+            ParallelException failure;
 #pragma omp parallel for num_threads(num_threads) schedule(dynamic)
             for (size_t id = start; id < end; ++id) {
-                const size_t worker = static_cast<size_t>(omp_get_thread_num());
-                if (!errors[worker]) {
-                    try {
-                        fn(id, worker);
-                    } catch (...) { errors[worker] = std::current_exception(); }
-                }
-            }
-            for (const auto& error : errors) {
-                if (error) {
-                    std::rethrow_exception(error);
+                if (failure.failed()) {
+                    continue;
                 }
+                try {
+                    fn(id, static_cast<size_t>(omp_get_thread_num()));
+                } catch (...) { failure.capture(); }
             }
+            failure.rethrow();
             return;
         }
 #endif
@@ -113,4 +110,24 @@ inline void parallel_for(
     }
 }
 
+#if RABITQ_USE_OPENMP
+// Runs fn(i, worker) for i in [0, count) on an OpenMP team with dynamic scheduling.
+// After an iteration throws, the remaining ones are skipped; the first exception is
+// rethrown once the team has joined.
+template <class Function>
+inline void omp_for_dynamic(size_t count, size_t num_workers, Function fn) {
+    ParallelException failure;
+#pragma omp parallel for num_threads(num_workers) schedule(dynamic)
+    for (std::ptrdiff_t index = 0; index < static_cast<std::ptrdiff_t>(count); ++index) {
+        if (failure.failed()) {
+            continue;
+        }
+        try {
+            fn(static_cast<size_t>(index), static_cast<size_t>(omp_get_thread_num()));
+        } catch (...) { failure.capture(); }
+    }
+    failure.rethrow();
+}
+#endif
+
 }  // namespace rabitqlib::index_detail
diff --git a/include/rabitqlib/index/hnsw/hnsw.hpp b/include/rabitqlib/index/hnsw/hnsw.hpp
index 5d9554a6..a72ccae4 100644
--- a/include/rabitqlib/index/hnsw/hnsw.hpp
+++ b/include/rabitqlib/index/hnsw/hnsw.hpp
@@ -438,7 +438,12 @@ class HierarchicalNSW {
 
     void add_point(PID, PID, const quant::RabitqConfig&);
 
-    maxheap<std::pair<float, PID>> search_base_layer(PID, PID, int);
+    // Insertion steps shared by raw and quantized scoring. distance(id) scores an
+    // existing node against the node being inserted; prefetch(id) may warm its data.
+    template <class Distance>
+    PID greedy_descend(PID, int from_level, int to_level, Distance);
+    template <class Distance, class Prefetch>
+    maxheap<std::pair<float, PID>> search_layer(PID, int, Distance, Prefetch);
 
     PID mutually_connect_new_element(PID, maxheap<std::pair<float, PID>>&, int);
 
@@ -451,9 +456,6 @@ class HierarchicalNSW {
     void reconstruct_rotated(PID, float*) const;
     float get_quant_dist(PID target, PID query, detail::QuantQuery&) const;
     void add_point_quant(PID, int, LinkStorage&, detail::QuantQuery&);
-    maxheap<std::pair<float, PID>> search_base_layer_quant(
-        PID, PID, int, detail::QuantQuery&
-    );
     PID mutually_connect_quant(
         PID, maxheap<std::pair<float, PID>>&, int, detail::QuantQuery&
     );
diff --git a/include/rabitqlib/utils/space.hpp b/include/rabitqlib/utils/space.hpp
index 8c6ec87b..c74fc471 100644
--- a/include/rabitqlib/utils/space.hpp
+++ b/include/rabitqlib/utils/space.hpp
@@ -17,10 +17,6 @@
 #include "rabitqlib/utils/tools.hpp"
 #include "rabitqlib/utils/warmup_space.hpp"
 
-#if RABITQ_USE_OPENMP
-#include <limits>
-#endif
-
 namespace rabitqlib {
 
 inline uint64_t load_unaligned_u64(const void* data) {
@@ -191,14 +187,7 @@ inline std::vector<T> compute_centroid(
     const T* data, size_t num_points, size_t dim, size_t num_threads
 ) {
 #if RABITQ_USE_OPENMP
-    const auto thread_count = static_cast<int>(std::max<size_t>(
-        1,
-        std::min(
-            {resolve_num_threads(num_threads),
-             std::max<size_t>(num_points, 1),
-             static_cast<size_t>(std::numeric_limits<int>::max())}
-        )
-    ));
+    const int thread_count = detail::openmp_team_size(num_threads, num_points);
 #else
     (void)num_threads;
     constexpr int thread_count = 1;
diff --git a/include/rabitqlib/utils/tools.hpp b/include/rabitqlib/utils/tools.hpp
index a36a847e..cbe4dab5 100644
--- a/include/rabitqlib/utils/tools.hpp
+++ b/include/rabitqlib/utils/tools.hpp
@@ -5,6 +5,7 @@
 #include <cstddef>
 #include <cstdint>
 #include <functional>
+#include <limits>
 #include <queue>
 #include <random>
 #include <thread>
@@ -93,6 +94,20 @@ inline size_t resolve_num_threads(size_t requested) {
     return requested == 0 ? available : std::min(requested, available);
 }
 
+namespace detail {
+// OpenMP team size for work_items items: at least one and representable as int.
+inline int openmp_team_size(size_t requested, size_t work_items) {
+    return static_cast<int>(std::max<size_t>(
+        1,
+        std::min(
+            {resolve_num_threads(requested),
+             work_items,
+             static_cast<size_t>(std::numeric_limits<int>::max())}
+        )
+    ));
+}
+}  // namespace detail
+
 template <typename T, typename TP>
 float distance_ratio(
     const T* data,
diff --git a/python_bindings/bindings_common.hpp b/python_bindings/bindings_common.hpp
index 85dc1783..123198a6 100644
--- a/python_bindings/bindings_common.hpp
+++ b/python_bindings/bindings_common.hpp
@@ -4,6 +4,7 @@
 #include <pybind11/numpy.h>
 #include <pybind11/pytypes.h>
 
+#include <algorithm>
 #include <cstddef>
 #include <cstdint>
 #include <cstring>
@@ -14,7 +15,6 @@
 #include <vector>
 
 #include "rabitqlib/core/metric.hpp"
-#include "rabitqlib/core/rotator.hpp"
 #include "rabitqlib/index/types.hpp"
 #include "rabitqlib/utils/packed.hpp"
 
@@ -139,19 +139,9 @@ inline std::string metric_to_string(rabitqlib::MetricType metric) {
     return metric == rabitqlib::METRIC_IP ? "ip" : "l2";
 }
 
-inline rabitqlib::RotatorType rotator_from_string(const std::string& method) {
-    if (method == "matrix") {
-        return rabitqlib::RotatorType::MatrixRotator;
-    }
-    if (method == "fht_kac" || method == "fht") {
-        return rabitqlib::RotatorType::FhtKacRotator;
-    }
-    throw std::invalid_argument("Unsupported rotator method. Use 'fht_kac' or 'matrix'.");
-}
-
 template <typename T>
-inline py::array_t<T, py::array::c_style | py::array::forcecast> ensure_2d_array(
-    py::handle value, const char* name
+inline py::array_t<T, py::array::c_style | py::array::forcecast> ensure_array_rank(
+    py::handle value, const char* name, py::ssize_t rank
 ) {
     auto array = py::array_t<T, py::array::c_style | py::array::forcecast>::ensure(value);
     if (!array) {
@@ -161,11 +151,20 @@ inline py::array_t<T, py::array::c_style | py::array::forcecast> ensure_2d_array
             array_input_type(value)
         );
     }
-    if (array.ndim() != 2) {
+    if (array.ndim() != rank) {
         throw std::invalid_argument(
-            std::string(name) + " must be a 2D NumPy array; got shape " + array_shape(array)
+            std::string(name) + " must be a " + std::to_string(rank) +
+            "D NumPy array; got shape " + array_shape(array)
         );
     }
+    return array;
+}
+
+template <typename T>
+inline py::array_t<T, py::array::c_style | py::array::forcecast> ensure_2d_array(
+    py::handle value, const char* name
+) {
+    auto array = ensure_array_rank<T>(value, name, 2);
     // C-contiguous NumPy views may start at an unaligned byte offset. Native
     // scalar accesses require alignof(T), independently of SIMD load alignment.
     if (reinterpret_cast<std::uintptr_t>(array.data()) % alignof(T) != 0) {
@@ -186,20 +185,28 @@ template <typename T>
 inline py::array_t<T, py::array::c_style | py::array::forcecast> ensure_1d_array(
     py::handle value, const char* name
 ) {
-    auto array = py::array_t<T, py::array::c_style | py::array::forcecast>::ensure(value);
-    if (!array) {
-        throw std::invalid_argument(
-            std::string(name) + " must be convertible to a NumPy array of dtype " +
-            py::str(py::dtype::of<T>()).cast<std::string>() + "; got " +
-            array_input_type(value)
-        );
-    }
-    if (array.ndim() != 1) {
-        throw std::invalid_argument(
-            std::string(name) + " must be a 1D NumPy array; got shape " + array_shape(array)
-        );
-    }
-    return array;
+    return ensure_array_rank<T>(value, name, 1);
+}
+
+struct SearchResults {
+    py::array_t<PID> ids;
+    py::array_t<float> distances;
+};
+
+// (num_queries, k) outputs. Slots without a neighbor keep kPidMax and +inf.
+inline SearchResults unfilled_search_results(size_t num_queries, size_t k) {
+    const auto shape = std::vector<py::ssize_t>{
+        static_cast<py::ssize_t>(num_queries), static_cast<py::ssize_t>(k)};
+    SearchResults results{py::array_t<PID>(shape), py::array_t<float>(shape)};
+    std::fill(
+        results.ids.mutable_data(), results.ids.mutable_data() + results.ids.size(), kPidMax
+    );
+    std::fill(
+        results.distances.mutable_data(),
+        results.distances.mutable_data() + results.distances.size(),
+        std::numeric_limits<float>::infinity()
+    );
+    return results;
 }
 
 // Like ensure_1d_array<int64_t>, but refuses input whose dtype is not an integer. forcecast
diff --git a/python_bindings/hnsw_bindings.cpp b/python_bindings/hnsw_bindings.cpp
index 4d38191a..889e356f 100644
--- a/python_bindings/hnsw_bindings.cpp
+++ b/python_bindings/hnsw_bindings.cpp
@@ -7,7 +7,6 @@
 
 #include <algorithm>
 #include <cstddef>
-#include <limits>
 #include <memory>
 #include <stdexcept>
 #include <string>
@@ -30,8 +29,8 @@ class HnswIndex {
         size_t M,
         size_t ef_construction,
         size_t nbits,
-        const std::string& metric = "l2",
-        size_t random_seed = 100
+        const std::string& metric,
+        size_t random_seed
     )
         : index_(
               std::make_unique<rabitqlib::hnsw::HierarchicalNSW>(
@@ -49,8 +48,8 @@ class HnswIndex {
         py::handle data,
         py::handle centroids,
         py::handle cluster_ids,
-        size_t num_threads = 1,
-        bool fast_quantization = false
+        size_t num_threads,
+        bool fast_quantization
     ) {
         const auto operation = access_.write();
         auto data_array = ensure_2d_array<float>(data, "data");
@@ -104,8 +103,8 @@ class HnswIndex {
     py::array_t<rabitqlib::PID> add(
         py::handle data,
         const py::object& cluster_ids,
-        size_t num_threads = 1,
-        bool fast_quantization = false
+        size_t num_threads,
+        bool fast_quantization
     ) {
         const auto operation = access_.write();
         auto data_array = ensure_2d_array<float>(data, "data");
@@ -168,7 +167,7 @@ class HnswIndex {
         return built_ ? index_->num_stored_points() : 0;
     }
 
-    py::tuple search(py::handle queries, size_t k, size_t ef = 0, size_t num_threads = 1) {
+    py::tuple search(py::handle queries, size_t k, size_t ef, size_t num_threads) {
         const auto operation = access_.read();
         auto query_array = ensure_2d_array<float>(queries, "queries");
         if (!built_) {
@@ -186,22 +185,9 @@ class HnswIndex {
             throw std::invalid_argument("k must be between 1 and max_elements");
         }
 
-        const auto shape = std::vector<py::ssize_t>{
-            static_cast<py::ssize_t>(query_array.shape(0)), static_cast<py::ssize_t>(k)};
-        auto result_ids = py::array_t<rabitqlib::PID>(shape);
-        auto distances = py::array_t<float>(shape);
-        std::fill(
-            result_ids.mutable_data(),
-            result_ids.mutable_data() + result_ids.size(),
-            rabitqlib::kPidMax
-        );
-        std::fill(
-            distances.mutable_data(),
-            distances.mutable_data() + distances.size(),
-            std::numeric_limits<float>::infinity()
-        );
-        auto result_ids_view = result_ids.mutable_unchecked<2>();
-        auto distances_view = distances.mutable_unchecked<2>();
+        auto output = unfilled_search_results(static_cast<size_t>(query_array.shape(0)), k);
+        auto result_ids_view = output.ids.mutable_unchecked<2>();
+        auto distances_view = output.distances.mutable_unchecked<2>();
 
         const auto* queries_data = query_array.data();
         const auto num_queries = static_cast<size_t>(query_array.shape(0));
@@ -222,7 +208,7 @@ class HnswIndex {
                     results[static_cast<size_t>(i)][static_cast<size_t>(j)].first;
             }
         }
-        return py::make_tuple(result_ids, distances);
+        return py::make_tuple(output.ids, output.distances);
     }
 
     void save(py::handle path) const {
diff --git a/python_bindings/ivf_bindings.cpp b/python_bindings/ivf_bindings.cpp
index ae51719e..cecb3acb 100644
--- a/python_bindings/ivf_bindings.cpp
+++ b/python_bindings/ivf_bindings.cpp
@@ -6,10 +6,8 @@
 #include <pybind11/pytypes.h>
 #include <pybind11/stl.h>  // IWYU pragma: keep; registers std::optional casters
 
-#include <algorithm>
 #include <cstddef>
 #include <cstdint>
-#include <limits>
 #include <memory>
 #include <optional>
 #include <stdexcept>
@@ -63,8 +61,8 @@ class IvfIndex {
         size_t max_elements,
         size_t num_clusters,
         size_t nbits,
-        const std::string& metric = "l2",
-        const std::string& initializer = "auto"
+        const std::string& metric,
+        const std::string& initializer
     ) {
         const auto metric_type = metric_from_string(metric);
         index_ = std::make_unique<rabitqlib::ivf::IVF>(
@@ -82,8 +80,8 @@ class IvfIndex {
         py::handle data,
         py::handle centroids,
         py::handle cluster_ids,
-        size_t num_threads = 1,
-        bool fast_quantization = false
+        size_t num_threads,
+        bool fast_quantization
     ) {
         const auto operation = access_.write();
         auto data_array = ensure_2d_array<float>(data, "data");
@@ -118,29 +116,21 @@ class IvfIndex {
         const auto* data_ptr = data_array.data();
         const auto* centroids_ptr = centroids_array.data();
         const auto* cluster_ids_ptr = cluster_ids_vec.data();
-        try {
-            py::gil_scoped_release release;
-            index_->construct(
-                data_ptr, centroids_ptr, cluster_ids_ptr, fast_quantization, num_threads
-            );
-        } catch (...) {
-            // Early failures preserve the old index; allocation can invalidate it.
-            // Stored count remains nonzero even when every point was removed.
-            built_ = index_->num_stored_points() != 0;
-            throw;
-        }
-        built_ = true;
+        py::gil_scoped_release release;
+        index_->construct(
+            data_ptr, centroids_ptr, cluster_ids_ptr, fast_quantization, num_threads
+        );
     }
 
     py::array_t<rabitqlib::PID> add(
         py::handle data,
         const py::object& cluster_ids,
-        size_t num_threads = 1,
-        bool fast_quantization = false
+        size_t num_threads,
+        bool fast_quantization
     ) {
         const auto operation = access_.write();
         auto data_array = ensure_2d_array<float>(data, "data");
-        if (!built_) {
+        if (!built()) {
             throw std::runtime_error("IvfIndex must be built or loaded before add");
         }
         ensure_dimension(
@@ -177,7 +167,7 @@ class IvfIndex {
 
     void reserve(int64_t cluster_id, int64_t num_points) {
         const auto operation = access_.write();
-        if (!built_) {
+        if (!built()) {
             throw std::runtime_error("IvfIndex must be built or loaded before reserve");
         }
         if (cluster_id < 0 || static_cast<uint64_t>(cluster_id) >= index_->num_clusters()) {
@@ -195,7 +185,7 @@ class IvfIndex {
     size_t remove(py::handle ids) {
         const auto operation = access_.write();
         auto ids_array = ensure_1d_integer_array(ids, "ids");
-        if (!built_) {
+        if (!built()) {
             throw std::runtime_error("IvfIndex must be built or loaded before remove");
         }
         auto to_remove = checked_removal_ids(ids_array, index_->num_stored_points());
@@ -207,12 +197,12 @@ class IvfIndex {
         py::handle queries,
         size_t k,
         size_t nprobe,
-        std::optional<bool> high_accuracy = std::nullopt,
-        size_t num_threads = 1
+        std::optional<bool> high_accuracy,
+        size_t num_threads
     ) {
         const auto operation = access_.read();
         auto query_array = ensure_2d_array<float>(queries, "queries");
-        if (!built_) {
+        if (!built()) {
             throw std::runtime_error("IvfIndex must be built or loaded before search");
         }
         ensure_dimension(
@@ -226,18 +216,9 @@ class IvfIndex {
         }
 
         const size_t num_queries = static_cast<size_t>(query_array.shape(0));
-        const auto shape = std::vector<py::ssize_t>{
-            static_cast<py::ssize_t>(num_queries), static_cast<py::ssize_t>(k)};
-        auto result_ids = py::array_t<rabitqlib::PID>(shape);
-        auto distances = py::array_t<float>(shape);
-        auto* ids_data = result_ids.mutable_data();
-        auto* distances_data = distances.mutable_data();
-        std::fill(ids_data, ids_data + result_ids.size(), rabitqlib::kPidMax);
-        std::fill(
-            distances_data,
-            distances_data + distances.size(),
-            std::numeric_limits<float>::infinity()
-        );
+        auto output = unfilled_search_results(num_queries, k);
+        auto* ids_data = output.ids.mutable_data();
+        auto* distances_data = output.distances.mutable_data();
 
         const auto* queries_data = query_array.data();
         {
@@ -254,13 +235,13 @@ class IvfIndex {
             );
         }
 
-        return py::make_tuple(result_ids, distances);
+        return py::make_tuple(output.ids, output.distances);
     }
 
     void save(py::handle path) const {
         const auto filename = path_string(path);
         const auto operation = access_.read();
-        if (!built_) {
+        if (!built()) {
             throw std::runtime_error("IvfIndex must be built or loaded before save");
         }
         py::gil_scoped_release release;
@@ -275,7 +256,6 @@ class IvfIndex {
             py::gil_scoped_release release;
             wrapper.index_->load(filename.c_str());
         }
-        wrapper.built_ = true;
         return wrapper;
     }
 
@@ -289,7 +269,7 @@ class IvfIndex {
                ", max_elements=" + std::to_string(index_->max_elements()) +
                ", num_clusters=" + std::to_string(index_->num_clusters()) +
                ", initializer='" + initializer_to_string(index_->initializer_type()) +
-               "', is_built=" + (built_ ? "True" : "False") + ")";
+               "', is_built=" + (built() ? "True" : "False") + ")";
     }
 
     [[nodiscard]] size_t dim() const {
@@ -318,7 +298,7 @@ class IvfIndex {
     }
     [[nodiscard]] bool is_built() const {
         const auto operation = access_.read();
-        return built_;
+        return built();
     }
     [[nodiscard]] std::string initializer() const {
         const auto operation = access_.read();
@@ -332,7 +312,10 @@ class IvfIndex {
    private:
     IvfIndex() = default;
 
-    bool built_ = false;
+    // Construction or loading succeeded. A failed construct keeps the old index, or
+    // leaves it empty if allocation had begun. Removed points remain stored.
+    [[nodiscard]] bool built() const { return index_->num_stored_points() != 0; }
+
     mutable IndexAccess access_;
     std::unique_ptr<rabitqlib::ivf::IVF> index_;
 };
diff --git a/python_bindings/kmeans_bindings.cpp b/python_bindings/kmeans_bindings.cpp
index 25410985..e27c9aa3 100644
--- a/python_bindings/kmeans_bindings.cpp
+++ b/python_bindings/kmeans_bindings.cpp
@@ -6,7 +6,6 @@
 #include <pybind11/pytypes.h>
 #include <pybind11/stl.h>  // IWYU pragma: keep (vector conversion for iteration_stats)
 
-#include <cstring>
 #include <stdexcept>
 #include <string>
 #include <vector>
@@ -19,13 +18,10 @@ namespace py = pybind11;
 
 namespace rabitqlib::python_bindings {
 
+// Without a base object, pybind11 copies the data into a new array.
 template <typename T>
 py::array_t<T> vector_to_array(const std::vector<T>& values) {
-    py::array_t<T> result(static_cast<py::ssize_t>(values.size()));
-    if (!values.empty()) {
-        std::memcpy(result.mutable_data(), values.data(), values.size() * sizeof(T));
-    }
-    return result;
+    return py::array_t<T>(static_cast<py::ssize_t>(values.size()), values.data());
 }
 
 template <typename Clustering>
@@ -35,13 +31,7 @@ py::array_t<float> centroids_to_array(const Clustering& clustering, const char*
     }
     const auto shape = std::vector<py::ssize_t>{
         static_cast<py::ssize_t>(clustering.k), static_cast<py::ssize_t>(clustering.d)};
-    py::array_t<float> result(shape);
-    std::memcpy(
-        result.mutable_data(),
-        clustering.centroids.data(),
-        clustering.centroids.size() * sizeof(float)
-    );
-    return result;
+    return py::array_t<float>(shape, clustering.centroids.data());
 }
 
 template <typename Clustering>
@@ -126,12 +116,14 @@ void bind_clustering(py::module_& m, const char* native_name, const char* public
         .def_property_readonly(
             "obj",
             [](const Clustering& clustering) {
-                std::vector<double> values;
-                values.reserve(clustering.iteration_stats.size());
+                py::array_t<double> values(
+                    static_cast<py::ssize_t>(clustering.iteration_stats.size())
+                );
+                double* out = values.mutable_data();
                 for (const auto& stats : clustering.iteration_stats) {
-                    values.push_back(stats.obj);
+                    *out++ = stats.obj;
                 }
-                return vector_to_array(values);
+                return values;
             }
         )
         .def_readonly("iteration_stats", &Clustering::iteration_stats);
diff --git a/python_bindings/symqg_bindings.cpp b/python_bindings/symqg_bindings.cpp
index 7f49e0f1..a3a298fa 100644
--- a/python_bindings/symqg_bindings.cpp
+++ b/python_bindings/symqg_bindings.cpp
@@ -29,10 +29,7 @@ namespace rabitqlib::python_bindings {
 class SymqgIndex {
    public:
     SymqgIndex(
-        size_t dim,
-        size_t max_degree,
-        const std::string& metric = "l2",
-        size_t quantization_bits = 0
+        size_t dim, size_t max_degree, const std::string& metric, size_t quantization_bits
     )
         : dim_(dim)
         , max_degree_(max_degree)
@@ -47,10 +44,7 @@ class SymqgIndex {
     }
 
     void build(
-        py::handle data,
-        size_t ef_construction,
-        size_t num_threads = 1,
-        const std::string& init = "pipnn"
+        py::handle data, size_t ef_construction, size_t num_threads, const std::string& init
     ) {
         const auto operation = access_.write();
         if (init != "random" && init != "pipnn") {
@@ -88,13 +82,12 @@ class SymqgIndex {
             builder.build();
         }
         index_ = std::move(index);
-        built_ = true;
     }
 
-    py::tuple search(py::handle queries, size_t k, size_t ef, size_t num_threads = 1) {
+    py::tuple search(py::handle queries, size_t k, size_t ef, size_t num_threads) {
         const auto operation = access_.read();
         auto query_array = ensure_2d_array<float>(queries, "queries");
-        if (!built_) {
+        if (index_ == nullptr) {
             throw std::runtime_error("SymqgIndex must be built or loaded before search");
         }
         if (k == 0 || k > index_->num_points()) {
@@ -135,7 +128,7 @@ class SymqgIndex {
     void save(py::handle path) const {
         const auto filename = path_string(path);
         const auto operation = access_.read();
-        if (!built_) {
+        if (index_ == nullptr) {
             throw std::runtime_error("SymqgIndex must be built before save");
         }
         py::gil_scoped_release release;
@@ -154,7 +147,6 @@ class SymqgIndex {
         wrapper.max_degree_ = wrapper.index_->degree_bound();
         wrapper.metric_type_ = wrapper.index_->metric_type();
         wrapper.quantization_bits_ = wrapper.index_->quantization_bits();
-        wrapper.built_ = true;
         return wrapper;
     }
 
@@ -165,7 +157,7 @@ class SymqgIndex {
                "', quantization_bits=" + std::to_string(quantization_bits_) +
                ", num_points=" + std::to_string(index_ ? index_->num_points() : 0) +
                ", max_degree=" + std::to_string(max_degree_) +
-               ", is_built=" + (built_ ? "True" : "False") + ")";
+               ", is_built=" + (index_ != nullptr ? "True" : "False") + ")";
     }
 
     [[nodiscard]] size_t dim() const {
@@ -190,7 +182,7 @@ class SymqgIndex {
     }
     [[nodiscard]] bool is_built() const {
         const auto operation = access_.read();
-        return built_;
+        return index_ != nullptr;
     }
     [[nodiscard]] std::string metric() const {
         const auto operation = access_.read();
@@ -208,7 +200,6 @@ class SymqgIndex {
     size_t max_degree_ = 0;
     rabitqlib::MetricType metric_type_ = rabitqlib::METRIC_L2;
     size_t quantization_bits_ = 0;
-    bool built_ = false;
     mutable IndexAccess access_;
     std::unique_ptr<rabitqlib::symqg::QuantizedGraph<float>> index_;
 };
diff --git a/src/core/api.cpp b/src/core/api.cpp
index a23a12a6..c3497eba 100644
--- a/src/core/api.cpp
+++ b/src/core/api.cpp
@@ -16,7 +16,6 @@
 #include "rabitqlib/core/encoder.hpp"
 #include "rabitqlib/core/metric.hpp"
 #include "rabitqlib/core/pack_excode.hpp"
-#include "rabitqlib/utils/bitops.hpp"
 #include "rabitqlib/utils/space.hpp"
 
 namespace rabitqlib::quant::v1 {
@@ -191,21 +190,14 @@ void read_codes(const Parameters& p, CodeView input, ScratchStorage& s) {
         }
         return;
     }
+    rabitq_impl::unpack_sign_bits(bytes, p.dim, s.codes);
     const size_t ex_bits = p.total_bits - 1;
     if (ex_bits > 0) {
         rabitq_impl::ex_bits::unpacking_rabitqplus_code(
             static_cast<const uint8_t*>(input.extra.data), s.extra, p.dim, ex_bits
         );
-    }
-    // Sign words retain the existing full-64 / final-32 coordinate convention.
-    for (size_t start = 0; start < p.dim; start += 64) {
-        const size_t width = std::min(size_t{64}, p.dim - start);
-        const uint64_t word = bitops::load_word(bytes + start / 8, width / 8);
-        for (size_t i = 0; i < width; ++i) {
-            const auto sign = static_cast<uint16_t>((word >> (width - 1 - i)) & 1U);
-            s.codes[start + i] = static_cast<uint16_t>(
-                (sign << ex_bits) | (ex_bits == 0 ? 0 : s.extra[start + i])
-            );
+        for (size_t i = 0; i < p.dim; ++i) {
+            s.codes[i] = static_cast<uint16_t>((s.codes[i] << ex_bits) | s.extra[i]);
         }
     }
 }
diff --git a/src/index/hnsw.cpp b/src/index/hnsw.cpp
index 8570eeb8..c51f7918 100644
--- a/src/index/hnsw.cpp
+++ b/src/index/hnsw.cpp
@@ -587,6 +587,96 @@ void HierarchicalNSW::construct(
     construction_data_ = nullptr;
 }
 
+template <class Distance>
+PID HierarchicalNSW::greedy_descend(
+    PID current_id, int from_level, int to_level, Distance distance
+) {
+    float current_distance = distance(current_id);
+    for (int level = from_level; level > to_level; level--) {
+        bool changed = true;
+        while (changed) {
+            changed = false;
+            std::unique_lock<std::mutex> lock(link_list_locks_[current_id]);
+            const PID* data = get_linklist(current_id, level);
+            const size_t size = get_list_count(data);
+            const PID* datal = data + 1;
+            for (size_t i = 0; i < size; i++) {
+                PID cand = datal[i];
+                assert(cand < num_stored_points_);
+                float d = distance(cand);
+                if (d < current_distance) {
+                    current_distance = d;
+                    current_id = cand;
+                    changed = true;
+                }
+            }
+        }
+    }
+    return current_id;
+}
+
+template <class Distance, class Prefetch>
+maxheap<std::pair<float, PID>> HierarchicalNSW::search_layer(
+    PID entry_point_id, int layer, Distance distance, Prefetch prefetch
+) {
+    auto visited = visited_list_pool_->get_free_vislist();
+
+    maxheap<std::pair<float, PID>> top_candidates;
+    minheap<std::pair<float, PID>> candidate_set;
+
+    float lower_bound = distance(entry_point_id);
+    top_candidates.emplace(lower_bound, entry_point_id);
+    candidate_set.emplace(lower_bound, entry_point_id);
+    visited->set(entry_point_id);
+
+    while (!candidate_set.empty()) {
+        std::pair<float, PID> curr_el_pair = candidate_set.top();
+        if (curr_el_pair.first > lower_bound && top_candidates.size() == ef_construction_) {
+            break;
+        }
+        candidate_set.pop();
+
+        PID current_id = curr_el_pair.second;
+        std::unique_lock<std::mutex> lock(link_list_locks_[current_id]);
+
+        const PID* data =
+            (layer == 0) ? get_linklist0(current_id) : get_linklist(current_id, layer);
+        const size_t size = get_list_count(data);
+        const PID* datal = data + 1;
+
+        prefetch(*datal);
+        prefetch(*(datal + 1));
+
+        for (size_t j = 0; j < size; j++) {
+            PID candidate_id = datal[j];
+            assert(candidate_id < num_stored_points_);
+            if (visited->get(candidate_id)) {
+                continue;
+            }
+            visited->set(candidate_id);
+
+            if (j < size - 1) {
+                prefetch(*(datal + j + 1));
+            }
+
+            float candidate_distance = distance(candidate_id);
+            if (top_candidates.size() < ef_construction_ ||
+                lower_bound > candidate_distance) {
+                candidate_set.emplace(candidate_distance, candidate_id);
+                top_candidates.emplace(candidate_distance, candidate_id);
+                if (top_candidates.size() > ef_construction_) {
+                    top_candidates.pop();
+                }
+                if (!top_candidates.empty()) {
+                    lower_bound = top_candidates.top().first;
+                }
+            }
+        }
+    }
+    visited_list_pool_->release_vis_list(std::move(visited));
+    return top_candidates;
+}
+
 void HierarchicalNSW::add_point(
     PID label, PID cluster_id, const quant::RabitqConfig& config
 ) {
@@ -656,34 +746,25 @@ void HierarchicalNSW::add_point(
     }
 
     if (static_cast<signed>(current_id) != -1) {
+        const auto distance = [this, new_node_id](PID id) {
+            return get_data_dist(id, new_node_id);
+        };
+        const auto prefetch = [this](PID id) {
+            rabitqlib::memory::mem_prefetch_l1(
+                reinterpret_cast<const char*>(
+                    construction_data_ + (get_external_label(id) * dim_)
+                ),
+                padded_dim_ / 16
+            );
+        };
         if (node_level < previous_max_level) {
-            float current_distance = get_data_dist(current_id, new_node_id);
-            for (int level = previous_max_level; level > node_level; level--) {
-                bool changed = true;
-                while (changed) {
-                    changed = false;
-                    std::unique_lock<std::mutex> lock(link_list_locks_[current_id]);
-                    const PID* data = get_linklist(current_id, level);
-                    const size_t size = get_list_count(data);
-
-                    const PID* datal = data + 1;
-                    for (size_t i = 0; i < size; i++) {
-                        PID cand = datal[i];
-                        assert(cand < num_stored_points_);
-                        float d = get_data_dist(cand, new_node_id);
-                        if (d < current_distance) {
-                            current_distance = d;
-                            current_id = cand;
-                            changed = true;
-                        }
-                    }
-                }
-            }
+            current_id =
+                greedy_descend(current_id, previous_max_level, node_level, distance);
         }
 
         for (int level = std::min(node_level, previous_max_level); level >= 0; level--) {
             maxheap<std::pair<float, PID>> top_candidates =
-                search_base_layer(current_id, new_node_id, level);
+                search_layer(current_id, level, distance, prefetch);
             current_id = mutually_connect_new_element(new_node_id, top_candidates, level);
         }
     }
@@ -695,83 +776,6 @@ void HierarchicalNSW::add_point(
     }
 }
 
-maxheap<std::pair<float, PID>> HierarchicalNSW::search_base_layer(
-    PID entry_point_id, PID new_node_id, int layer
-) {
-    auto visited = visited_list_pool_->get_free_vislist();
-
-    maxheap<std::pair<float, PID>> top_candidates;
-    minheap<std::pair<float, PID>> candidate_set;
-
-    float lower_bound = get_data_dist(entry_point_id, new_node_id);
-    top_candidates.emplace(lower_bound, entry_point_id);
-    candidate_set.emplace(lower_bound, entry_point_id);
-    visited->set(entry_point_id);
-
-    while (!candidate_set.empty()) {
-        std::pair<float, PID> curr_el_pair = candidate_set.top();
-        if (curr_el_pair.first > lower_bound && top_candidates.size() == ef_construction_) {
-            break;
-        }
-        candidate_set.pop();
-
-        PID current_id = curr_el_pair.second;
-
-        std::unique_lock<std::mutex> lock(link_list_locks_[current_id]);
-
-        const PID* data =
-            (layer == 0) ? get_linklist0(current_id) : get_linklist(current_id, layer);
-        const size_t size = get_list_count(data);
-        const PID* datal = data + 1;
-
-        rabitqlib::memory::mem_prefetch_l1(
-            reinterpret_cast<const char*>(
-                construction_data_ + (get_external_label(*datal) * dim_)
-            ),
-            padded_dim_ / 16
-        );
-
-        rabitqlib::memory::mem_prefetch_l1(
-            reinterpret_cast<const char*>(
-                construction_data_ + (get_external_label(*(datal + 1)) * dim_)
-            ),
-            padded_dim_ / 16
-        );
-
-        for (size_t j = 0; j < size; j++) {
-            PID candidate_id = *(datal + j);
-            if (visited->get(candidate_id)) {
-                continue;
-            }
-            visited->set(candidate_id);
-
-            if (j < size - 1) {
-                rabitqlib::memory::mem_prefetch_l1(
-                    reinterpret_cast<const char*>(
-                        construction_data_ + (get_external_label(*(datal + j + 1)) * dim_)
-                    ),
-                    padded_dim_ / 16
-                );
-            }
-
-            float candidate_distance = get_data_dist(candidate_id, new_node_id);
-            if (top_candidates.size() < ef_construction_ ||
-                lower_bound > candidate_distance) {
-                candidate_set.emplace(candidate_distance, candidate_id);
-                top_candidates.emplace(candidate_distance, candidate_id);
-                if (top_candidates.size() > ef_construction_) {
-                    top_candidates.pop();
-                }
-                if (!top_candidates.empty()) {
-                    lower_bound = top_candidates.top().first;
-                }
-            }
-        }
-    }
-    visited_list_pool_->release_vis_list(std::move(visited));
-    return top_candidates;
-}
-
 namespace detail {
 template <class Distance>
 void prune_neighbors(
@@ -1004,16 +1008,7 @@ void HierarchicalNSW::reconstruct_rotated(PID id, float* out) const {
     combined.resize(padded_dim_);
 
     const ConstBinDataMap<float> bin(get_bindata_by_internalid(id), padded_dim_);
-    const uint8_t* sign_code = bin.bin_code();
-    for (size_t word = 0; word * 64 < padded_dim_; ++word) {
-        uint64_t bits = 0;
-        const size_t base = word * 64;
-        const size_t span = std::min<size_t>(64, padded_dim_ - base);
-        std::memcpy(&bits, sign_code + (word * sizeof(uint64_t)), span / 8);
-        for (size_t j = 0; j < span; ++j) {
-            combined[base + j] = static_cast<uint16_t>((bits >> (span - 1 - j)) & 1ULL);
-        }
-    }
+    quant::rabitq_impl::unpack_sign_bits(bin.bin_code(), padded_dim_, combined.data());
 
     const auto* centroid = reinterpret_cast<const float*>(centroids_memory_.get()) +
                            (get_clusterid_by_internalid(id) * padded_dim_);
@@ -1090,60 +1085,6 @@ float HierarchicalNSW::get_quant_dist(
         .distance;
 }
 
-maxheap<std::pair<float, PID>> HierarchicalNSW::search_base_layer_quant(
-    PID entry_point_id, PID new_node_id, int layer, detail::QuantQuery& query
-) {
-    auto visited = visited_list_pool_->get_free_vislist();
-
-    maxheap<std::pair<float, PID>> top_candidates;
-    minheap<std::pair<float, PID>> candidate_set;
-
-    float lower_bound = get_quant_dist(entry_point_id, new_node_id, query);
-    top_candidates.emplace(lower_bound, entry_point_id);
-    candidate_set.emplace(lower_bound, entry_point_id);
-    visited->set(entry_point_id);
-
-    while (!candidate_set.empty()) {
-        std::pair<float, PID> curr_el_pair = candidate_set.top();
-        if (curr_el_pair.first > lower_bound && top_candidates.size() == ef_construction_) {
-            break;
-        }
-        candidate_set.pop();
-
-        PID current_id = curr_el_pair.second;
-        std::unique_lock<std::mutex> lock(link_list_locks_[current_id]);
-
-        const PID* data =
-            (layer == 0) ? get_linklist0(current_id) : get_linklist(current_id, layer);
-        const size_t size = get_list_count(data);
-        const PID* datal = data + 1;
-
-        for (size_t j = 0; j < size; j++) {
-            PID candidate_id = datal[j];
-            assert(candidate_id < num_stored_points_);
-            if (visited->get(candidate_id)) {
-                continue;
-            }
-            visited->set(candidate_id);
-
-            float candidate_distance = get_quant_dist(candidate_id, new_node_id, query);
-            if (top_candidates.size() < ef_construction_ ||
-                lower_bound > candidate_distance) {
-                candidate_set.emplace(candidate_distance, candidate_id);
-                top_candidates.emplace(candidate_distance, candidate_id);
-                if (top_candidates.size() > ef_construction_) {
-                    top_candidates.pop();
-                }
-                if (!top_candidates.empty()) {
-                    lower_bound = top_candidates.top().first;
-                }
-            }
-        }
-    }
-    visited_list_pool_->release_vis_list(std::move(visited));
-    return top_candidates;
-}
-
 PID HierarchicalNSW::mutually_connect_quant(
     PID new_node_id,
     maxheap<std::pair<float, PID>>& top_candidates,
@@ -1184,34 +1125,18 @@ void HierarchicalNSW::add_point_quant(
 
     try {
         if (static_cast<signed>(current_id) != -1) {
+            const auto distance = [this, new_node_id, &query](PID id) {
+                return get_quant_dist(id, new_node_id, query);
+            };
             if (node_level < previous_max_level) {
-                float current_distance = get_quant_dist(current_id, new_node_id, query);
-                for (int level = previous_max_level; level > node_level; level--) {
-                    bool changed = true;
-                    while (changed) {
-                        changed = false;
-                        std::unique_lock<std::mutex> lock(link_list_locks_[current_id]);
-                        const PID* data = get_linklist(current_id, level);
-                        const size_t size = get_list_count(data);
-                        const PID* datal = data + 1;
-                        for (size_t i = 0; i < size; i++) {
-                            PID cand = datal[i];
-                            assert(cand < num_stored_points_);
-                            float d = get_quant_dist(cand, new_node_id, query);
-                            if (d < current_distance) {
-                                current_distance = d;
-                                current_id = cand;
-                                changed = true;
-                            }
-                        }
-                    }
-                }
+                current_id =
+                    greedy_descend(current_id, previous_max_level, node_level, distance);
             }
 
             for (int level = std::min(node_level, previous_max_level); level >= 0;
                  level--) {
                 maxheap<std::pair<float, PID>> top_candidates =
-                    search_base_layer_quant(current_id, new_node_id, level, query);
+                    search_layer(current_id, level, distance, [](PID) {});
                 current_id =
                     mutually_connect_quant(new_node_id, top_candidates, level, query);
             }
diff --git a/src/index/ivf.cpp b/src/index/ivf.cpp
index 642cb5eb..e2237a0b 100644
--- a/src/index/ivf.cpp
+++ b/src/index/ivf.cpp
@@ -234,15 +234,8 @@ void IVF::construct(
     }
     /* Quantize each cluster */
     {
-        index_detail::ParallelException failure;
-#pragma omp parallel for schedule(dynamic) num_threads(num_threads)
-        for (std::ptrdiff_t index = 0; index < static_cast<std::ptrdiff_t>(num_clusters_);
-             ++index) {
-            if (failure.failed()) {
-                continue;
-            }
-            try {
-                const size_t i = static_cast<size_t>(index);
+        index_detail::omp_for_dynamic(
+            num_clusters_, num_threads, [&](size_t i, size_t thread) {
                 const float* cur_centroid = centroids + (i * dim_);
                 float* cur_rotated_c = &rotated_centroids[i * padded_dim_];
                 Cluster& cp = clusters_[i];
@@ -252,12 +245,11 @@ void IVF::construct(
                     data,
                     cur_centroid,
                     cur_rotated_c,
-                    rotated_blocks[static_cast<size_t>(omp_get_thread_num())].data(),
+                    rotated_blocks[thread].data(),
                     config
                 );
-            } catch (...) { failure.capture(); }
-        }
-        failure.rethrow();
+            }
+        );
     }
 
     this->initializer_->add_vectors(rotated_centroids.data(), num_threads);
@@ -1078,16 +1070,8 @@ void IVF::search_batch(
         }
         return;
     }
-    index_detail::ParallelException failure;
-#pragma omp parallel for num_threads(num_workers) schedule(dynamic)
-    for (std::ptrdiff_t query_index = 0;
-         query_index < static_cast<std::ptrdiff_t>(num_queries);
-         ++query_index) {
-        if (failure.failed()) {
-            continue;
-        }
-        const size_t i = static_cast<size_t>(query_index);
-        try {
+    index_detail::omp_for_dynamic(
+        num_queries, num_workers, [&](size_t i, size_t /*thread*/) {
             search(
                 queries + i * dim_,
                 k,
@@ -1096,9 +1080,8 @@ void IVF::search_batch(
                 distances == nullptr ? nullptr : distances + i * k,
                 hacc
             );
-        } catch (...) { failure.capture(); }
-    }
-    failure.rethrow();
+        }
+    );
 }
 
 inline void IVF::validate_search(
diff --git a/src/index/pipnn.cpp b/src/index/pipnn.cpp
index 9080f34d..a00c6cbd 100644
--- a/src/index/pipnn.cpp
+++ b/src/index/pipnn.cpp
@@ -21,6 +21,7 @@
 #include "rabitqlib/core/metric.hpp"
 #include "rabitqlib/index/detail/buffer.hpp"
 #include "rabitqlib/index/detail/parallel_exception.hpp"
+#include "rabitqlib/index/detail/parallel_for.hpp"
 #include "rabitqlib/index/types.hpp"
 #include "rabitqlib/utils/space.hpp"
 #include "rabitqlib/utils/tools.hpp"
@@ -412,20 +413,11 @@ std::vector<Bucket> cluster(
             // nested, forcing the runtime to create workers again on every reset.
             process_bucket(0);
         } else {
-            {
-                index_detail::ParallelException failure;
-#pragma omp parallel for num_threads(std::min(num_workers, active.size())) schedule(dynamic)
-                for (std::ptrdiff_t i = 0; i < static_cast<std::ptrdiff_t>(active.size());
-                     ++i) {
-                    if (failure.failed()) {
-                        continue;
-                    }
-                    try {
-                        process_bucket(static_cast<size_t>(i));
-                    } catch (...) { failure.capture(); }
-                }
-                failure.rethrow();
-            }
+            index_detail::omp_for_dynamic(
+                active.size(),
+                std::min(num_workers, active.size()),
+                [&](size_t i, size_t /*thread*/) { process_bucket(i); }
+            );
         }
         std::vector<BucketJob> new_active;
         for (size_t i = 0; i < active.size(); ++i) {
diff --git a/src/index/qg.cpp b/src/index/qg.cpp
index 518f751f..1615f609 100644
--- a/src/index/qg.cpp
+++ b/src/index/qg.cpp
@@ -498,11 +498,7 @@ void QuantizedGraph<float>::search_batch_with_ef(
         return;
     }
 
-    const auto num_workers = static_cast<int>(std::min(
-        {resolve_num_threads(num_threads),
-         num_queries,
-         static_cast<size_t>(std::numeric_limits<int>::max())}
-    ));
+    const int num_workers = rabitqlib::detail::openmp_team_size(num_threads, num_queries);
     index_detail::ParallelException failure;
 #pragma omp parallel num_threads(num_workers) if (num_workers > 1)
     {
diff --git a/src/index/qg_builder.cpp b/src/index/qg_builder.cpp
index 9c92c673..9d361a41 100644
--- a/src/index/qg_builder.cpp
+++ b/src/index/qg_builder.cpp
@@ -20,6 +20,7 @@
 #include "rabitqlib/index/detail/buffer.hpp"
 #include "rabitqlib/index/detail/exact_nn.hpp"
 #include "rabitqlib/index/detail/parallel_exception.hpp"
+#include "rabitqlib/index/detail/parallel_for.hpp"
 #include "rabitqlib/index/detail/visited_set.hpp"
 #include "rabitqlib/index/symqg/detail/pipnn.hpp"
 #include "rabitqlib/index/symqg/qg.hpp"
@@ -463,17 +464,10 @@ void QGBuilder::heuristic_prune(
  */
 void QGBuilder::search_new_neighbors() {
     {
-        index_detail::ParallelException failure;
-#pragma omp parallel for num_threads(num_workers_) schedule(dynamic)
-        for (std::ptrdiff_t index = 0; index < static_cast<std::ptrdiff_t>(num_points_);
-             ++index) {
-            if (failure.failed()) {
-                continue;
-            }
-            try {
-                const size_t i = static_cast<size_t>(index);
+        index_detail::omp_for_dynamic(
+            num_points_, num_workers_, [&](size_t i, size_t thread) {
                 PID point_id = i;
-                auto& worker = workers_[static_cast<size_t>(omp_get_thread_num())];
+                auto& worker = workers_[thread];
                 CandidateList& candidates = worker.candidates;
                 candidates.clear();
                 VisitedSet& visited = worker.visited;
@@ -513,9 +507,8 @@ void QGBuilder::search_new_neighbors() {
                 // prune and update qg
                 new_neighbors_[point_id].reserve(max_degree_);
                 heuristic_prune(candidates, new_neighbors_[point_id], /*recorded=*/nullptr);
-            } catch (...) { failure.capture(); }
-        }
-        failure.rethrow();
+            }
+        );
     }
 }
 
@@ -528,15 +521,9 @@ void QGBuilder::add_reverse_edges(bool final_pass) {
     // Keep new_neighbors_ read-only while reverse candidates are collected. Mutating a
     // destination row here races with another worker reading that row as its source.
     {
-        index_detail::ParallelException failure;
-#pragma omp parallel for num_threads(num_workers_) schedule(dynamic)
-        for (std::ptrdiff_t index = 0; index < static_cast<std::ptrdiff_t>(num_points_);
-             ++index) {
-            if (failure.failed()) {
-                continue;
-            }
-            try {
-                const PID data_id = static_cast<PID>(index);
+        index_detail::omp_for_dynamic(
+            num_points_, num_workers_, [&](size_t i, size_t /*thread*/) {
+                const PID data_id = static_cast<PID>(i);
                 for (const auto& nei : new_neighbors_[data_id]) {
                     const PID destination = nei.id;
                     const CandidateList& destination_neighbors =
@@ -560,21 +547,14 @@ void QGBuilder::add_reverse_edges(bool final_pass) {
                         reverse_buffer_[destination].emplace_back(data_id, nei.distance);
                     }
                 }
-            } catch (...) { failure.capture(); }
-        }
-        failure.rethrow();
+            }
+        );
     }
 
     {
-        index_detail::ParallelException failure;
-#pragma omp parallel for num_threads(num_workers_) schedule(dynamic)
-        for (std::ptrdiff_t index = 0; index < static_cast<std::ptrdiff_t>(num_points_);
-             ++index) {
-            if (failure.failed()) {
-                continue;
-            }
-            try {
-                const PID data_id = static_cast<PID>(index);
+        index_detail::omp_for_dynamic(
+            num_points_, num_workers_, [&](size_t i, size_t thread) {
+                const PID data_id = static_cast<PID>(i);
                 CandidateList& candidate_pool = reverse_buffer_[data_id];
                 if (qg_.is_quantized() && !candidate_pool.empty()) {
                     // RaBitQ estimates are directional: score destination -> source afresh.
@@ -591,7 +571,7 @@ void QGBuilder::add_reverse_edges(bool final_pass) {
                     new_neighbors_[data_id].end()
                 );
                 std::sort(candidate_pool.begin(), candidate_pool.end());
-                auto& worker = workers_[static_cast<size_t>(omp_get_thread_num())];
+                auto& worker = workers_[thread];
                 worker.pruned.clear();
                 if (final_pass) {
                     worker.pruned.reserve(kMaxPrunedSize);
@@ -606,9 +586,8 @@ void QGBuilder::add_reverse_edges(bool final_pass) {
                 if (final_pass) {
                     complete_degree(data_id, worker);
                 }
-            } catch (...) { failure.capture(); }
-        }
-        failure.rethrow();
+            }
+        );
     }
 }
 
@@ -616,15 +595,8 @@ void QGBuilder::random_init() {
     const PID min_id = 0;
     const PID max_id = num_points_ - 1;
     {
-        index_detail::ParallelException failure;
-#pragma omp parallel for num_threads(num_workers_) schedule(dynamic)
-        for (std::ptrdiff_t index = 0; index < static_cast<std::ptrdiff_t>(num_points_);
-             ++index) {
-            if (failure.failed()) {
-                continue;
-            }
-            try {
-                const size_t i = static_cast<size_t>(index);
+        index_detail::omp_for_dynamic(
+            num_points_, num_workers_, [&](size_t i, size_t thread) {
                 new_neighbors_[i].clear();
                 std::mt19937 rng(seed_ ^ static_cast<uint32_t>(i));
                 std::uniform_int_distribution<PID> distribution(min_id, max_id);
@@ -644,14 +616,9 @@ void QGBuilder::random_init() {
                     new_neighbors_[i].emplace_back(cur_neigh, query.distance_to(cur_neigh));
                 }
 
-                write_row(
-                    static_cast<PID>(i),
-                    new_neighbors_[i],
-                    workers_[static_cast<size_t>(omp_get_thread_num())]
-                );
-            } catch (...) { failure.capture(); }
-        }
-        failure.rethrow();
+                write_row(static_cast<PID>(i), new_neighbors_[i], workers_[thread]);
+            }
+        );
     }
 }
 
@@ -728,23 +695,11 @@ void QGBuilder::run_pass(bool final_pass) {
 
     // update qg
     {
-        index_detail::ParallelException failure;
-#pragma omp parallel for num_threads(num_workers_) schedule(dynamic)
-        for (std::ptrdiff_t index = 0; index < static_cast<std::ptrdiff_t>(num_points_);
-             ++index) {
-            if (failure.failed()) {
-                continue;
+        index_detail::omp_for_dynamic(
+            num_points_, num_workers_, [&](size_t i, size_t thread) {
+                write_row(static_cast<PID>(i), new_neighbors_[i], workers_[thread]);
             }
-            try {
-                const size_t i = static_cast<size_t>(index);
-                write_row(
-                    static_cast<PID>(i),
-                    new_neighbors_[i],
-                    workers_[static_cast<size_t>(omp_get_thread_num())]
-                );
-            } catch (...) { failure.capture(); }
-        }
-        failure.rethrow();
+        );
     }
 }
 void QGBuilder::BuildQuery::reset(PID id) {
diff --git a/tests/unit/rabitqlib/utils/thread_count_test.cpp b/tests/unit/rabitqlib/utils/thread_count_test.cpp
index e945269f..dd8f81d3 100644
--- a/tests/unit/rabitqlib/utils/thread_count_test.cpp
+++ b/tests/unit/rabitqlib/utils/thread_count_test.cpp
@@ -226,3 +226,53 @@ TEST(ThreadCount, NestedParallelForPreservesBindingAndWorkerIds) {
     GTEST_SKIP() << "OpenMP binding requires OpenMP 4.0";
 #endif
 }
+
+TEST(ThreadCount, OmpForDynamicVisitsEachItemOnceAndPropagatesErrors) {
+    constexpr size_t kCount = 257;
+    for (size_t workers : {1U, 4U}) {
+        SCOPED_TRACE(workers);
+        std::vector<std::atomic<size_t>> visits(kCount);
+        for (auto& visited : visits) {
+            visited.store(0);
+        }
+        rabitqlib::index_detail::omp_for_dynamic(
+            kCount, workers, [&](size_t i, size_t worker) {
+                EXPECT_LT(worker, workers);
+                visits[i].fetch_add(1);
+            }
+        );
+        for (const auto& visited : visits) {
+            EXPECT_EQ(visited.load(), 1U);
+        }
+
+        std::atomic<size_t> calls{0};
+        EXPECT_THROW(
+            rabitqlib::index_detail::omp_for_dynamic(
+                kCount,
+                workers,
+                [&](size_t i, size_t) {
+                    calls.fetch_add(1);
+                    if (i == 3) {
+                        throw std::runtime_error("failed");
+                    }
+                }
+            ),
+            std::runtime_error
+        );
+        // Items after a failure are skipped; the serial order stops at the first one.
+        if (workers == 1) {
+            EXPECT_EQ(calls.load(), 4U);
+        }
+    }
+}
+
+TEST(ThreadCount, OpenMpTeamSizeIsPositiveAndBounded) {
+    const size_t available = rabitqlib::total_threads();
+    EXPECT_EQ(rabitqlib::detail::openmp_team_size(0, 0), 1);
+    EXPECT_EQ(rabitqlib::detail::openmp_team_size(1, 100), 1);
+    EXPECT_EQ(rabitqlib::detail::openmp_team_size(0, 1), 1);
+    EXPECT_EQ(
+        rabitqlib::detail::openmp_team_size(0, std::numeric_limits<size_t>::max()),
+        static_cast<int>(std::min<size_t>(available, std::numeric_limits<int>::max()))
+    );
+}