perf(ann): unit-vector dot product and a reusable visited set

- Cosine distance was 1 - dot/(|a||b|), re-deriving both norms on every
  evaluation in the innermost loop of build and search. The index now stores
  unit vectors (prepared at build, insert, graph load and HDF5 load; the query
  once per search) and uses 1 - dot. Zero vectors stay zero, giving distance 1
  as before. Returned distances are unchanged.
- search_layer allocated a HashSet of visited nodes per call. It is now an
  epoch-stamped u32 array in thread-local scratch, reused across calls, so
  search(&self) stays shareable between threads.
- clawhdf5-accel caches the detected SIMD backend in a OnceLock.

Recall is identical. Build 2.75 -> 1.89 s (10K), ~38 -> 21 s (100K); QPS at
ef=64 22.7K -> 39K (10K), 10.4K -> 14K (100K).

Co-Authored-By: Claude Fable 5.1 <[email protected]>

crates/clawhdf5-ann/src/hnsw.rs

@@ -1,6 +1,6 @@
 //! HNSW index implementation with HDF5 serialization.
 
-use std::collections::{BinaryHeap, HashSet};
+use std::collections::BinaryHeap;
 
 use clawhdf5_format::attribute::extract_attributes_full;
 use clawhdf5_format::data_layout::DataLayout;
@@ -51,10 +51,30 @@ impl DistanceMetric {
 fn compute_distance(a: &[f32], b: &[f32], metric: DistanceMetric) -> f32 {
     match metric {
         DistanceMetric::L2 => clawhdf5_accel::l2_distance(a, b),
-        DistanceMetric::Cosine => 1.0 - clawhdf5_accel::cosine_similarity(a, b),
+        // Both sides are unit length (see `prepare`), so cosine similarity is
+        // the plain dot product. Computing it as dot / (|a| * |b|) re-derived
+        // both norms on every call — three reductions instead of one, in the
+        // innermost loop of both build and search.
+        DistanceMetric::Cosine => 1.0 - clawhdf5_accel::dot_product(a, b),
     }
 }
 
+/// Put a vector in the form the index stores and compares: unit length for the
+/// cosine metric, unchanged for L2. A zero vector stays zero, giving distance 1
+/// to everything — what the cosine kernel reports for a degenerate input.
+fn prepare(mut v: Vec<f32>, metric: DistanceMetric) -> Vec<f32> {
+    if metric == DistanceMetric::Cosine {
+        let norm = clawhdf5_accel::vector_norm(&v);
+        if norm > f32::EPSILON {
+            let inv = 1.0 / norm;
+            v.iter_mut().for_each(|x| *x *= inv);
+        } else {
+            v.iter_mut().for_each(|x| *x = 0.0);
+        }
+    }
+    v
+}
+
 /// Assign a random level to a new node based on the HNSW probability distribution.
 ///
 /// Uses a deterministic approach based on the node index for reproducibility.
@@ -204,6 +224,8 @@ impl HnswIndex {
 
         let m_max0 = m * 2;
         let n = vectors.len();
+        let prepared: Vec<Vec<f32>> = vectors.iter().map(|v| prepare(v.clone(), metric)).collect();
+        let vectors: &[Vec<f32>] = &prepared;
 
         // Assign levels to all nodes
         let mut node_levels = Vec::with_capacity(n);
@@ -288,7 +310,7 @@ impl HnswIndex {
         }
 
         Self {
-            vectors: vectors.to_vec(),
+            vectors: prepared,
             graph,
             deleted: vec![false; n],
             entry_point,
@@ -327,6 +349,7 @@ impl HnswIndex {
     /// # Panics
     /// Panics if `vector`'s dimension does not match the existing vectors.
     pub fn insert(&mut self, vector: Vec<f32>) -> usize {
+        let vector = prepare(vector, self.metric);
         let id = self.vectors.len();
 
         // Seed an empty index.
@@ -483,6 +506,8 @@ impl HnswIndex {
             "query dimension mismatch"
         );
         let ef = ef.max(k);
+        let prepared_query = prepare(query.to_vec(), self.metric);
+        let query = prepared_query.as_slice();
 
         let mut ep = self.entry_point;
         let top_layer = self.graph.len().saturating_sub(1);
@@ -633,7 +658,9 @@ impl HnswIndex {
                     actual: flat_vectors.len(),
                 });
             }
-            vectors.push(flat_vectors[start..end].to_vec());
+            // Files written before vectors were stored unit-length hold the
+            // raw ones; preparing is idempotent, so this handles both.
+            vectors.push(prepare(flat_vectors[start..end].to_vec(), metric));
         }
 
         // Read graph layers
@@ -835,7 +862,7 @@ impl HnswIndex {
         }
 
         Ok(Self {
-            vectors,
+            vectors: vectors.into_iter().map(|v| prepare(v, metric)).collect(),
             graph,
             deleted,
             entry_point,
@@ -937,9 +964,62 @@ fn search_layer(
         distance: ep_dist,
     });
 
-    let mut visited = HashSet::new();
-    visited.insert(ep);
+    VISITED.with_borrow_mut(|visited| {
+        visited.begin(vectors.len());
+        visited.insert(ep);
+        search_layer_visit(
+            vectors, layer, query, ef, metric, visited, candidates, results,
+        )
+    })
+}
 
+/// Which nodes a layer search has already seen. A `HashSet` allocated per call
+/// was the hottest non-arithmetic cost in both build and query; this is one
+/// `u32` stamp per node, reused across calls: a node is visited iff its stamp
+/// equals the current epoch, so "clearing" is just bumping the epoch.
+#[derive(Default)]
+struct Visited {
+    stamps: Vec<u32>,
+    epoch: u32,
+}
+
+impl Visited {
+    fn begin(&mut self, n: usize) {
+        if self.stamps.len() < n {
+            self.stamps.resize(n, 0);
+        }
+        self.epoch = self.epoch.wrapping_add(1);
+        if self.epoch == 0 {
+            // Wrapped: stale stamps could collide with the new epoch.
+            self.stamps.iter_mut().for_each(|s| *s = 0);
+            self.epoch = 1;
+        }
+    }
+
+    /// Mark `id` visited; `true` if it was not already.
+    fn insert(&mut self, id: usize) -> bool {
+        let seen = self.stamps[id] == self.epoch;
+        self.stamps[id] = self.epoch;
+        !seen
+    }
+}
+
+thread_local! {
+    /// Per-thread scratch, so `search(&self)` stays shareable across threads.
+    static VISITED: std::cell::RefCell<Visited> = std::cell::RefCell::new(Visited::default());
+}
+
+#[allow(clippy::too_many_arguments)]
+fn search_layer_visit(
+    vectors: &[Vec<f32>],
+    layer: &[Vec<usize>],
+    query: &[f32],
+    ef: usize,
+    metric: DistanceMetric,
+    visited: &mut Visited,
+    mut candidates: BinaryHeap<Candidate>,
+    mut results: BinaryHeap<FarCandidate>,
+) -> Vec<Candidate> {
     while let Some(closest) = candidates.pop() {
         let furthest_dist = results.peek().map_or(f32::MAX, |f| f.distance);
         if closest.distance > furthest_dist && results.len() >= ef {
@@ -947,10 +1027,9 @@ fn search_layer(
         }
 
         for &neighbor in &layer[closest.id] {
-            if visited.contains(&neighbor) {
+            if !visited.insert(neighbor) {
                 continue;
             }
-            visited.insert(neighbor);
 
             let d = compute_distance(query, &vectors[neighbor], metric);
             let furthest_dist = results.peek().map_or(f32::MAX, |f| f.distance);
@@ -1236,6 +1315,7 @@ fn get_attr_string(attrs: &[(String, AttrValue)], name: &str) -> Result<String,
 #[cfg(test)]
 mod tests {
     use super::*;
+    use std::collections::HashSet;
 
     /// Tight, well-separated clusters — the shape real embeddings have, and
     /// the case plain closest-M neighbour selection fails on: each cluster