taladb 0.12.0

//! Exact vector storage and scoring. Persistent HNSW lives in `vector_graph`.
//!
//! Storage layout
//! ──────────────
//! Meta table  : `meta::vector_indexes`
//!   key       : `<collection>::<field>`   (UTF-8)
//!   value     : postcard-serialised `VectorDef`
//!
//! Vector table: `vec::<collection>::<field>`
//!   key       : ULID bytes (16 B, big-endian)
//!   value     : raw f32 LE bytes  (dimensions × 4 B)
//!
use std::collections::HashMap;
use std::sync::{Arc, Mutex};

use serde::{Deserialize, Serialize};

use crate::document::{Document, Value};

// ---------------------------------------------------------------------------
// Decoded-vector cache (flat search)
// ---------------------------------------------------------------------------
//
// The flat search path streams the `vec::` table and keeps only top-k scoring
// state. When the table fits the byte budget, this cache retains one contiguous
// decoded block so repeated queries are memory-bound rather than storage-bound.
//
// Each entry is tagged with the persisted vector-table revision read in the
// query's transaction. Embedding mutations bump it in the same write transaction;
// unrelated metadata writes preserve it. Flat and HNSW caches both reject stale
// revisions, including when another database handle performed the write.

/// Decoded vectors for one `collection::field`, tagged with its vector revision.
pub struct CachedVectors {
    pub generation: u64,
    pub vectors: Arc<VectorBlock>,
}

/// Contiguous decoded vectors; IDs and values avoid one allocation per vector.
#[derive(Default)]
pub struct VectorBlock {
    pub ids: Vec<ulid::Ulid>,
    pub values: Vec<f32>,
    pub dimensions: usize,
}
impl VectorBlock {
    pub fn bytes(&self) -> usize {
        self.ids.capacity() * std::mem::size_of::<ulid::Ulid>() + self.values.capacity() * 4
    }
}
/// Maximum retained decoded vector bytes per database. Oversized indexes stream.
pub const DEFAULT_VECTOR_CACHE_BYTES: usize = crate::search_cache::DEFAULT_SEARCH_CACHE_BYTES;

#[derive(Clone, Copy, PartialEq)]
pub(crate) struct Candidate(pub ulid::Ulid, pub f32);
impl Eq for Candidate {}
impl PartialOrd for Candidate {
    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
        Some(self.cmp(other))
    }
}
impl Ord for Candidate {
    fn cmp(&self, other: &Self) -> std::cmp::Ordering {
        self.1
            .total_cmp(&other.1)
            .then_with(|| other.0.cmp(&self.0))
    }
}

pub type VectorCacheMap = HashMap<String, CachedVectors>;
/// Shared across every `Collection` handle from the same `Database`, so a
/// write through one handle invalidates reads through another.
pub type SharedVectorCache = Arc<Mutex<VectorCacheMap>>;

pub fn new_shared_vector_cache() -> SharedVectorCache {
    Arc::new(Mutex::new(HashMap::new()))
}

// ---------------------------------------------------------------------------
// Constants
// ---------------------------------------------------------------------------

pub const META_VECTOR_TABLE: &str = "meta::vector_indexes";
/// Stores HNSW options: key = `"{collection}::{field}"`, value = postcard `HnswOptions`.
pub const META_HNSW_TABLE: &str = "meta::hnsw_indexes";
/// HNSW graph blob table: `hnsw::{collection}::{field}`, single key = `b"graph"`.
pub fn hnsw_table_name(collection: &str, field: &str) -> String {
    format!("hnsw::{collection}::{field}")
}

// ---------------------------------------------------------------------------
// VectorMetric
// ---------------------------------------------------------------------------

/// Similarity metric used when searching a vector index.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize, Default)]
pub enum VectorMetric {
    /// Cosine similarity — angle between vectors, range [-1, 1].
    /// Best for text embeddings where magnitude is not meaningful.
    #[default]
    Cosine,
    /// Raw dot product — magnitude-sensitive.
    Dot,
    /// Euclidean distance converted to similarity via `1 / (1 + dist)`.
    /// Range (0, 1]; identical vectors score 1.0.
    Euclidean,
}

// ---------------------------------------------------------------------------
// VectorDef — index metadata persisted in redb
// ---------------------------------------------------------------------------

#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct VectorDef {
    pub collection: String,
    pub field: String,
    /// Expected dimensionality; enforced on insert and search.
    pub dimensions: usize,
    pub metric: VectorMetric,
}

// ---------------------------------------------------------------------------
// HnswOptions — HNSW index configuration (stored in META_HNSW_TABLE)
// ---------------------------------------------------------------------------

/// Configuration for an HNSW vector index.
///
/// Only relevant when the `vector-hnsw` feature is enabled.  Storing these
/// options in a separate table keeps `VectorDef` backward-compatible: flat
/// indexes have no entry in `META_HNSW_TABLE`, HNSW indexes do.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct HnswOptions {
    /// Number of bi-directional links per node (connectivity).
    /// Supported range: 2 through 128.
    pub m: u32,
    /// Build-time quality parameter (ef during construction).
    /// Higher values produce a better graph at the cost of slower builds.
    /// Must be ≥ `m`. Default: 200.
    pub ef_construction: u32,
}

impl Default for HnswOptions {
    fn default() -> Self {
        Self {
            m: 32,
            ef_construction: 200,
        }
    }
}

// ---------------------------------------------------------------------------
// Table / key naming helpers
// ---------------------------------------------------------------------------

pub fn vec_meta_key(collection: &str, field: &str) -> String {
    format!("{collection}::{field}")
}

pub fn vec_table_name(collection: &str, field: &str) -> String {
    format!("vec::{collection}::{field}")
}

// ---------------------------------------------------------------------------
// Encoding / decoding
// ---------------------------------------------------------------------------

/// Encode a `&[f32]` as little-endian bytes.
pub fn encode_f32_vec(v: &[f32]) -> Vec<u8> {
    let mut buf = Vec::with_capacity(v.len() * 4);
    for f in v {
        buf.extend_from_slice(&f.to_le_bytes());
    }
    buf
}

/// Decode little-endian bytes back to `Vec<f32>`.
/// Returns `None` if `bytes.len()` is not a multiple of 4.
pub fn decode_f32_vec(bytes: &[u8]) -> Option<Vec<f32>> {
    // `as_chunks` splits and length-checks in one step: a non-empty remainder
    // *is* "not a multiple of 4", so the guard cannot drift from the split it
    // guards. The chunks come back as `[u8; 4]`, which `from_le_bytes` takes
    // directly — no per-element indexing to bounds-check.
    let (quads, rest) = bytes.as_chunks::<4>();
    if !rest.is_empty() {
        return None;
    }
    Some(quads.iter().copied().map(f32::from_le_bytes).collect())
}

/// Extract a float vector from a document field value.
/// Accepts `Value::Array` whose elements are all numeric (`Float` or `Int`).
/// Returns `None` if the value is not a numeric array.
pub fn value_to_f32_vec(v: &Value) -> Option<Vec<f32>> {
    match v {
        Value::Array(arr) => arr
            .iter()
            .map(|item| match item {
                Value::Float(f) => Some(*f as f32),
                Value::Int(n) => Some(*n as f32),
                _ => None,
            })
            .collect(),
        _ => None,
    }
}

// ---------------------------------------------------------------------------
// Similarity functions
// ---------------------------------------------------------------------------

pub fn cosine_similarity(a: &[f32], b: &[f32]) -> f32 {
    let dot: f32 = a.iter().zip(b.iter()).map(|(x, y)| x * y).sum();
    let norm_a: f32 = a.iter().map(|x| x * x).sum::<f32>().sqrt();
    let norm_b: f32 = b.iter().map(|x| x * x).sum::<f32>().sqrt();
    if norm_a == 0.0 || norm_b == 0.0 {
        0.0
    } else {
        dot / (norm_a * norm_b)
    }
}

pub fn dot_similarity(a: &[f32], b: &[f32]) -> f32 {
    // Lane-parallel for the same reason as `dot_and_norm_sq` — see `LANES`.
    let mut acc = [0.0f32; LANES];
    let (a_lanes, a_rest) = a.as_chunks::<LANES>();
    let (b_lanes, b_rest) = b.as_chunks::<LANES>();
    for (x, y) in a_lanes.iter().zip(b_lanes) {
        for i in 0..LANES {
            acc[i] += x[i] * y[i];
        }
    }
    let mut total: f32 = acc.iter().sum();
    for (x, y) in a_rest.iter().zip(b_rest) {
        total += x * y;
    }
    total
}

/// The norm half of [`dot_and_norm_sq`], on its own.
///
/// Lane-for-lane identical to the fused version — same accumulator count, same
/// summation order — so a norm computed once and cached is *bit-equal* to one
/// recomputed inside the fused loop. That matters more than it sounds: scores
/// order the result set, and a norm that differs in the last ulp can swap two
/// near-tied neighbours and move a document in or out of top-k.
#[inline]
pub fn norm_sq(v: &[f32]) -> f32 {
    let mut acc = [0.0f32; LANES];
    let (lanes, rest) = v.as_chunks::<LANES>();
    for s in lanes {
        for i in 0..LANES {
            acc[i] += s[i] * s[i];
        }
    }
    let mut total: f32 = acc.iter().sum();
    for s in rest {
        total += s * s;
    }
    total
}

/// [`score_with_query_norm`] with the *stored* vector's squared norm also known.
///
/// Cosine is the only metric that needs it, and it is a property of the stored
/// vector alone — invariant across every query it will ever be compared to. The
/// graph traversal scores the same cached nodes repeatedly, so recomputing it
/// per comparison was a full extra pass over the vector each time.
#[inline]
pub fn score_with_norms(
    metric: &VectorMetric,
    query: &[f32],
    query_norm: f32,
    stored: &[f32],
    stored_norm_sq: f32,
) -> f32 {
    match metric {
        VectorMetric::Cosine => {
            let norm_b = stored_norm_sq.sqrt();
            if query_norm == 0.0 || norm_b == 0.0 {
                0.0
            } else {
                dot_similarity(query, stored) / (query_norm * norm_b)
            }
        }
        _ => score_with_query_norm(metric, query, query_norm, stored),
    }
}

/// Converts Euclidean distance to a similarity score in `(0, 1]`.
/// Identical vectors → 1.0; the further apart, the closer to 0.
pub fn euclidean_similarity(a: &[f32], b: &[f32]) -> f32 {
    let mut acc = [0.0f32; LANES];
    let (a_lanes, a_rest) = a.as_chunks::<LANES>();
    let (b_lanes, b_rest) = b.as_chunks::<LANES>();
    for (x, y) in a_lanes.iter().zip(b_lanes) {
        for i in 0..LANES {
            let d = x[i] - y[i];
            acc[i] += d * d;
        }
    }
    let mut dist_sq: f32 = acc.iter().sum();
    for (x, y) in a_rest.iter().zip(b_rest) {
        let d = x - y;
        dist_sq += d * d;
    }
    1.0 / (1.0 + dist_sq.sqrt())
}

pub fn compute_similarity(metric: &VectorMetric, a: &[f32], b: &[f32]) -> f32 {
    match metric {
        VectorMetric::Cosine => cosine_similarity(a, b),
        VectorMetric::Dot => dot_similarity(a, b),
        VectorMetric::Euclidean => euclidean_similarity(a, b),
    }
}

/// L2 (Euclidean) norm of a vector. Hoist this out of a scan loop for cosine —
/// the query's norm is constant across every candidate, so recomputing it per
/// candidate (as `cosine_similarity` does) wastes one dot+sqrt per stored
/// vector.
#[inline]
pub fn l2_norm(a: &[f32]) -> f32 {
    a.iter().map(|x| x * x).sum::<f32>().sqrt()
}

/// [`compute_similarity`] with the query's L2 norm supplied by the caller.
///
/// This is the form a flat scan wants: `query_norm` is invariant across the
/// whole table, so computing it once at the top of the loop removes a dot
/// product and a square root per stored vector. `query_norm` is ignored for
/// metrics other than [`VectorMetric::Cosine`].
///
/// `query` and `stored` are expected to have the same length. They are walked
/// in lockstep and the shorter one ends the reduction, so a mismatch yields a
/// score computed over the common prefix rather than an error — callers that
/// can encounter untrusted or corrupt vectors should length-check first, as
/// `Collection::find_nearest` does.
#[inline]
pub fn score_with_query_norm(
    metric: &VectorMetric,
    query: &[f32],
    query_norm: f32,
    stored: &[f32],
) -> f32 {
    match metric {
        VectorMetric::Cosine => {
            let (dot, norm_sq) = dot_and_norm_sq(query, stored);
            let norm_b = norm_sq.sqrt();
            if query_norm == 0.0 || norm_b == 0.0 {
                0.0
            } else {
                dot / (query_norm * norm_b)
            }
        }
        VectorMetric::Dot => dot_similarity(query, stored),
        VectorMetric::Euclidean => euclidean_similarity(query, stored),
    }
}

/// Number of independent accumulators used by the scoring reductions.
///
/// Floating-point addition is not associative, so LLVM may not reorder a `+=`
/// chain into vector lanes: every element waits on the previous partial sum.
/// Splitting the reduction into `LANES` independent accumulators authorises
/// that reordering *in the source*, which is both portable — x86 AVX, aarch64
/// NEON, and the wasm128 SIMD every shipping browser has — and free of
/// `core::arch` intrinsics, `unsafe`, and nightly features.
///
/// Eight covers a 256-bit vector of `f32` and splits cleanly into two 128-bit
/// operations where that is all the target has. Embedding dimensions are
/// multiples of 8 in practice (384, 768, 1536), so the scalar remainder loop is
/// usually empty.
const LANES: usize = 8;

/// Dot product and the stored vector's squared L2 norm.
///
/// Cosine needs both. These were once accumulated in a single fused loop, to
/// read each element once and keep the pass memory-bound. That fusion cost far
/// more than the second walk saved: two `[f32; LANES]` accumulators live in one
/// loop body exceed what LLVM's SLP vectoriser will keep in registers, so it
/// gave up and emitted the whole reduction scalar. Measured with
/// `--emit=asm -C opt-level=3`, counting vector instructions in this function:
///
/// | target                | fused | split |
/// |-----------------------|-------|-------|
/// | wasm32 `+simd128`     |     0 |    74 |
/// | aarch64 (NEON)        |     2 |    26 |
/// | x86_64 (SSE2 baseline)|   114 |   132 |
///
/// x86_64 vectorised either way, which is why the regression went unseen: the
/// benchmark machine is x86. The browser and every Apple target — where this is
/// the flat scan's innermost loop — were running it scalar.
///
/// Both passes are zipped against `q_lanes`, not just the norm's own chunks, so
/// a length mismatch still ends the reduction at the shorter side exactly as the
/// fused loop did. Per-lane accumulation order is unchanged, so both totals are
/// bit-identical to the fused version — see the test at the bottom of this file.
#[inline]
fn dot_and_norm_sq(query: &[f32], stored: &[f32]) -> (f32, f32) {
    // `as_chunks` rather than `chunks_exact`: the chunks arrive as
    // `&[[f32; LANES]]`, so `q[i]` for `i in 0..LANES` is in bounds by type and
    // needs no bounds check at all. With `chunks_exact` each chunk is a
    // runtime-length `&[f32]` and the checks only vanish if the optimiser
    // happens to prove the range — which is a thin thing to rest a hot loop on.
    let (q_lanes, q_rest) = query.as_chunks::<LANES>();
    let (s_lanes, s_rest) = stored.as_chunks::<LANES>();

    let mut dot = [0.0f32; LANES];
    for (q, s) in q_lanes.iter().zip(s_lanes) {
        for i in 0..LANES {
            dot[i] += q[i] * s[i];
        }
    }

    // `_q` is unused, but zipping against it keeps this loop's element set
    // identical to the one above when the two slices differ in length.
    let mut norm_sq = [0.0f32; LANES];
    for (_q, s) in q_lanes.iter().zip(s_lanes) {
        for i in 0..LANES {
            norm_sq[i] += s[i] * s[i];
        }
    }

    let mut dot_total: f32 = dot.iter().sum();
    let mut norm_total: f32 = norm_sq.iter().sum();
    for (q, s) in q_rest.iter().zip(s_rest) {
        dot_total += q * s;
        norm_total += s * s;
    }
    (dot_total, norm_total)
}

// NOTE: a `score_from_bytes` variant that scored directly against the stored
// little-endian bytes used to live here. It had no callers: `find_nearest`
// decodes the whole vector table once into a generation-keyed cache
// (`Collection::find_nearest`, step 3b) and scores against `&[f32]` via
// `score_with_query_norm`, so nothing ever reaches the raw bytes in the hot
// loop. It was removed rather than wired up — its own doc comment claimed to be
// the flat scan's fast path, which sent at least one reviewer optimising a
// function that never ran.

// ---------------------------------------------------------------------------
// Search result
// ---------------------------------------------------------------------------

/// A single result returned by `Collection::find_nearest`.
#[derive(Debug, Clone)]
pub struct VectorSearchResult {
    /// The matched document (all fields, including `_id`).
    pub document: Document,
    /// Similarity score — higher is more similar.
    /// Range depends on the metric: cosine ∈ [-1,1], dot ∈ ℝ, euclidean ∈ (0,1].
    pub score: f32,
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn encode_decode_roundtrip() {
        let v = vec![0.1f32, -0.5, 1.0, 99.9];
        let decoded = decode_f32_vec(&encode_f32_vec(&v)).unwrap();
        for (a, b) in v.iter().zip(decoded.iter()) {
            assert!((a - b).abs() < 1e-6, "{a} != {b}");
        }
    }

    #[test]
    fn decode_rejects_odd_length() {
        assert!(decode_f32_vec(&[1, 2, 3]).is_none());
    }

    #[test]
    fn cosine_identical() {
        let a = vec![1.0f32, 2.0, 3.0];
        assert!((cosine_similarity(&a, &a) - 1.0).abs() < 1e-6);
    }

    #[test]
    fn cosine_orthogonal() {
        let a = vec![1.0f32, 0.0];
        let b = vec![0.0f32, 1.0];
        assert!(cosine_similarity(&a, &b).abs() < 1e-6);
    }

    #[test]
    fn euclidean_identical() {
        let a = vec![1.0f32, 2.0, 3.0];
        assert!((euclidean_similarity(&a, &a) - 1.0).abs() < 1e-6);
    }

    #[test]
    fn value_extraction_floats() {
        let v = Value::Array(vec![Value::Float(0.5), Value::Int(1), Value::Float(-0.25)]);
        let r = value_to_f32_vec(&v).unwrap();
        assert_eq!(r.len(), 3);
        assert!((r[0] - 0.5f32).abs() < 1e-6);
        assert_eq!(r[1], 1.0f32);
        assert!((r[2] - (-0.25f32)).abs() < 1e-6);
    }

    #[test]
    fn value_extraction_rejects_mixed() {
        let v = Value::Array(vec![Value::Float(0.5), Value::Str("x".into())]);
        assert!(value_to_f32_vec(&v).is_none());
    }

    #[test]
    fn value_extraction_rejects_non_array() {
        assert!(value_to_f32_vec(&Value::Str("vec".into())).is_none());
    }

    // -----------------------------------------------------------------------
    // Lane-parallel reductions
    //
    // The scoring loops split their sums across `LANES` independent
    // accumulators so the compiler can vectorise them. That changes the order
    // of floating-point additions, so these tests pin the results against a
    // straightforward scalar reference — including at lengths that are not a
    // multiple of `LANES`, where the remainder loop has to pick up the tail.
    // -----------------------------------------------------------------------

    fn scalar_dot(a: &[f32], b: &[f32]) -> f32 {
        a.iter().zip(b).map(|(x, y)| x * y).sum()
    }

    fn scalar_euclidean(a: &[f32], b: &[f32]) -> f32 {
        let d: f32 = a.iter().zip(b).map(|(x, y)| (x - y) * (x - y)).sum();
        1.0 / (1.0 + d.sqrt())
    }

    fn ramp(len: usize, offset: f32) -> Vec<f32> {
        (0..len)
            .map(|i| (i as f32).mul_add(0.125, offset))
            .collect()
    }

    #[test]
    fn lane_reductions_match_scalar_at_every_remainder_length() {
        // 0 and 1..=17 covers empty, sub-LANES, exactly LANES, and two full
        // chunks plus a tail.
        for len in std::iter::once(0).chain(1..=17) {
            let a = ramp(len, 0.5);
            let b = ramp(len, -0.25);

            let dot = dot_similarity(&a, &b);
            assert!(
                (dot - scalar_dot(&a, &b)).abs() < 1e-4,
                "dot mismatch at len {len}: {dot} vs {}",
                scalar_dot(&a, &b)
            );

            let euc = euclidean_similarity(&a, &b);
            assert!(
                (euc - scalar_euclidean(&a, &b)).abs() < 1e-6,
                "euclidean mismatch at len {len}"
            );

            // Cosine goes through the fused dot + norm pass.
            let norm = l2_norm(&a);
            let cos = score_with_query_norm(&VectorMetric::Cosine, &a, norm, &b);
            let expected = cosine_similarity(&a, &b);
            assert!(
                (cos - expected).abs() < 1e-5,
                "cosine mismatch at len {len}: {cos} vs {expected}"
            );
        }
    }

    #[test]
    fn identical_vectors_still_score_perfectly_after_lane_split() {
        // A reassociated sum must not drift far enough to break the invariant
        // callers rely on: a vector is maximally similar to itself.
        for len in [4usize, 8, 12, 384] {
            let v = ramp(len, 1.0);
            let cos = score_with_query_norm(&VectorMetric::Cosine, &v, l2_norm(&v), &v);
            assert!(
                (cos - 1.0).abs() < 1e-5,
                "cosine self-similarity at len {len}: {cos}"
            );
            let euc = euclidean_similarity(&v, &v);
            assert!(
                (euc - 1.0).abs() < 1e-6,
                "euclidean self-similarity at len {len}"
            );
        }
    }
    /// `dot_and_norm_sq` was refactored from `chunks_exact(LANES)` to
    /// `as_chunks::<LANES>()` to clear `clippy::chunks_exact_to_as_chunks`. It is
    /// a hot float reduction whose *lane structure is the point*, so the bar is
    /// bit-identical output, not "close enough" — a reordered accumulation would
    /// silently shift every cosine score in the database.
    ///
    /// This lives in-module because `dot_and_norm_sq` is private; the public
    /// `dot_similarity` / `euclidean_similarity` / `decode_f32_vec` were checked
    /// the same way.
    #[test]
    fn dot_and_norm_sq_matches_the_chunks_exact_original_bit_for_bit() {
        // The whole point of `original` is to be the pre-refactor code, so the
        // lint that prompted the refactor must not be applied to it. Rewriting
        // it to `as_chunks` would make it a copy of the new implementation and
        // the test would compare something to itself.
        #[allow(clippy::chunks_exact_to_as_chunks)]
        fn original(query: &[f32], stored: &[f32]) -> (f32, f32) {
            let mut dot = [0.0f32; LANES];
            let mut norm_sq = [0.0f32; LANES];
            let mut q_chunks = query.chunks_exact(LANES);
            let mut s_chunks = stored.chunks_exact(LANES);
            for (q, s) in q_chunks.by_ref().zip(s_chunks.by_ref()) {
                for i in 0..LANES {
                    dot[i] += q[i] * s[i];
                    norm_sq[i] += s[i] * s[i];
                }
            }
            let mut dot_total: f32 = dot.iter().sum();
            let mut norm_total: f32 = norm_sq.iter().sum();
            for (q, s) in q_chunks.remainder().iter().zip(s_chunks.remainder()) {
                dot_total += q * s;
                norm_total += s * s;
            }
            (dot_total, norm_total)
        }

        // Deterministic, and spread across magnitudes so summation order matters.
        let mut seed = 0xC0FF_EE12_u64;
        let mut next = || {
            seed = seed
                .wrapping_mul(6364136223846793005)
                .wrapping_add(1442695040888963407);
            let bits = (seed >> 33) as u32;
            (f32::from_bits((bits & 0x007F_FFFF) | 0x3F80_0000) - 1.5) * 1e3
        };

        // Lengths either side of LANES, and real embedding sizes, so both the
        // lane loop and the remainder tail are exercised.
        for len in [0usize, 1, 7, 8, 9, 15, 16, 17, 31, 384, 385, 768, 1536] {
            let q: Vec<f32> = (0..len).map(|_| next()).collect();
            let s: Vec<f32> = (0..len).map(|_| next()).collect();
            let (new_dot, new_norm) = dot_and_norm_sq(&q, &s);
            let (old_dot, old_norm) = original(&q, &s);
            assert_eq!(
                new_dot.to_bits(),
                old_dot.to_bits(),
                "dot diverged at len={len}"
            );
            assert_eq!(
                new_norm.to_bits(),
                old_norm.to_bits(),
                "norm_sq diverged at len={len}"
            );
        }
    }

    /// Splitting the fused pass into two loops introduced a way to get this
    /// wrong that the equal-length test above cannot see: if the norm loop
    /// walked the stored vector's own chunks instead of zipping against the
    /// query's, a shorter query would leave the norm accumulating over stored
    /// elements the dot product never reached.
    ///
    /// The reference is the fused original, not a hand-derived "common prefix":
    /// when the two slices land in different numbers of lane chunks the
    /// remainder zip pairs elements at different offsets (at `ql=7, sl=9` it
    /// multiplies `q[0]` by `s[8]`). That is odd, but it predates this split and
    /// is only reachable on a length mismatch, which `find_nearest` rejects
    /// before scoring. What must not change is that the split agrees with the
    /// fused loop bit for bit.
    #[test]
    fn mismatched_lengths_match_the_fused_original_bit_for_bit() {
        #[allow(clippy::chunks_exact_to_as_chunks)]
        fn fused(query: &[f32], stored: &[f32]) -> (f32, f32) {
            let mut dot = [0.0f32; LANES];
            let mut norm_sq = [0.0f32; LANES];
            let mut q_chunks = query.chunks_exact(LANES);
            let mut s_chunks = stored.chunks_exact(LANES);
            for (q, s) in q_chunks.by_ref().zip(s_chunks.by_ref()) {
                for i in 0..LANES {
                    dot[i] += q[i] * s[i];
                    norm_sq[i] += s[i] * s[i];
                }
            }
            let mut dot_total: f32 = dot.iter().sum();
            let mut norm_total: f32 = norm_sq.iter().sum();
            for (q, s) in q_chunks.remainder().iter().zip(s_chunks.remainder()) {
                dot_total += q * s;
                norm_total += s * s;
            }
            (dot_total, norm_total)
        }

        let mut seed = 0x5DEE_CE66_D125_u64;
        let mut next = move || {
            seed = seed
                .wrapping_mul(6364136223846793005)
                .wrapping_add(1442695040888963407);
            let bits = (seed >> 33) as u32;
            (f32::from_bits((bits & 0x007F_FFFF) | 0x3F80_0000) - 1.5) * 1e3
        };

        // Both orderings, mismatches inside a lane chunk and across a chunk
        // boundary, and the empty-side cases.
        for (ql, sl) in [
            (8usize, 16usize),
            (16, 8),
            (7, 9),
            (9, 7),
            (384, 385),
            (385, 384),
            (0, 8),
            (8, 0),
        ] {
            let q: Vec<f32> = (0..ql).map(|_| next()).collect();
            let s: Vec<f32> = (0..sl).map(|_| next()).collect();
            let (dot, norm) = dot_and_norm_sq(&q, &s);
            let (fd, fnm) = fused(&q, &s);
            assert_eq!(dot.to_bits(), fd.to_bits(), "dot at ql={ql} sl={sl}");
            assert_eq!(norm.to_bits(), fnm.to_bits(), "norm at ql={ql} sl={sl}");
        }
    }
}