diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index c874e1b361..30ee4ae606 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -354,6 +354,8 @@ jobs:
               --package diskann-wide \
               --package diskann-vector \
               --package diskann-quantization \
+              --package diskann \
+              --features pipnn \
               -- --skip compile_tests \
                  --skip pivots::tests::run_test_happy_path
 
@@ -416,6 +418,8 @@ jobs:
             --package diskann-wide \
             --package diskann-vector \
             --package diskann-quantization \
+            --package diskann \
+            --features pipnn \
             -- --skip compile_tests
 
   test-workspace:
@@ -459,6 +463,7 @@ jobs:
         os:
           - windows-latest
           - ubuntu-latest
+          - ubuntu-24.04-arm
 
     steps:
       - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1
diff --git a/Cargo.lock b/Cargo.lock
index 2588eee92a..40240e77c5 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -438,6 +438,7 @@ dependencies = [
  "anyhow",
  "bytemuck",
  "dashmap",
+ "diskann-linalg",
  "diskann-utils",
  "diskann-vector",
  "diskann-wide",
@@ -448,6 +449,7 @@ dependencies = [
  "pin-project",
  "rand",
  "relative-path 2.0.1",
+ "rstest",
  "serde",
  "serde_json",
  "thiserror 2.0.17",
diff --git a/diskann-linalg/src/faer.rs b/diskann-linalg/src/faer.rs
index 700396de4e..f4f498a13b 100644
--- a/diskann-linalg/src/faer.rs
+++ b/diskann-linalg/src/faer.rs
@@ -53,6 +53,33 @@ pub(super) fn sgemm_impl(
     faer::linalg::matmul::matmul(c, beta, a, b, alpha, Par::Seq)
 }
 
+/// Compute the lower triangle of `A * Aᵀ` with Faer.
+///
+/// Leaf selection reads each symmetric pair once. It does not read the upper
+/// triangle. `BlockStructure::TriangularLower` prevents writes to that triangle.
+///
+/// `sgemm_aat_lower` checks both slice lengths and both size products. Therefore,
+/// the Faer matrix views stay inside their backing slices.
+pub(super) fn sgemm_aat_lower_impl(m: usize, k: usize, a: &[f32], c: &mut [f32]) {
+    use faer::linalg::matmul::triangular::{matmul, BlockStructure};
+
+    let a = faer::mat::MatRef::from_row_major_slice(a, m, k);
+    let at = a.transpose();
+    let c = faer::mat::MatMut::from_row_major_slice_mut(c, m, m);
+
+    matmul(
+        c,
+        BlockStructure::TriangularLower,
+        faer::Accum::Replace,
+        a,
+        BlockStructure::Rectangular,
+        at,
+        BlockStructure::Rectangular,
+        1.0,
+        Par::Seq,
+    );
+}
+
 /// See the documentation for `svd_into`.
 ///
 /// The implementation may assume the the specified invariants hold for the sizes of the
diff --git a/diskann-linalg/src/lib.rs b/diskann-linalg/src/lib.rs
index 7ee59d6b60..778fb94cf8 100644
--- a/diskann-linalg/src/lib.rs
+++ b/diskann-linalg/src/lib.rs
@@ -82,6 +82,30 @@ impl fmt::Display for SgemmError {
 
 impl std::error::Error for SgemmError {}
 
+fn check_matrix(
+    matrix_name: MatrixName,
+    actual_len: usize,
+    rows: usize,
+    cols: usize,
+) -> Result<(), SgemmError> {
+    let expected_len = rows
+        .checked_mul(cols)
+        .ok_or(SgemmError::DimensionOverflow {
+            matrix_name,
+            rows,
+            cols,
+        })?;
+    if actual_len != expected_len {
+        return Err(SgemmError::InvalidMatrixDimensions {
+            matrix_name,
+            expected_rows: rows,
+            expected_cols: cols,
+            actual_len,
+        });
+    }
+    Ok(())
+}
+
 // Make the reference implementation available for internal testing.
 #[cfg(test)]
 mod reference;
@@ -156,57 +180,33 @@ pub fn sgemm(
     beta: Option<f32>,
     c: &mut [f32],
 ) -> Result<(), SgemmError> {
-    // Check size requirements with overflow protection.
-    let expected_a_len = m.checked_mul(k).ok_or(SgemmError::DimensionOverflow {
-        matrix_name: MatrixName::A,
-        rows: m,
-        cols: k,
-    })?;
-
-    if a.len() != expected_a_len {
-        return Err(SgemmError::InvalidMatrixDimensions {
-            matrix_name: MatrixName::A,
-            expected_rows: m,
-            expected_cols: k,
-            actual_len: a.len(),
-        });
-    }
-
-    let expected_b_len = k.checked_mul(n).ok_or(SgemmError::DimensionOverflow {
-        matrix_name: MatrixName::B,
-        rows: k,
-        cols: n,
-    })?;
-
-    if b.len() != expected_b_len {
-        return Err(SgemmError::InvalidMatrixDimensions {
-            matrix_name: MatrixName::B,
-            expected_rows: k,
-            expected_cols: n,
-            actual_len: b.len(),
-        });
-    }
-
-    let expected_c_len = m.checked_mul(n).ok_or(SgemmError::DimensionOverflow {
-        matrix_name: MatrixName::C,
-        rows: m,
-        cols: n,
-    })?;
-
-    if c.len() != expected_c_len {
-        return Err(SgemmError::InvalidMatrixDimensions {
-            matrix_name: MatrixName::C,
-            expected_rows: m,
-            expected_cols: n,
-            actual_len: c.len(),
-        });
-    }
+    check_matrix(MatrixName::A, a.len(), m, k)?;
+    check_matrix(MatrixName::B, b.len(), k, n)?;
+    check_matrix(MatrixName::C, c.len(), m, n)?;
 
     // Invoke the actual implementation.
     sgemm_impl(atranspose, btranspose, m, n, k, alpha, a, b, beta, c);
     Ok(())
 }
 
+/// Compute the lower triangle of $C = A A^\mathsf{T}$.
+///
+/// `A` is a dense row-major $m \times k$ matrix. The function overwrites the
+/// lower triangle of `C`, including its diagonal. It does not change the upper
+/// triangle.
+///
+/// # Errors
+///
+/// Returns an error if a size product overflows. It also returns an error if a
+/// slice length does not match its declared matrix shape.
+pub fn sgemm_aat_lower(m: usize, k: usize, a: &[f32], c: &mut [f32]) -> Result<(), SgemmError> {
+    check_matrix(MatrixName::A, a.len(), m, k)?;
+    check_matrix(MatrixName::C, c.len(), m, m)?;
+
+    faer::sgemm_aat_lower_impl(m, k, a, c);
+    Ok(())
+}
+
 /// Compute the SVD of the provided matrix implicit row-major matrix `data`.
 ///
 /// * `m`: The number of rows in `a`.
@@ -690,3 +690,115 @@ mod tests {
         }
     }
 }
+#[cfg(test)]
+#[allow(
+    clippy::expect_used,
+    clippy::unwrap_used,
+    reason = "deterministic test fixture construction must abort on invalid setup"
+)]
+mod sgemm_aat_lower_tests {
+    use super::{sgemm_aat_lower, MatrixName, SgemmError};
+
+    #[test]
+    fn computes_lower_triangle_and_preserves_upper_triangle() {
+        #[rustfmt::skip]
+        let a = [
+            1.0, 2.0,
+            3.0, 4.0,
+            5.0, 6.0,
+        ];
+        let untouched = -123.0;
+        let mut c = [untouched; 9];
+
+        sgemm_aat_lower(3, 2, &a, &mut c).unwrap();
+
+        #[rustfmt::skip]
+        assert_eq!(c, [
+             5.0, untouched, untouched,
+            11.0,     25.0, untouched,
+            17.0,     39.0,      61.0,
+        ]);
+    }
+
+    #[test]
+    fn accepts_a_matrix_with_no_rows() {
+        sgemm_aat_lower(0, 3, &[], &mut []).unwrap();
+    }
+
+    #[test]
+    fn zero_inner_dimension_zeros_only_the_lower_triangle() {
+        let untouched = -123.0;
+        let mut c = [untouched; 9];
+
+        sgemm_aat_lower(3, 0, &[], &mut c).unwrap();
+
+        #[rustfmt::skip]
+        assert_eq!(c, [
+            0.0, untouched, untouched,
+            0.0,       0.0, untouched,
+            0.0,       0.0,       0.0,
+        ]);
+    }
+
+    #[test]
+    fn rejects_invalid_input_dimensions() {
+        let mut c = [0.0; 4];
+
+        let error = sgemm_aat_lower(2, 2, &[0.0; 3], &mut c).unwrap_err();
+
+        assert_eq!(
+            error,
+            SgemmError::InvalidMatrixDimensions {
+                matrix_name: MatrixName::A,
+                expected_rows: 2,
+                expected_cols: 2,
+                actual_len: 3,
+            }
+        );
+    }
+
+    #[test]
+    fn rejects_invalid_output_dimensions() {
+        let mut c = [0.0; 3];
+
+        let error = sgemm_aat_lower(2, 2, &[0.0; 4], &mut c).unwrap_err();
+
+        assert_eq!(
+            error,
+            SgemmError::InvalidMatrixDimensions {
+                matrix_name: MatrixName::C,
+                expected_rows: 2,
+                expected_cols: 2,
+                actual_len: 3,
+            }
+        );
+    }
+
+    #[test]
+    fn rejects_input_size_overflow() {
+        let error = sgemm_aat_lower(usize::MAX, 2, &[], &mut []).unwrap_err();
+
+        assert_eq!(
+            error,
+            SgemmError::DimensionOverflow {
+                matrix_name: MatrixName::A,
+                rows: usize::MAX,
+                cols: 2,
+            }
+        );
+    }
+
+    #[test]
+    fn rejects_output_size_overflow() {
+        let error = sgemm_aat_lower(usize::MAX, 0, &[], &mut []).unwrap_err();
+
+        assert_eq!(
+            error,
+            SgemmError::DimensionOverflow {
+                matrix_name: MatrixName::C,
+                rows: usize::MAX,
+                cols: usize::MAX,
+            }
+        );
+    }
+}
diff --git a/diskann-utils/src/views.rs b/diskann-utils/src/views.rs
index a9352918c9..e38b5b6c58 100644
--- a/diskann-utils/src/views.rs
+++ b/diskann-utils/src/views.rs
@@ -211,7 +211,7 @@ where
     /// The length of the base must be equal to `nrows * ncols`.
     pub fn try_from(data: T, nrows: usize, ncols: usize) -> Result<Self, TryFromError<T>> {
         let len = data.as_slice().len();
-        if len != nrows * ncols {
+        if nrows.checked_mul(ncols) != Some(len) {
             Err(TryFromError { data, nrows, ncols })
         } else {
             Ok(Self { data, nrows, ncols })
@@ -1053,6 +1053,8 @@ mod tests {
             m.unwrap_err().to_string(),
             "tried to construct a matrix view with 5 rows and 4 columns over a slice of length 12"
         );
+
+        assert!(MatrixView::try_from(&[] as &[usize], usize::MAX / 2 + 1, 2).is_err());
     }
 
     #[test]
diff --git a/diskann-wide/src/arch/aarch64/f32x2_.rs b/diskann-wide/src/arch/aarch64/f32x2_.rs
index 318227ca5d..f0b7431a87 100644
--- a/diskann-wide/src/arch/aarch64/f32x2_.rs
+++ b/diskann-wide/src/arch/aarch64/f32x2_.rs
@@ -31,6 +31,7 @@ macros::aarch64_define_loadstore!(f32x2, vld1_f32, internal::load_first::f32x2,
 helpers::unsafe_map_binary_op!(f32x2, std::ops::Add, add, vadd_f32, "neon");
 helpers::unsafe_map_binary_op!(f32x2, std::ops::Sub, sub, vsub_f32, "neon");
 helpers::unsafe_map_binary_op!(f32x2, std::ops::Mul, mul, vmul_f32, "neon");
+helpers::unsafe_map_binary_op!(f32x2, std::ops::Div, div, vdiv_f32, "neon");
 macros::aarch64_define_fma!(f32x2, vfma_f32);
 
 macros::aarch64_define_cmp!(
@@ -90,6 +91,7 @@ mod tests {
     test_utils::ops::test_add!(f32x2, 0xcd7a8fea9a3fb727, test_neon());
     test_utils::ops::test_sub!(f32x2, 0x3f6562c94c923238, test_neon());
     test_utils::ops::test_mul!(f32x2, 0x07e48666c0fc564c, test_neon());
+    test_utils::ops::test_div!(f32x2, 0xa0352efeb9bc5ca5, test_neon());
     test_utils::ops::test_fma!(f32x2, 0xcfde9d031302cf2c, test_neon());
 
     test_utils::ops::test_cmp!(f32x2, 0xc4f468b224622326, test_neon());
diff --git a/diskann-wide/src/arch/aarch64/f32x4_.rs b/diskann-wide/src/arch/aarch64/f32x4_.rs
index 82cf391076..83779dacf0 100644
--- a/diskann-wide/src/arch/aarch64/f32x4_.rs
+++ b/diskann-wide/src/arch/aarch64/f32x4_.rs
@@ -32,6 +32,7 @@ macros::aarch64_splitjoin!(f32x4, f32x2, vget_low_f32, vget_high_f32, vcombine_f
 helpers::unsafe_map_binary_op!(f32x4, std::ops::Add, add, vaddq_f32, "neon");
 helpers::unsafe_map_binary_op!(f32x4, std::ops::Sub, sub, vsubq_f32, "neon");
 helpers::unsafe_map_binary_op!(f32x4, std::ops::Mul, mul, vmulq_f32, "neon");
+helpers::unsafe_map_binary_op!(f32x4, std::ops::Div, div, vdivq_f32, "neon");
 helpers::unsafe_map_unary_op!(f32x4, SIMDAbs, abs_simd, vabsq_f32, "neon");
 macros::aarch64_define_fma!(f32x4, vfmaq_f32);
 
@@ -187,6 +188,7 @@ mod tests {
     test_utils::ops::test_add!(f32x4, 0xcd7a8fea9a3fb727, test_neon());
     test_utils::ops::test_sub!(f32x4, 0x3f6562c94c923238, test_neon());
     test_utils::ops::test_mul!(f32x4, 0x07e48666c0fc564c, test_neon());
+    test_utils::ops::test_div!(f32x4, 0xa0352efeb9bc5ca5, test_neon());
     test_utils::ops::test_fma!(f32x4, 0xcfde9d031302cf2c, test_neon());
     test_utils::ops::test_abs!(f32x4, 0xb8f702ba85375041, test_neon());
     test_utils::ops::test_minmax!(f32x4, 0x6d7fc8ed6d852187, test_neon());
diff --git a/diskann-wide/src/arch/x86_64/v3/f32x16_.rs b/diskann-wide/src/arch/x86_64/v3/f32x16_.rs
index 836193a452..b93c861e7a 100644
--- a/diskann-wide/src/arch/x86_64/v3/f32x16_.rs
+++ b/diskann-wide/src/arch/x86_64/v3/f32x16_.rs
@@ -54,6 +54,7 @@ mod test_x86_f32 {
     test_utils::ops::test_add!(f32x16, 0xa8989b97ca888d11, V3::new_checked_uncached());
     test_utils::ops::test_sub!(f32x16, 0xb2554fc13fdc1182, V3::new_checked_uncached());
     test_utils::ops::test_mul!(f32x16, 0x23becaa968b0cd71, V3::new_checked_uncached());
+    test_utils::ops::test_div!(f32x16, 0x6fd16af08fa1f498, V3::new_checked_uncached());
     test_utils::ops::test_fma!(f32x16, 0x32a814070a93df4e, V3::new_checked_uncached());
     test_utils::ops::test_minmax!(f32x16, 0x6d7fc8ed6d852187, V3::new_checked_uncached());
     test_utils::ops::test_abs!(f32x16, 0x6799e60873a2efe2, V3::new_checked_uncached());
diff --git a/diskann-wide/src/arch/x86_64/v3/f32x4_.rs b/diskann-wide/src/arch/x86_64/v3/f32x4_.rs
index 60ffa4477c..45cc64a4a2 100644
--- a/diskann-wide/src/arch/x86_64/v3/f32x4_.rs
+++ b/diskann-wide/src/arch/x86_64/v3/f32x4_.rs
@@ -31,6 +31,7 @@ macros::x86_define_default!(f32x4, _mm_setzero_ps, "sse");
 helpers::unsafe_map_binary_op!(f32x4, std::ops::Add, add, _mm_add_ps, "sse");
 helpers::unsafe_map_binary_op!(f32x4, std::ops::Sub, sub, _mm_sub_ps, "sse");
 helpers::unsafe_map_binary_op!(f32x4, std::ops::Mul, mul, _mm_mul_ps, "sse");
+helpers::unsafe_map_binary_op!(f32x4, std::ops::Div, div, _mm_div_ps, "sse");
 
 impl f32x4 {
     #[inline(always)]
@@ -253,6 +254,7 @@ mod test_x86_f32 {
     test_utils::ops::test_add!(f32x4, 0xcd7a8fea9a3fb727, V3::new_checked_uncached());
     test_utils::ops::test_sub!(f32x4, 0x3f6562c94c923238, V3::new_checked_uncached());
     test_utils::ops::test_mul!(f32x4, 0x07e48666c0fc564c, V3::new_checked_uncached());
+    test_utils::ops::test_div!(f32x4, 0xa0352efeb9bc5ca5, V3::new_checked_uncached());
     test_utils::ops::test_fma!(f32x4, 0xcfde9d031302cf2c, V3::new_checked_uncached());
     test_utils::ops::test_minmax!(f32x4, 0x6d7fc8ed6d852187, V3::new_checked_uncached());
     test_utils::ops::test_abs!(f32x4, 0x8e6d9944c9c43a74, V3::new_checked_uncached());
diff --git a/diskann-wide/src/arch/x86_64/v3/f32x8_.rs b/diskann-wide/src/arch/x86_64/v3/f32x8_.rs
index 054b249e8f..48ecea7aee 100644
--- a/diskann-wide/src/arch/x86_64/v3/f32x8_.rs
+++ b/diskann-wide/src/arch/x86_64/v3/f32x8_.rs
@@ -33,6 +33,7 @@ macros::x86_splitjoin!(f32x8, f32x4, _mm256_extractf128_ps, _mm256_set_m128, "av
 helpers::unsafe_map_binary_op!(f32x8, std::ops::Add, add, _mm256_add_ps, "avx");
 helpers::unsafe_map_binary_op!(f32x8, std::ops::Sub, sub, _mm256_sub_ps, "avx");
 helpers::unsafe_map_binary_op!(f32x8, std::ops::Mul, mul, _mm256_mul_ps, "avx");
+helpers::unsafe_map_binary_op!(f32x8, std::ops::Div, div, _mm256_div_ps, "avx");
 
 impl f32x8 {
     #[inline(always)]
@@ -266,6 +267,7 @@ mod test_x86_f32 {
     test_utils::ops::test_add!(f32x8, 0x3824379d4a43a416, V3::new_checked_uncached());
     test_utils::ops::test_sub!(f32x8, 0x548fc74c07ba425d, V3::new_checked_uncached());
     test_utils::ops::test_mul!(f32x8, 0x6d340672ff91b256, V3::new_checked_uncached());
+    test_utils::ops::test_div!(f32x8, 0x776f54898c62dd0b, V3::new_checked_uncached());
     test_utils::ops::test_fma!(f32x8, 0x5f566d8968d4d201, V3::new_checked_uncached());
     test_utils::ops::test_minmax!(f32x8, 0x6d7fc8ed6d852187, V3::new_checked_uncached());
     test_utils::ops::test_abs!(f32x8, 0x2a4a9651d8ebe912, V3::new_checked_uncached());
diff --git a/diskann-wide/src/arch/x86_64/v4/f32x16_.rs b/diskann-wide/src/arch/x86_64/v4/f32x16_.rs
index d38465f906..ff59119992 100644
--- a/diskann-wide/src/arch/x86_64/v4/f32x16_.rs
+++ b/diskann-wide/src/arch/x86_64/v4/f32x16_.rs
@@ -57,6 +57,7 @@ impl crate::SplitJoin for f32x16 {
 helpers::unsafe_map_binary_op!(f32x16, std::ops::Add, add, _mm512_add_ps, "avx512f");
 helpers::unsafe_map_binary_op!(f32x16, std::ops::Sub, sub, _mm512_sub_ps, "avx512f");
 helpers::unsafe_map_binary_op!(f32x16, std::ops::Mul, mul, _mm512_mul_ps, "avx512f");
+helpers::unsafe_map_binary_op!(f32x16, std::ops::Div, div, _mm512_div_ps, "avx512f");
 
 impl f32x16 {
     #[inline(always)]
@@ -240,6 +241,7 @@ mod test_x86_f32 {
     test_utils::ops::test_add!(f32x16, 0xa8989b97ca888d11, V4::new_checked_uncached());
     test_utils::ops::test_sub!(f32x16, 0xb2554fc13fdc1182, V4::new_checked_uncached());
     test_utils::ops::test_mul!(f32x16, 0x23becaa968b0cd71, V4::new_checked_uncached());
+    test_utils::ops::test_div!(f32x16, 0x6fd16af08fa1f498, V4::new_checked_uncached());
     test_utils::ops::test_fma!(f32x16, 0x32a814070a93df4e, V4::new_checked_uncached());
     test_utils::ops::test_minmax!(f32x16, 0x6d7fc8ed6d852187, V4::new_checked_uncached());
     test_utils::ops::test_abs!(f32x16, 0x6799e60873a2efe2, V4::new_checked_uncached());
diff --git a/diskann-wide/src/arch/x86_64/v4/f32x4_.rs b/diskann-wide/src/arch/x86_64/v4/f32x4_.rs
index 328dba4d26..7028b2fdcf 100644
--- a/diskann-wide/src/arch/x86_64/v4/f32x4_.rs
+++ b/diskann-wide/src/arch/x86_64/v4/f32x4_.rs
@@ -33,6 +33,7 @@ macros::x86_retarget!(f32x4 => v3::f32x4);
 helpers::unsafe_map_binary_op!(f32x4, std::ops::Add, add, _mm_add_ps, "sse");
 helpers::unsafe_map_binary_op!(f32x4, std::ops::Sub, sub, _mm_sub_ps, "sse");
 helpers::unsafe_map_binary_op!(f32x4, std::ops::Mul, mul, _mm_mul_ps, "sse");
+helpers::unsafe_map_binary_op!(f32x4, std::ops::Div, div, _mm_div_ps, "sse");
 
 impl f32x4 {
     #[inline(always)]
@@ -210,6 +211,7 @@ mod test_x86_f32 {
     test_utils::ops::test_add!(f32x4, 0xcd7a8fea9a3fb727, V4::new_checked_uncached());
     test_utils::ops::test_sub!(f32x4, 0x3f6562c94c923238, V4::new_checked_uncached());
     test_utils::ops::test_mul!(f32x4, 0x07e48666c0fc564c, V4::new_checked_uncached());
+    test_utils::ops::test_div!(f32x4, 0xa0352efeb9bc5ca5, V4::new_checked_uncached());
     test_utils::ops::test_fma!(f32x4, 0xcfde9d031302cf2c, V4::new_checked_uncached());
     test_utils::ops::test_minmax!(f32x4, 0x6d7fc8ed6d852187, V4::new_checked_uncached());
     test_utils::ops::test_abs!(f32x4, 0x8e6d9944c9c43a74, V4::new_checked_uncached());
diff --git a/diskann-wide/src/arch/x86_64/v4/f32x8_.rs b/diskann-wide/src/arch/x86_64/v4/f32x8_.rs
index 3158ffc1dd..d38de49de8 100644
--- a/diskann-wide/src/arch/x86_64/v4/f32x8_.rs
+++ b/diskann-wide/src/arch/x86_64/v4/f32x8_.rs
@@ -36,6 +36,7 @@ macros::x86_retarget!(f32x8 => v3::f32x8);
 helpers::unsafe_map_binary_op!(f32x8, std::ops::Add, add, _mm256_add_ps, "avx");
 helpers::unsafe_map_binary_op!(f32x8, std::ops::Sub, sub, _mm256_sub_ps, "avx");
 helpers::unsafe_map_binary_op!(f32x8, std::ops::Mul, mul, _mm256_mul_ps, "avx");
+helpers::unsafe_map_binary_op!(f32x8, std::ops::Div, div, _mm256_div_ps, "avx");
 
 impl f32x8 {
     #[inline(always)]
@@ -206,6 +207,7 @@ mod test_x86_f32 {
     test_utils::ops::test_add!(f32x8, 0x3824379d4a43a416, V4::new_checked_uncached());
     test_utils::ops::test_sub!(f32x8, 0x548fc74c07ba425d, V4::new_checked_uncached());
     test_utils::ops::test_mul!(f32x8, 0x6d340672ff91b256, V4::new_checked_uncached());
+    test_utils::ops::test_div!(f32x8, 0x776f54898c62dd0b, V4::new_checked_uncached());
     test_utils::ops::test_fma!(f32x8, 0x5f566d8968d4d201, V4::new_checked_uncached());
     test_utils::ops::test_minmax!(f32x8, 0x6d7fc8ed6d852187, V4::new_checked_uncached());
     test_utils::ops::test_abs!(f32x8, 0x2a4a9651d8ebe912, V4::new_checked_uncached());
diff --git a/diskann-wide/src/doubled.rs b/diskann-wide/src/doubled.rs
index 30d08e6cb1..d6adcb7b13 100644
--- a/diskann-wide/src/doubled.rs
+++ b/diskann-wide/src/doubled.rs
@@ -205,6 +205,15 @@ impl<T: std::ops::Mul<Output = T>> std::ops::Mul for Doubled<T> {
     }
 }
 
+impl<T: std::ops::Div<Output = T>> std::ops::Div for Doubled<T> {
+    type Output = Self;
+
+    #[inline(always)]
+    fn div(self, rhs: Self) -> Self {
+        Self(self.0 / rhs.0, self.1 / rhs.1)
+    }
+}
+
 impl<T: std::ops::BitAnd<Output = T>> std::ops::BitAnd for Doubled<T> {
     type Output = Self;
     #[inline(always)]
diff --git a/diskann-wide/src/emulated.rs b/diskann-wide/src/emulated.rs
index 507dc1288f..1e8d707769 100644
--- a/diskann-wide/src/emulated.rs
+++ b/diskann-wide/src/emulated.rs
@@ -199,6 +199,15 @@ where
     }
 }
 
+impl<const N: usize, A> std::ops::Div for Emulated<f32, N, A> {
+    type Output = Self;
+
+    #[inline(always)]
+    fn div(self, rhs: Self) -> Self {
+        Self::from_arch_fn(self.1, |i| self.0[i] / rhs.0[i])
+    }
+}
+
 /// MulAdd
 impl<T, const N: usize, A> SIMDMulAdd for Emulated<T, N, A>
 where
@@ -902,6 +911,10 @@ mod test_emulated {
     test_emulated!(f32, 4);
     test_emulated!(f32, 8);
     test_emulated!(f32, 16);
+    test_utils::ops::test_div!(Emulated<f32, 1>, 0x32f0d2991be50f13, SC);
+    test_utils::ops::test_div!(Emulated<f32, 4>, 0xf65f08475f5e30c9, SC);
+    test_utils::ops::test_div!(Emulated<f32, 8>, 0x31e044b2369bf812, SC);
+    test_utils::ops::test_div!(Emulated<f32, 16>, 0x87f74cf00a528a2d, SC);
     // test_emulated!(f64, 8);
 
     // unsigned integer
diff --git a/diskann-wide/src/test_utils/ops.rs b/diskann-wide/src/test_utils/ops.rs
index fc15661f1e..69bba4bc41 100644
--- a/diskann-wide/src/test_utils/ops.rs
+++ b/diskann-wide/src/test_utils/ops.rs
@@ -425,6 +425,38 @@ macro_rules! test_mul {
     };
 }
 
+macro_rules! test_div {
+    ($wide:ident $(< $($ps:tt),+ >)?, $seed:literal, $arch:expr) => {
+        paste::paste! {
+            #[test]
+            fn [<div_ $wide:lower $(_$($ps )x+)?>]() {
+                use $crate::SIMDVector;
+                type T = $wide $(< $($ps),+>)?;
+                type Scalar = <T as SIMDVector>::Scalar;
+
+                if let Some(arch) = $arch {
+                    let f = move |a: &[Scalar], b: &[Scalar]| {
+                        let got = (
+                            <T>::from_array(arch, a.try_into().unwrap()) /
+                            <T>::from_array(arch, b.try_into().unwrap())
+                        ).to_array();
+                        test_utils::test_binary_op(
+                            &a,
+                            &b,
+                            &got,
+                            &|l: Scalar, r: Scalar| { l / r },
+                            "binary division",
+                        )
+                    };
+
+                    let n = T::LANES;
+                    $crate::test_utils::driver::drive_binary(&f, (n, n), $seed);
+                }
+            }
+        }
+    };
+}
+
 macro_rules! test_fma {
     ($wide:ident $(< $($ps:tt),+ >)?, $seed:literal, $arch:expr) => {
         paste::paste! {
@@ -1141,6 +1173,7 @@ pub(crate) use test_add;
 pub(crate) use test_bitops;
 pub(crate) use test_cast;
 pub(crate) use test_cmp;
+pub(crate) use test_div;
 pub(crate) use test_fma;
 pub(crate) use test_lossless_convert;
 pub(crate) use test_minmax;
diff --git a/diskann-wide/src/traits.rs b/diskann-wide/src/traits.rs
index 09150f0c7d..a233622615 100644
--- a/diskann-wide/src/traits.rs
+++ b/diskann-wide/src/traits.rs
@@ -28,7 +28,7 @@ use super::{
 /// - <https://users.rust-lang.org/t/limitation-of-associated-const-in-traits/73491/2>
 /// - <https://github.com/rust-lang/rust/issues/76560>
 pub trait ArrayType<T>: SupportedLaneCount {
-    type Type;
+    type Type: AsRef<[T]> + AsMut<[T]>;
 }
 
 /// Map scalar + lengths to arrays.
diff --git a/diskann/Cargo.toml b/diskann/Cargo.toml
index 14f6ba0746..c0016416dc 100644
--- a/diskann/Cargo.toml
+++ b/diskann/Cargo.toml
@@ -30,12 +30,14 @@ diskann-wide = { workspace = true }
 
 # Optional Dependencies
 dashmap = { workspace = true, optional = true }
+diskann-linalg = { workspace = true, optional = true }
 
 [dev-dependencies]
 futures-util = { workspace = true, default-features = false  }
 pin-project.workspace = true
 rand.workspace = true
 relative-path = "2.0.1"
+rstest.workspace = true
 serde = { workspace = true, features = ["derive"] }
 serde_json = { workspace = true }
 tokio = { workspace = true, features = ["macros", "sync"] }
@@ -56,6 +58,9 @@ panic = "warn"
 [features]
 default = ["tracing"]
 
+# Enable PiPNN numerical kernels.
+pipnn = ["dep:diskann-linalg"]
+
 # Enable "tracing" diagnostics.
 tracing = ["dep:tracing"]
 
diff --git a/diskann/src/graph/mod.rs b/diskann/src/graph/mod.rs
index 374efc6443..b9a9fc2345 100644
--- a/diskann/src/graph/mod.rs
+++ b/diskann/src/graph/mod.rs
@@ -8,6 +8,9 @@ pub use search_output_buffer::{
     BufferState, IdDistance, IdDistanceAssociatedData, SearchOutputBuffer,
 };
 
+#[cfg(feature = "pipnn")]
+pub mod pipnn;
+
 pub mod adjacencylist;
 pub use adjacencylist::AdjacencyList;
 
diff --git a/diskann/src/graph/pipnn/kernel_metric.rs b/diskann/src/graph/pipnn/kernel_metric.rs
new file mode 100644
index 0000000000..91f40fd69b
--- /dev/null
+++ b/diskann/src/graph/pipnn/kernel_metric.rs
@@ -0,0 +1,144 @@
+/*
+ * Copyright (c) Microsoft Corporation.
+ * Licensed under the MIT license.
+ */
+
+//! This module provides metric markers and shared numerical functions.
+
+mod leaf;
+mod partition;
+
+pub(super) use leaf::LeafMetric;
+pub(super) use partition::PartitionMetric;
+
+use super::simd::PiPNNSIMDVector;
+
+pub(super) struct L2;
+pub(super) struct Cosine;
+pub(super) struct CosineNormalized;
+pub(super) struct InnerProduct;
+
+/// Prepared norms for one point stripe and its sampled leaders.
+#[derive(Clone, Copy, Debug)]
+pub(super) struct PartitionNorms<'a> {
+    pub(super) point_norms: &'a [f32],
+    pub(super) leader_norms: &'a [f32],
+}
+
+/// Compute SIMD cosine distance with the DiskANN zero-norm and NaN rules.
+///
+/// Each lane contains one point pair. A zero norm produces zero similarity.
+/// Finite similarity is clamped to the cosine range before distance conversion.
+#[inline(always)]
+pub(super) fn cosine_distance_simd<F>(arch: F::Arch, dot: F, source_norm: F, target_norm: F) -> F
+where
+    F: PiPNNSIMDVector,
+{
+    let zero = F::default(arch);
+    let one = F::splat(arch, 1.0);
+    let minimum_norm = F::splat(arch, f32::MIN_POSITIVE.sqrt());
+    let source_zero = source_norm.lt_simd(minimum_norm);
+    let target_zero = target_norm.lt_simd(minimum_norm);
+    let denominator = source_norm * target_norm;
+    let safe_denominator = F::select(source_zero, one, F::select(target_zero, one, denominator));
+    let cosine = F::select(
+        source_zero,
+        zero,
+        F::select(target_zero, zero, dot / safe_denominator),
+    );
+    let negative_one = F::splat(arch, -1.0);
+    one - negative_one.max_simd(cosine.min_simd(one))
+}
+
+/// Compute one cosine distance with the DiskANN zero-norm and NaN rules.
+#[inline(always)]
+pub(super) fn cosine_distance_single(dot: f32, source_norm: f32, target_norm: f32) -> f32 {
+    if source_norm < f32::MIN_POSITIVE.sqrt() || target_norm < f32::MIN_POSITIVE.sqrt() {
+        1.0
+    } else {
+        let cosine = dot / (source_norm * target_norm);
+        1.0 - (-1.0_f32).max(1.0_f32.min(cosine))
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::cosine_distance_single;
+
+    mod cosine_distance_single_tests {
+        use super::cosine_distance_single;
+
+        #[test]
+        fn cosine_zero_source_norm_produces_unit_distance() {
+            // Given
+            let source_norm = 0.0;
+            let zero_norm_similarity = 0.0;
+            let expected_one_minus_zero_similarity = 1.0 - zero_norm_similarity;
+
+            // When
+            let actual_distance = cosine_distance_single(100.0, source_norm, 2.0);
+
+            // Then
+            assert_eq!(actual_distance, expected_one_minus_zero_similarity);
+        }
+
+        #[test]
+        fn cosine_zero_target_norm_produces_unit_distance() {
+            // Given
+            let target_norm = 0.0;
+            let zero_norm_similarity = 0.0;
+            let expected_one_minus_zero_similarity = 1.0 - zero_norm_similarity;
+
+            // When
+            let actual_distance = cosine_distance_single(-100.0, 2.0, target_norm);
+
+            // Then
+            assert_eq!(actual_distance, expected_one_minus_zero_similarity);
+        }
+
+        #[test]
+        fn cosine_similarity_above_one_clamps_to_zero_distance() {
+            // Given
+            let dot_product_above_valid_similarity = 5.0;
+            let maximum_cosine_similarity = 1.0;
+            let expected_one_minus_maximum_similarity = 1.0 - maximum_cosine_similarity;
+
+            // When
+            let actual_distance =
+                cosine_distance_single(dot_product_above_valid_similarity, 2.0, 2.0);
+
+            // Then
+            assert_eq!(actual_distance, expected_one_minus_maximum_similarity);
+        }
+
+        #[test]
+        fn cosine_similarity_below_negative_one_clamps_to_distance_two() {
+            // Given
+            let dot_product_below_valid_similarity = -5.0;
+            let minimum_cosine_similarity = -1.0;
+            let expected_one_minus_minimum_similarity = 1.0 - minimum_cosine_similarity;
+
+            // When
+            let actual_distance =
+                cosine_distance_single(dot_product_below_valid_similarity, 2.0, 2.0);
+
+            // Then
+            assert_eq!(actual_distance, expected_one_minus_minimum_similarity);
+        }
+
+        #[test]
+        fn cosine_nan_similarity_follows_the_diskann_min_max_rule() {
+            // Given
+            let nan_dot_product = f32::NAN;
+            // `f32::min` keeps its finite operand when the other operand is NaN.
+            let finite_operand_selected_by_min = 1.0;
+            let expected_one_minus_selected_operand = 1.0 - finite_operand_selected_by_min;
+
+            // When
+            let actual_distance = cosine_distance_single(nan_dot_product, 1.0, 1.0);
+
+            // Then
+            assert_eq!(actual_distance, expected_one_minus_selected_operand);
+        }
+    }
+}
diff --git a/diskann/src/graph/pipnn/kernel_metric/leaf.rs b/diskann/src/graph/pipnn/kernel_metric/leaf.rs
new file mode 100644
index 0000000000..844e84b401
--- /dev/null
+++ b/diskann/src/graph/pipnn/kernel_metric/leaf.rs
@@ -0,0 +1,348 @@
+/*
+ * Copyright (c) Microsoft Corporation.
+ * Licensed under the MIT license.
+ */
+
+use diskann_utils::views::MatrixView;
+use diskann_wide::{SIMDMinMax, SIMDMulAdd, SIMDVector};
+
+use super::super::simd::{PiPNNSIMDSchema, PiPNNSIMDVector};
+use super::{
+    Cosine, CosineNormalized, InnerProduct, L2, cosine_distance_simd, cosine_distance_single,
+};
+
+/// Compute leaf distances for one concrete metric.
+pub(in super::super) trait LeafMetric: Send + Sync + 'static {
+    /// SIMD representation for leaf distance scores.
+    type Simd<A>: PiPNNSIMDVector<Arch = A>
+    where
+        A: PiPNNSIMDSchema;
+
+    /// Prepare one contiguous metric-specific norm for each leaf-local point.
+    fn prepare_leaf_norms(_dots: MatrixView<'_, f32>, norms: &mut Vec<f32>) {
+        norms.clear();
+    }
+
+    /// Prepare one source norm for reuse across SIMD target groups.
+    #[inline(always)]
+    fn source_simd<A>(arch: A, _norms: &[f32], _source: usize) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema,
+    {
+        Self::Simd::<A>::default(arch)
+    }
+
+    /// Prepare one source norm for reuse across single target values.
+    #[inline(always)]
+    fn source_single(_norms: &[f32], _source: usize) -> f32 {
+        0.0
+    }
+
+    /// Compute distances for one complete SIMD group.
+    fn distances_simd<A>(
+        arch: A,
+        norms: &[f32],
+        source_norms: Self::Simd<A>,
+        dot_products: Self::Simd<A>,
+        first_target: usize,
+    ) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema;
+
+    /// Compute one distance outside the complete SIMD prefix.
+    fn distance_single(norms: &[f32], source_norm: f32, dot_product: f32, target: usize) -> f32;
+}
+
+/// Load one complete SIMD group of prepared norms.
+#[inline(always)]
+fn load_norms_simd<F>(arch: F::Arch, norms: &[f32], first_norm: usize) -> F
+where
+    F: PiPNNSIMDVector,
+{
+    let last_norm = first_norm + F::LANES;
+    let norm_group = &norms[first_norm..last_norm];
+
+    // SAFETY: `norm_group` contains one complete SIMD group.
+    unsafe { F::load_simd(arch, norm_group.as_ptr()) }
+}
+
+impl LeafMetric for L2 {
+    type Simd<A>
+        = A::LeafScore
+    where
+        A: PiPNNSIMDSchema;
+
+    fn prepare_leaf_norms(dots: MatrixView<'_, f32>, norms: &mut Vec<f32>) {
+        norms.resize(dots.nrows(), 0.0);
+        for (point, norm) in norms.iter_mut().enumerate() {
+            *norm = dots[(point, point)];
+        }
+    }
+
+    #[inline(always)]
+    fn source_simd<A>(arch: A, norms: &[f32], source: usize) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema,
+    {
+        Self::Simd::<A>::splat(arch, norms[source])
+    }
+
+    #[inline(always)]
+    fn source_single(norms: &[f32], source: usize) -> f32 {
+        norms[source]
+    }
+
+    #[inline(always)]
+    fn distances_simd<A>(
+        arch: A,
+        norms: &[f32],
+        source_norms: Self::Simd<A>,
+        dot_products: Self::Simd<A>,
+        first_target: usize,
+    ) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema,
+    {
+        let target_norms = load_norms_simd::<Self::Simd<A>>(arch, norms, first_target);
+        (Self::Simd::<A>::splat(arch, -2.0).mul_add_simd(dot_products, source_norms) + target_norms)
+            .max_simd(Self::Simd::<A>::default(arch))
+    }
+
+    #[inline(always)]
+    fn distance_single(norms: &[f32], source_norm: f32, dot_product: f32, target: usize) -> f32 {
+        ((-2.0_f32).mul_add(dot_product, source_norm) + norms[target]).max(0.0)
+    }
+}
+
+impl LeafMetric for Cosine {
+    type Simd<A>
+        = A::LeafScore
+    where
+        A: PiPNNSIMDSchema;
+
+    fn prepare_leaf_norms(dots: MatrixView<'_, f32>, norms: &mut Vec<f32>) {
+        norms.resize(dots.nrows(), 0.0);
+        for (point, norm) in norms.iter_mut().enumerate() {
+            *norm = dots[(point, point)].sqrt();
+        }
+    }
+
+    #[inline(always)]
+    fn source_simd<A>(arch: A, norms: &[f32], source: usize) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema,
+    {
+        Self::Simd::<A>::splat(arch, norms[source])
+    }
+
+    #[inline(always)]
+    fn source_single(norms: &[f32], source: usize) -> f32 {
+        norms[source]
+    }
+
+    #[inline(always)]
+    fn distances_simd<A>(
+        arch: A,
+        norms: &[f32],
+        source_norms: Self::Simd<A>,
+        dot_products: Self::Simd<A>,
+        first_target: usize,
+    ) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema,
+    {
+        let target_norms = load_norms_simd::<Self::Simd<A>>(arch, norms, first_target);
+        cosine_distance_simd(arch, dot_products, source_norms, target_norms)
+            .max_simd(Self::Simd::<A>::default(arch))
+    }
+
+    #[inline(always)]
+    fn distance_single(norms: &[f32], source_norm: f32, dot_product: f32, target: usize) -> f32 {
+        cosine_distance_single(dot_product, source_norm, norms[target]).max(0.0)
+    }
+}
+
+impl LeafMetric for CosineNormalized {
+    type Simd<A>
+        = A::LeafScore
+    where
+        A: PiPNNSIMDSchema;
+
+    #[inline(always)]
+    fn distances_simd<A>(
+        arch: A,
+        _norms: &[f32],
+        _source_norms: Self::Simd<A>,
+        dot_products: Self::Simd<A>,
+        _first_target: usize,
+    ) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema,
+    {
+        Self::Simd::<A>::splat(arch, 1.0) - dot_products
+    }
+
+    #[inline(always)]
+    fn distance_single(_norms: &[f32], _source_norm: f32, dot_product: f32, _target: usize) -> f32 {
+        1.0 - dot_product
+    }
+}
+
+impl LeafMetric for InnerProduct {
+    type Simd<A>
+        = A::LeafScore
+    where
+        A: PiPNNSIMDSchema;
+
+    #[inline(always)]
+    fn distances_simd<A>(
+        arch: A,
+        _norms: &[f32],
+        _source_norms: Self::Simd<A>,
+        dot_products: Self::Simd<A>,
+        _first_target: usize,
+    ) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema,
+    {
+        Self::Simd::<A>::default(arch) - dot_products
+    }
+
+    #[inline(always)]
+    fn distance_single(_norms: &[f32], _source_norm: f32, dot_product: f32, _target: usize) -> f32 {
+        -dot_product
+    }
+}
+
+#[cfg(test)]
+#[allow(
+    clippy::unwrap_used,
+    reason = "deterministic test matrices must abort on invalid setup"
+)]
+mod tests {
+    use super::*;
+
+    mod prepare_leaf_norms_tests {
+        use super::*;
+
+        #[test]
+        fn l2_leaf_norms_equal_the_gram_diagonal() {
+            // Given
+            let first_squared_norm = 4.0_f32;
+            let second_squared_norm = 9.0_f32;
+            let lower_gram_values = [first_squared_norm, 0.0, 0.0, second_squared_norm];
+            let lower_gram = MatrixView::try_from(&lower_gram_values[..], 2, 2).unwrap();
+            let expected_gram_diagonal = [first_squared_norm, second_squared_norm];
+            let mut actual_norms = Vec::new();
+
+            // When
+            L2::prepare_leaf_norms(lower_gram, &mut actual_norms);
+
+            // Then
+            assert_eq!(actual_norms, expected_gram_diagonal);
+        }
+
+        #[test]
+        fn cosine_leaf_norms_equal_square_roots_of_the_gram_diagonal() {
+            // Given
+            let first_squared_norm = 4.0_f32;
+            let second_squared_norm = 9.0_f32;
+            let lower_gram_values = [first_squared_norm, 0.0, 0.0, second_squared_norm];
+            let lower_gram = MatrixView::try_from(&lower_gram_values[..], 2, 2).unwrap();
+            let expected_square_roots_of_diagonal =
+                [first_squared_norm.sqrt(), second_squared_norm.sqrt()];
+            let mut actual_norms = Vec::new();
+
+            // When
+            Cosine::prepare_leaf_norms(lower_gram, &mut actual_norms);
+
+            // Then
+            assert_eq!(actual_norms, expected_square_roots_of_diagonal);
+        }
+    }
+
+    mod distance_single_tests {
+        use super::*;
+
+        #[test]
+        fn l2_distance_equals_squared_norm_sum_minus_twice_the_dot_product() {
+            // Given
+            let source_squared_norm = 4.0;
+            let target_squared_norm = 9.0;
+            let dot_product = 6.0;
+            let squared_norms = [source_squared_norm, target_squared_norm];
+            let expected_squared_l2_distance =
+                source_squared_norm + target_squared_norm - 2.0 * dot_product;
+
+            // When
+            let actual_distance =
+                L2::distance_single(&squared_norms, source_squared_norm, dot_product, 1);
+
+            // Then
+            assert_eq!(actual_distance, expected_squared_l2_distance);
+        }
+
+        #[test]
+        fn l2_clamps_negative_roundoff_to_zero() {
+            // Given
+            let source_squared_norm = 1.0;
+            let target_squared_norm = 1.0;
+            let dot_product_above_exact_norm = 1.000_001;
+            let squared_norms = [source_squared_norm, target_squared_norm];
+            let expected_non_negative_distance = 0.0;
+
+            // When
+            let actual_distance = L2::distance_single(
+                &squared_norms,
+                source_squared_norm,
+                dot_product_above_exact_norm,
+                1,
+            );
+
+            // Then
+            assert_eq!(actual_distance, expected_non_negative_distance);
+        }
+
+        #[test]
+        fn cosine_distance_equals_one_minus_dot_over_norm_product() {
+            // Given
+            let source_norm = 2.0;
+            let target_norm = 4.0;
+            let dot_product = 4.0;
+            let norms = [source_norm, target_norm];
+            let expected_one_minus_normalized_dot = 1.0 - dot_product / (source_norm * target_norm);
+
+            // When
+            let actual_distance = Cosine::distance_single(&norms, source_norm, dot_product, 1);
+
+            // Then
+            assert_eq!(actual_distance, expected_one_minus_normalized_dot);
+        }
+
+        #[test]
+        fn normalized_cosine_distance_is_one_minus_the_dot_product() {
+            // Given
+            let dot_product = 0.25;
+            let expected_one_minus_dot = 1.0 - dot_product;
+
+            // When
+            let actual_distance = CosineNormalized::distance_single(&[], 0.0, dot_product, 0);
+
+            // Then
+            assert_eq!(actual_distance, expected_one_minus_dot);
+        }
+
+        #[test]
+        fn inner_product_distance_is_the_negative_dot_product() {
+            // Given
+            let dot_product = 3.0;
+            let expected_negative_dot = -dot_product;
+
+            // When
+            let actual_distance = InnerProduct::distance_single(&[], 0.0, dot_product, 0);
+
+            // Then
+            assert_eq!(actual_distance, expected_negative_dot);
+        }
+    }
+}
diff --git a/diskann/src/graph/pipnn/kernel_metric/partition.rs b/diskann/src/graph/pipnn/kernel_metric/partition.rs
new file mode 100644
index 0000000000..05b54fe9b3
--- /dev/null
+++ b/diskann/src/graph/pipnn/kernel_metric/partition.rs
@@ -0,0 +1,400 @@
+/*
+ * Copyright (c) Microsoft Corporation.
+ * Licensed under the MIT license.
+ */
+
+use diskann_utils::views::MatrixView;
+use diskann_vector::{Norm, norm::FastL2NormSquared};
+use diskann_wide::{SIMDMulAdd, SIMDVector};
+
+use super::super::simd::{PiPNNSIMDSchema, PiPNNSIMDVector};
+use super::{
+    Cosine, CosineNormalized, InnerProduct, L2, PartitionNorms, cosine_distance_simd,
+    cosine_distance_single,
+};
+
+/// Compute partition rankings for one concrete metric.
+pub(in super::super) trait PartitionMetric: Send + Sync + 'static {
+    /// SIMD representation for partition ranking scores.
+    type Simd<A>: PiPNNSIMDVector<Arch = A>
+    where
+        A: PiPNNSIMDSchema;
+
+    /// Prepare one norm value for each point in the active stripe.
+    fn prepare_point_norms(_points: MatrixView<'_, f32>, norms: &mut Vec<f32>) {
+        norms.clear();
+    }
+
+    /// Prepare one norm value for each sampled leader.
+    fn prepare_leader_norms(_leaders: MatrixView<'_, f32>, norms: &mut Vec<f32>) {
+        norms.clear();
+    }
+
+    /// Prepare one point norm for reuse across SIMD leader groups.
+    #[inline(always)]
+    fn point_simd<A>(arch: A, _norms: PartitionNorms<'_>, _point: usize) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema,
+    {
+        Self::Simd::<A>::default(arch)
+    }
+
+    /// Prepare one point norm for reuse across single leader values.
+    #[inline(always)]
+    fn point_single(_norms: PartitionNorms<'_>, _point: usize) -> f32 {
+        0.0
+    }
+
+    /// Compute rankings for one complete SIMD group.
+    fn rankings_simd<A>(
+        arch: A,
+        norms: PartitionNorms<'_>,
+        point_norms: Self::Simd<A>,
+        dot_products: Self::Simd<A>,
+        first_leader: usize,
+    ) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema;
+
+    /// Compute one ranking outside the complete SIMD prefix.
+    fn ranking_single(
+        norms: PartitionNorms<'_>,
+        point_norm: f32,
+        dot_product: f32,
+        leader: usize,
+    ) -> f32;
+}
+
+/// Load one complete SIMD group of prepared norms.
+#[inline(always)]
+fn load_norms_simd<F>(arch: F::Arch, norms: &[f32], first_norm: usize) -> F
+where
+    F: PiPNNSIMDVector,
+{
+    let last_norm = first_norm + F::LANES;
+    let norm_group = &norms[first_norm..last_norm];
+
+    // SAFETY: `norm_group` contains one complete SIMD group.
+    unsafe { F::load_simd(arch, norm_group.as_ptr()) }
+}
+
+impl PartitionMetric for L2 {
+    type Simd<A>
+        = A::PartitionScore
+    where
+        A: PiPNNSIMDSchema;
+
+    fn prepare_leader_norms(leaders: MatrixView<'_, f32>, norms: &mut Vec<f32>) {
+        norms.resize(leaders.nrows(), 0.0);
+        for (norm, leader) in norms.iter_mut().zip(leaders.row_iter()) {
+            *norm = leader.iter().map(|value| value * value).sum();
+        }
+    }
+
+    #[inline(always)]
+    fn rankings_simd<A>(
+        arch: A,
+        norms: PartitionNorms<'_>,
+        _point_norms: Self::Simd<A>,
+        dot_products: Self::Simd<A>,
+        first_leader: usize,
+    ) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema,
+    {
+        let leader_norms = load_norms_simd::<Self::Simd<A>>(arch, norms.leader_norms, first_leader);
+        Self::Simd::<A>::splat(arch, -2.0).mul_add_simd(dot_products, leader_norms)
+    }
+
+    #[inline(always)]
+    fn ranking_single(
+        norms: PartitionNorms<'_>,
+        _point_norm: f32,
+        dot_product: f32,
+        leader: usize,
+    ) -> f32 {
+        (-2.0_f32).mul_add(dot_product, norms.leader_norms[leader])
+    }
+}
+
+impl PartitionMetric for Cosine {
+    type Simd<A>
+        = A::PartitionScore
+    where
+        A: PiPNNSIMDSchema;
+
+    fn prepare_point_norms(points: MatrixView<'_, f32>, norms: &mut Vec<f32>) {
+        norms.resize(points.nrows(), 0.0);
+        for (norm, point) in norms.iter_mut().zip(points.row_iter()) {
+            *norm = FastL2NormSquared.evaluate(point).sqrt();
+        }
+    }
+
+    fn prepare_leader_norms(leaders: MatrixView<'_, f32>, norms: &mut Vec<f32>) {
+        norms.resize(leaders.nrows(), 0.0);
+        for (norm, leader) in norms.iter_mut().zip(leaders.row_iter()) {
+            *norm = leader.iter().map(|value| value * value).sum::<f32>().sqrt();
+        }
+    }
+
+    #[inline(always)]
+    fn point_simd<A>(arch: A, norms: PartitionNorms<'_>, point: usize) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema,
+    {
+        Self::Simd::<A>::splat(arch, norms.point_norms[point])
+    }
+
+    #[inline(always)]
+    fn point_single(norms: PartitionNorms<'_>, point: usize) -> f32 {
+        norms.point_norms[point]
+    }
+
+    #[inline(always)]
+    fn rankings_simd<A>(
+        arch: A,
+        norms: PartitionNorms<'_>,
+        point_norms: Self::Simd<A>,
+        dot_products: Self::Simd<A>,
+        first_leader: usize,
+    ) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema,
+    {
+        let leader_norms = load_norms_simd::<Self::Simd<A>>(arch, norms.leader_norms, first_leader);
+        cosine_distance_simd(arch, dot_products, point_norms, leader_norms)
+    }
+
+    #[inline(always)]
+    fn ranking_single(
+        norms: PartitionNorms<'_>,
+        point_norm: f32,
+        dot_product: f32,
+        leader: usize,
+    ) -> f32 {
+        cosine_distance_single(dot_product, point_norm, norms.leader_norms[leader])
+    }
+}
+
+impl PartitionMetric for CosineNormalized {
+    type Simd<A>
+        = A::PartitionScore
+    where
+        A: PiPNNSIMDSchema;
+
+    #[inline(always)]
+    fn rankings_simd<A>(
+        arch: A,
+        _norms: PartitionNorms<'_>,
+        _point_norms: Self::Simd<A>,
+        dot_products: Self::Simd<A>,
+        _first_leader: usize,
+    ) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema,
+    {
+        Self::Simd::<A>::splat(arch, 1.0) - dot_products
+    }
+
+    #[inline(always)]
+    fn ranking_single(
+        _norms: PartitionNorms<'_>,
+        _point_norm: f32,
+        dot_product: f32,
+        _leader: usize,
+    ) -> f32 {
+        1.0 - dot_product
+    }
+}
+
+impl PartitionMetric for InnerProduct {
+    type Simd<A>
+        = A::PartitionScore
+    where
+        A: PiPNNSIMDSchema;
+
+    #[inline(always)]
+    fn rankings_simd<A>(
+        arch: A,
+        _norms: PartitionNorms<'_>,
+        _point_norms: Self::Simd<A>,
+        dot_products: Self::Simd<A>,
+        _first_leader: usize,
+    ) -> Self::Simd<A>
+    where
+        A: PiPNNSIMDSchema,
+    {
+        Self::Simd::<A>::default(arch) - dot_products
+    }
+
+    #[inline(always)]
+    fn ranking_single(
+        _norms: PartitionNorms<'_>,
+        _point_norm: f32,
+        dot_product: f32,
+        _leader: usize,
+    ) -> f32 {
+        -dot_product
+    }
+}
+
+#[cfg(test)]
+#[allow(
+    clippy::unwrap_used,
+    reason = "deterministic test matrices must abort on invalid setup"
+)]
+mod tests {
+    use super::*;
+
+    fn rank_single_leader<M: PartitionMetric>(
+        dot_product: f32,
+        point_norm: f32,
+        leader_norm: f32,
+    ) -> f32 {
+        let point_norms = [point_norm];
+        let leader_norms = [leader_norm];
+        let norms = PartitionNorms {
+            point_norms: &point_norms,
+            leader_norms: &leader_norms,
+        };
+        M::ranking_single(norms, M::point_single(norms, 0), dot_product, 0)
+    }
+
+    mod prepare_leader_norms_tests {
+        use super::*;
+
+        #[test]
+        fn l2_leader_norm_is_the_sum_of_squared_components() {
+            // Given
+            let first_leader = [1.0_f32, 2.0];
+            let second_leader = [3.0_f32, 4.0];
+            let leader_values = [
+                first_leader[0],
+                first_leader[1],
+                second_leader[0],
+                second_leader[1],
+            ];
+            let leaders = MatrixView::try_from(&leader_values[..], 2, 2).unwrap();
+            let expected_row_squared_norms = [
+                first_leader[0].powi(2) + first_leader[1].powi(2),
+                second_leader[0].powi(2) + second_leader[1].powi(2),
+            ];
+            let mut actual_norms = Vec::new();
+
+            // When
+            L2::prepare_leader_norms(leaders, &mut actual_norms);
+
+            // Then
+            assert_eq!(actual_norms, expected_row_squared_norms);
+        }
+
+        #[test]
+        fn cosine_leader_norm_is_the_square_root_of_the_squared_component_sum() {
+            // Given
+            let first_leader = [1.0_f32, 2.0];
+            let second_leader = [3.0_f32, 4.0];
+            let leader_values = [
+                first_leader[0],
+                first_leader[1],
+                second_leader[0],
+                second_leader[1],
+            ];
+            let leaders = MatrixView::try_from(&leader_values[..], 2, 2).unwrap();
+            let expected_row_norms = [
+                (first_leader[0].powi(2) + first_leader[1].powi(2)).sqrt(),
+                (second_leader[0].powi(2) + second_leader[1].powi(2)).sqrt(),
+            ];
+            let mut actual_norms = Vec::new();
+
+            // When
+            Cosine::prepare_leader_norms(leaders, &mut actual_norms);
+
+            // Then
+            assert_eq!(actual_norms, expected_row_norms);
+        }
+    }
+
+    #[test]
+    fn cosine_point_norm_is_the_square_root_of_the_squared_component_sum() {
+        // Given
+        let first_point = [1.0_f32, 2.0];
+        let second_point = [3.0_f32, 4.0];
+        let point_values = [
+            first_point[0],
+            first_point[1],
+            second_point[0],
+            second_point[1],
+        ];
+        let points = MatrixView::try_from(&point_values[..], 2, 2).unwrap();
+        let expected_row_norms = [
+            (first_point[0].powi(2) + first_point[1].powi(2)).sqrt(),
+            (second_point[0].powi(2) + second_point[1].powi(2)).sqrt(),
+        ];
+        let mut actual_norms = Vec::new();
+
+        // When
+        Cosine::prepare_point_norms(points, &mut actual_norms);
+
+        // Then
+        assert_eq!(actual_norms, expected_row_norms);
+    }
+
+    mod ranking_single_tests {
+        use super::*;
+
+        #[test]
+        fn l2_ranking_equals_leader_squared_norm_minus_twice_the_dot_product() {
+            // Given
+            let dot_product = 2.0;
+            let leader_squared_norm = 9.0;
+            let expected_leader_norm_minus_twice_dot = leader_squared_norm - 2.0 * dot_product;
+
+            // When
+            let actual_ranking = rank_single_leader::<L2>(dot_product, 0.0, leader_squared_norm);
+
+            // Then
+            assert_eq!(actual_ranking, expected_leader_norm_minus_twice_dot);
+        }
+
+        #[test]
+        fn cosine_ranking_equals_one_minus_dot_over_norm_product() {
+            // Given
+            let dot_product = 4.0;
+            let point_norm = 2.0;
+            let leader_norm = 4.0;
+            let expected_one_minus_normalized_dot = 1.0 - dot_product / (point_norm * leader_norm);
+
+            // When
+            let actual_ranking = rank_single_leader::<Cosine>(dot_product, point_norm, leader_norm);
+
+            // Then
+            assert_eq!(actual_ranking, expected_one_minus_normalized_dot);
+        }
+
+        #[test]
+        fn normalized_cosine_ranking_is_one_minus_the_dot_product() {
+            // Given
+            let dot_product = 0.25;
+            let expected_one_minus_dot = 1.0 - dot_product;
+
+            // When
+            let actual_ranking = rank_single_leader::<CosineNormalized>(dot_product, 0.0, 0.0);
+
+            // Then
+            assert_eq!(actual_ranking, expected_one_minus_dot);
+        }
+
+        #[test]
+        fn inner_product_ranking_is_the_negative_dot_product() {
+            // Given
+            let dot_product = 3.0;
+            let expected_negative_dot = -dot_product;
+
+            // When
+            let actual_ranking = rank_single_leader::<InnerProduct>(dot_product, 0.0, 0.0);
+
+            // Then
+            assert_eq!(actual_ranking, expected_negative_dot);
+        }
+    }
+}
diff --git a/diskann/src/graph/pipnn/leaf_kernel.rs b/diskann/src/graph/pipnn/leaf_kernel.rs
new file mode 100644
index 0000000000..49a4426631
--- /dev/null
+++ b/diskann/src/graph/pipnn/leaf_kernel.rs
@@ -0,0 +1,1283 @@
+/*
+ * Copyright (c) Microsoft Corporation.
+ * Licensed under the MIT license.
+ */
+
+//! Leaf-local top-k selection from packed `f32` point vectors.
+//!
+//! The kernel computes the lower-triangular Gram matrix and metric-specific
+//! norms. Its ranking loop reads each strict-lower point pair once and updates
+//! both points.
+//!
+//! The output is an `n × k` matrix of sorted [`LeafNeighbor`] values. Each target
+//! is a position in the leaf. Widths 1 through 3 use fixed insertion. Larger
+//! widths use the runtime insertion loop.
+//!
+//! Strict comparisons keep scan order for equal distances. They do not rank NaN.
+//! An unfilled output slot contains [`LeafNeighbor::default`]. All supported
+//! metrics use the same SIMD-group and single-value traversal.
+//!
+//! The caller supplies concrete architecture `A` and metric `M`. The private
+//! dot ranker receives the square matrix created by this module.
+//! [`LeafKernelWorkspace`] stores reusable numerical scratch.
+
+use crate::{ANNError, ANNResult};
+use diskann_utils::views::{MatrixView, MutMatrixView};
+use diskann_wide::{SIMDPartialOrd, SIMDVector};
+
+use super::{
+    kernel_metric::LeafMetric,
+    simd::{PiPNNSIMDSchema, PiPNNSIMDVector},
+};
+
+/// One leaf-local neighbor and its metric distance.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub(super) struct LeafNeighbor {
+    /// Target position in the leaf, not a dataset ID.
+    pub(super) target: u32,
+    /// Distance from the source point to `target`.
+    pub(super) distance: f32,
+}
+
+impl LeafNeighbor {
+    /// Construct a leaf-local neighbor.
+    ///
+    /// `target` is a position in the leaf. `distance` is its score relative to
+    /// the source of the output row.
+    pub(super) const fn new(target: u32, distance: f32) -> Self {
+        Self { target, distance }
+    }
+
+    /// Return true when this slot contains a rankable leaf-local target.
+    pub(super) const fn is_assigned(self) -> bool {
+        self.target != u32::MAX
+    }
+}
+
+impl Default for LeafNeighbor {
+    fn default() -> Self {
+        Self::new(u32::MAX, f32::INFINITY)
+    }
+}
+
+/// Reusable storage for one leaf numerical pipeline.
+#[derive(Debug, Default)]
+pub(super) struct LeafKernelWorkspace {
+    dot_scratch: Vec<f32>,
+    norm_scratch: Vec<f32>,
+    worst: Vec<f32>,
+}
+
+/// Validation error returned by the dot-ranking loop.
+#[derive(Clone, Copy, Debug, PartialEq, Eq, thiserror::Error)]
+pub(super) enum LeafKernelError {
+    /// A source requests more neighbors than the leaf or fixed kernel supports.
+    #[error("invalid leaf neighbor count {neighbors} for {points} points; maximum is {maximum}")]
+    InvalidNeighborCount {
+        points: usize,
+        neighbors: usize,
+        maximum: usize,
+    },
+}
+
+/// Return the non-self neighbor count for one leaf.
+///
+/// `points` is the number of points in the leaf. `requested_k` is the configured
+/// neighbor count. The result is `min(requested_k, points - 1)`.
+///
+pub(super) fn leaf_neighbor_count(points: usize, requested_k: usize) -> usize {
+    requested_k.min(points.saturating_sub(1))
+}
+
+/// Compute local nearest neighbors for one packed leaf matrix.
+///
+/// # Errors
+///
+/// Returns an error for invalid linear-algebra input or output width.
+pub(super) fn select_leaf_neighbors<A, M>(
+    arch: A,
+    points: MatrixView<'_, f32>,
+    output: MutMatrixView<'_, LeafNeighbor>,
+    workspace: &mut LeafKernelWorkspace,
+) -> ANNResult<()>
+where
+    A: PiPNNSIMDSchema,
+    M: LeafMetric,
+{
+    let point_count = points.nrows();
+    let dot_count = point_count * point_count;
+    let LeafKernelWorkspace {
+        dot_scratch,
+        norm_scratch,
+        worst,
+    } = workspace;
+    if dot_scratch.len() < dot_count {
+        dot_scratch.resize(dot_count, 0.0);
+    }
+    diskann_linalg::sgemm_aat_lower(
+        point_count,
+        points.ncols(),
+        points.as_slice(),
+        &mut dot_scratch[..dot_count],
+    )
+    .map_err(ANNError::new)?;
+    let dots = MatrixView::try_from(&dot_scratch[..dot_count], point_count, point_count)
+        .map_err(|error| ANNError::new(error.as_static()))?;
+    M::prepare_leaf_norms(dots, norm_scratch);
+    rank_leaf_dots::<A, M>(arch, dots, norm_scratch, output, worst).map_err(ANNError::new)
+}
+
+/// Rank a prepared lower-triangular Gram matrix.
+fn rank_leaf_dots<A, M>(
+    arch: A,
+    input: MatrixView<'_, f32>,
+    norms: &[f32],
+    mut output: MutMatrixView<'_, LeafNeighbor>,
+    worst: &mut Vec<f32>,
+) -> Result<(), LeafKernelError>
+where
+    A: PiPNNSIMDSchema,
+    M: LeafMetric,
+{
+    validate_neighbor_count(input, &output)?;
+    let neighbor_count = output.ncols();
+    if neighbor_count == 0 {
+        return Ok(());
+    }
+
+    worst.resize(input.nrows(), f32::INFINITY);
+    output.as_mut_slice().fill(LeafNeighbor::default());
+    worst.fill(f32::INFINITY);
+
+    match neighbor_count {
+        1 => scan_fixed_width::<A, M, 1>(arch, input, norms, output.as_mut_slice(), worst),
+        2 => scan_fixed_width::<A, M, 2>(arch, input, norms, output.as_mut_slice(), worst),
+        3 => scan_fixed_width::<A, M, 3>(arch, input, norms, output.as_mut_slice(), worst),
+        _ => scan_runtime_width::<A, M>(
+            arch,
+            input,
+            norms,
+            output.as_mut_slice(),
+            neighbor_count,
+            worst,
+        ),
+    }
+    Ok(())
+}
+
+/// Check the safety conditions for the SIMD kernel.
+///
+/// Check the output width against the number of non-self points.
+fn validate_neighbor_count(
+    input: MatrixView<'_, f32>,
+    output: &MutMatrixView<'_, LeafNeighbor>,
+) -> Result<(), LeafKernelError> {
+    let point_count = input.nrows();
+    let maximum_neighbors = point_count.saturating_sub(1);
+    let neighbor_count = output.ncols();
+    if neighbor_count > maximum_neighbors {
+        return Err(LeafKernelError::InvalidNeighborCount {
+            points: point_count,
+            neighbors: neighbor_count,
+            maximum: maximum_neighbors,
+        });
+    }
+    Ok(())
+}
+
+/// Select neighbors with a fixed output width.
+fn scan_fixed_width<A, M, const N: usize>(
+    arch: A,
+    input: MatrixView<'_, f32>,
+    norms: &[f32],
+    output: &mut [LeafNeighbor],
+    worst: &mut [f32],
+) where
+    A: PiPNNSIMDSchema,
+    M: LeafMetric,
+    [LeafNeighbor; N]: SortedInsert<LeafNeighbor>,
+{
+    let (rows, _) = output.as_chunks_mut::<N>();
+    scan_point_pairs::<A, M, _>(arch, input, norms, worst, |source, target, distance| {
+        insert_neighbor(&mut rows[source], target, distance)
+    });
+}
+
+/// Select neighbors with a runtime output width.
+fn scan_runtime_width<A, M>(
+    arch: A,
+    input: MatrixView<'_, f32>,
+    norms: &[f32],
+    output: &mut [LeafNeighbor],
+    width: usize,
+    worst: &mut [f32],
+) where
+    A: PiPNNSIMDSchema,
+    M: LeafMetric,
+{
+    scan_point_pairs::<A, M, _>(arch, input, norms, worst, |source, target, distance| {
+        let first = source * width;
+        insert_neighbor(&mut output[first..first + width], target, distance)
+    });
+}
+
+/// Select neighbors from all unordered point pairs in one leaf.
+///
+/// The function reads the strict lower triangle once. It offers each distance to
+/// both endpoint lists. SIMD groups and single values preserve pair scan order.
+#[inline(never)]
+fn scan_point_pairs<A, M, I>(
+    arch: A,
+    input: MatrixView<'_, f32>,
+    norms: &[f32],
+    worst: &mut [f32],
+    mut insert: I,
+) where
+    A: PiPNNSIMDSchema,
+    M: LeafMetric,
+    I: FnMut(usize, u32, f32) -> f32,
+{
+    let point_count = input.nrows();
+    let dots = input.as_slice();
+    let worst_ptr = worst.as_mut_ptr();
+
+    for source in 1..point_count {
+        let source_start = source * point_count;
+        let source_simd = M::source_simd(arch, norms, source);
+        let source_single = M::source_single(norms, source);
+        // SAFETY: `rank_leaf_dots` created one threshold for each point.
+        let mut source_worst = unsafe { *worst_ptr.add(source) };
+        let mut target = 0;
+        let simd_prefix = source - source % M::Simd::<A>::LANES;
+
+        while target < simd_prefix {
+            // SAFETY: This complete SIMD group is in the strict-lower prefix.
+            let dot_products =
+                unsafe { M::Simd::<A>::load_simd(arch, dots.as_ptr().add(source_start + target)) };
+            let distances = M::distances_simd(arch, norms, source_simd, dot_products, target);
+            let source_eligible = distances.lt_simd(M::Simd::<A>::splat(arch, source_worst));
+            // SAFETY: The complete target group is below `source < point_count`.
+            let target_worst = unsafe { M::Simd::<A>::load_simd(arch, worst_ptr.add(target)) };
+            let target_eligible = distances.lt_simd(target_worst);
+            let source_bits = M::Simd::<A>::active_lanes(source_eligible);
+            let target_bits = M::Simd::<A>::active_lanes(target_eligible);
+
+            if source_bits | target_bits != 0 {
+                let values = distances.to_array();
+                let values = values.as_ref();
+                let mut source_bits = source_bits;
+                while source_bits != 0 {
+                    let lane = source_bits.trailing_zeros() as usize;
+                    source_bits &= source_bits - 1;
+                    let distance = values[lane];
+                    if distance < source_worst {
+                        source_worst = insert(source, (target + lane) as u32, distance);
+                    }
+                }
+
+                let mut target_bits = target_bits;
+                while target_bits != 0 {
+                    let lane = target_bits.trailing_zeros() as usize;
+                    target_bits &= target_bits - 1;
+                    let target_source = target + lane;
+                    let new_worst = insert(target_source, source as u32, values[lane]);
+                    // SAFETY: `target_source < source < worst.len()`.
+                    unsafe { *worst_ptr.add(target_source) = new_worst };
+                }
+            }
+            target += M::Simd::<A>::LANES;
+        }
+
+        while target < source {
+            // SAFETY: The target is in this source's strict-lower prefix.
+            let dot_product = unsafe { *dots.get_unchecked(source_start + target) };
+            let distance = M::distance_single(norms, source_single, dot_product, target);
+            if distance < source_worst {
+                source_worst = insert(source, target as u32, distance);
+            }
+            // SAFETY: `target < source < worst.len()`.
+            let target_worst = unsafe { *worst_ptr.add(target) };
+            if distance < target_worst {
+                let new_worst = insert(target, source as u32, distance);
+                // SAFETY: `target < source < worst.len()`.
+                unsafe { *worst_ptr.add(target) = new_worst };
+            }
+            target += 1;
+        }
+        // SAFETY: `source < worst.len()`.
+        unsafe { *worst_ptr.add(source) = source_worst };
+    }
+}
+
+/// Insert one value that precedes the current last retained value.
+///
+/// The caller checks eligibility before insertion. The returned value is the new
+/// last retained value for the next eligibility check.
+trait SortedInsert<T: Copy> {
+    fn insert_sorted_by(&mut self, value: T, precedes: impl Fn(T, T) -> bool) -> T;
+}
+
+impl<T: Copy> SortedInsert<T> for [T; 1] {
+    #[inline(always)]
+    fn insert_sorted_by(&mut self, value: T, _precedes: impl Fn(T, T) -> bool) -> T {
+        self[0] = value;
+        value
+    }
+}
+
+impl<T: Copy> SortedInsert<T> for [T; 2] {
+    #[inline(always)]
+    fn insert_sorted_by(&mut self, value: T, precedes: impl Fn(T, T) -> bool) -> T {
+        let first = self[0];
+        if precedes(value, first) {
+            self[0] = value;
+            self[1] = first;
+            first
+        } else {
+            self[1] = value;
+            value
+        }
+    }
+}
+
+impl<T: Copy> SortedInsert<T> for [T; 3] {
+    #[inline(always)]
+    fn insert_sorted_by(&mut self, value: T, precedes: impl Fn(T, T) -> bool) -> T {
+        let (first, second) = (self[0], self[1]);
+        if precedes(value, first) {
+            self[0] = value;
+            self[1] = first;
+            self[2] = second;
+            second
+        } else if precedes(value, second) {
+            self[1] = value;
+            self[2] = second;
+            second
+        } else {
+            self[2] = value;
+            value
+        }
+    }
+}
+
+impl<T: Copy> SortedInsert<T> for [T] {
+    #[inline(always)]
+    fn insert_sorted_by(&mut self, value: T, precedes: impl Fn(T, T) -> bool) -> T {
+        let last = self.len() - 1;
+        let mut slot = last;
+        while slot > 0 && precedes(value, self[slot - 1]) {
+            self[slot] = self[slot - 1];
+            slot -= 1;
+        }
+        self[slot] = value;
+        self[last]
+    }
+}
+
+/// Insert one candidate that is nearer than the current farthest neighbor.
+///
+/// Return the new farthest retained distance for the next candidate check.
+#[inline(always)]
+fn insert_neighbor<R>(neighbors: &mut R, target: u32, distance: f32) -> f32
+where
+    R: SortedInsert<LeafNeighbor> + ?Sized,
+{
+    neighbors
+        .insert_sorted_by(
+            LeafNeighbor::new(target, distance),
+            |candidate, retained| candidate.distance < retained.distance,
+        )
+        .distance
+}
+
+#[cfg(test)]
+mod tests {
+    use std::cmp::Ordering;
+
+    use super::*;
+    use crate::graph::pipnn::kernel_metric::{Cosine, CosineNormalized, InnerProduct, L2};
+    use diskann_utils::views::{MatrixView, MutMatrixView};
+    use diskann_vector::distance::Metric;
+    use diskann_wide::arch::{self, Target1};
+
+    struct KernelCall<'a> {
+        input: MatrixView<'a, f32>,
+        norms: &'a [f32],
+        output: MutMatrixView<'a, LeafNeighbor>,
+        workspace: &'a mut LeafKernelWorkspace,
+    }
+
+    struct DispatchMetric(Metric);
+
+    impl<A> Target1<A, Result<(), LeafKernelError>, KernelCall<'_>> for DispatchMetric
+    where
+        A: PiPNNSIMDSchema,
+    {
+        fn run(self, arch: A, call: KernelCall<'_>) -> Result<(), LeafKernelError> {
+            match self.0 {
+                Metric::L2 => rank_leaf_dots::<A, L2>(
+                    arch,
+                    call.input,
+                    call.norms,
+                    call.output,
+                    &mut call.workspace.worst,
+                ),
+                Metric::Cosine => rank_leaf_dots::<A, Cosine>(
+                    arch,
+                    call.input,
+                    call.norms,
+                    call.output,
+                    &mut call.workspace.worst,
+                ),
+                Metric::CosineNormalized => rank_leaf_dots::<A, CosineNormalized>(
+                    arch,
+                    call.input,
+                    call.norms,
+                    call.output,
+                    &mut call.workspace.worst,
+                ),
+                Metric::InnerProduct => rank_leaf_dots::<A, InnerProduct>(
+                    arch,
+                    call.input,
+                    call.norms,
+                    call.output,
+                    &mut call.workspace.worst,
+                ),
+            }
+        }
+    }
+
+    fn lower_gram_view(dots: &[f32], points: usize) -> MatrixView<'_, f32> {
+        MatrixView::try_from(dots, points, points).unwrap()
+    }
+
+    fn metric_norms(metric: Metric, lower_gram: MatrixView<'_, f32>) -> Vec<f32> {
+        fn prepare<M: LeafMetric>(lower_gram: MatrixView<'_, f32>) -> Vec<f32> {
+            let mut norms = Vec::new();
+            M::prepare_leaf_norms(lower_gram, &mut norms);
+            norms
+        }
+
+        match metric {
+            Metric::L2 => prepare::<L2>(lower_gram),
+            Metric::Cosine => prepare::<Cosine>(lower_gram),
+            Metric::CosineNormalized => prepare::<CosineNormalized>(lower_gram),
+            Metric::InnerProduct => prepare::<InnerProduct>(lower_gram),
+        }
+    }
+
+    fn rank_neighbors_with_workspace(
+        metric: Metric,
+        dots: &[f32],
+        points: usize,
+        requested_k: usize,
+        workspace: &mut LeafKernelWorkspace,
+    ) -> (usize, Vec<LeafNeighbor>) {
+        let leaf_k = leaf_neighbor_count(points, requested_k);
+        let lower_gram = lower_gram_view(dots, points);
+        let norms = metric_norms(metric, lower_gram);
+        let mut output = vec![LeafNeighbor::default(); points * leaf_k];
+        arch::dispatch1_no_features(
+            DispatchMetric(metric),
+            KernelCall {
+                input: lower_gram,
+                norms: &norms,
+                output: MutMatrixView::try_from(output.as_mut_slice(), points, leaf_k).unwrap(),
+                workspace,
+            },
+        )
+        .unwrap();
+        (leaf_k, output)
+    }
+
+    fn rank_neighbors(
+        metric: Metric,
+        dots: &[f32],
+        points: usize,
+        requested_k: usize,
+    ) -> (usize, Vec<LeafNeighbor>) {
+        rank_neighbors_with_workspace(
+            metric,
+            dots,
+            points,
+            requested_k,
+            &mut LeafKernelWorkspace::default(),
+        )
+    }
+
+    fn reference_distance(
+        metric: Metric,
+        dot: f32,
+        source_diagonal: f32,
+        target_diagonal: f32,
+    ) -> f32 {
+        match metric {
+            Metric::L2 => ((-2.0_f32).mul_add(dot, source_diagonal) + target_diagonal).max(0.0),
+            Metric::CosineNormalized => 1.0 - dot,
+            Metric::InnerProduct => -dot,
+            Metric::Cosine => {
+                let source_norm = source_diagonal.sqrt();
+                let target_norm = target_diagonal.sqrt();
+                if source_norm < f32::MIN_POSITIVE.sqrt() || target_norm < f32::MIN_POSITIVE.sqrt()
+                {
+                    1.0
+                } else {
+                    let similarity = dot / (source_norm * target_norm);
+                    1.0 - similarity.clamp(-1.0, 1.0)
+                }
+            }
+        }
+    }
+
+    fn reference_neighbors(
+        metric: Metric,
+        dots: &[f32],
+        points: usize,
+        requested_k: usize,
+    ) -> Vec<LeafNeighbor> {
+        let leaf_k = requested_k.min(points.saturating_sub(1));
+        let mut output = vec![LeafNeighbor::default(); points * leaf_k];
+        for source in 0..points {
+            let mut candidates = Vec::with_capacity(points.saturating_sub(1));
+            for target in 0..points {
+                if source == target {
+                    continue;
+                }
+                let (row, column) = if source > target {
+                    (source, target)
+                } else {
+                    (target, source)
+                };
+                let distance = reference_distance(
+                    metric,
+                    dots[row * points + column],
+                    dots[source * points + source],
+                    dots[target * points + target],
+                );
+                if distance.partial_cmp(&f32::INFINITY) == Some(Ordering::Less) {
+                    candidates.push(LeafNeighbor::new(target as u32, distance));
+                }
+            }
+            candidates.sort_by(|left, right| left.distance.total_cmp(&right.distance));
+            let retained = candidates.len().min(leaf_k);
+            output[source * leaf_k..source * leaf_k + retained]
+                .copy_from_slice(&candidates[..retained]);
+        }
+        output
+    }
+
+    /// Build a unit-diagonal lower Gram matrix for the lane-boundary sweep.
+    /// Similarity decreases as the point-index separation increases.
+    fn index_distance_lower_gram(points: usize) -> Vec<f32> {
+        let mut dots = vec![f32::NAN; points * points];
+        for source in 0..points {
+            dots[source * points + source] = 1.0;
+            for target in 0..source {
+                let separation = (source - target) as f32;
+                dots[source * points + target] = 1.0 - separation / points as f32;
+            }
+        }
+        dots
+    }
+
+    fn square_matrix_with_constant_diagonal(points: usize, diagonal: f32) -> Vec<f32> {
+        let mut values = vec![0.0; points * points];
+        for point in 0..points {
+            values[point * points + point] = diagonal;
+        }
+        values
+    }
+
+    mod insert_neighbor_tests {
+        use super::*;
+
+        #[test]
+        fn one_slot_insertion_replaces_the_retained_neighbor() {
+            // Given
+            let retained_neighbor = LeafNeighbor::new(1, 4.0);
+            let nearer_candidate = LeafNeighbor::new(2, 2.0);
+            let expected_neighbors = [nearer_candidate];
+            let mut actual_neighbors = [retained_neighbor];
+
+            // When
+            insert_neighbor(
+                &mut actual_neighbors,
+                nearer_candidate.target,
+                nearer_candidate.distance,
+            );
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors);
+        }
+
+        #[test]
+        fn two_slot_insertion_places_a_nearer_candidate_first() {
+            // Given
+            let nearest = LeafNeighbor::new(1, 1.0);
+            let farthest = LeafNeighbor::new(2, 3.0);
+            let nearer_candidate = LeafNeighbor::new(3, 0.5);
+            let expected_neighbors = [nearer_candidate, nearest];
+            let mut actual_neighbors = [nearest, farthest];
+
+            // When
+            insert_neighbor(
+                &mut actual_neighbors,
+                nearer_candidate.target,
+                nearer_candidate.distance,
+            );
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors);
+        }
+
+        #[test]
+        fn two_slot_insertion_places_a_middle_distance_last() {
+            // Given
+            let nearest = LeafNeighbor::new(1, 1.0);
+            let farthest = LeafNeighbor::new(2, 3.0);
+            let eligible_candidate = LeafNeighbor::new(3, 2.0);
+            let expected_neighbors = [nearest, eligible_candidate];
+            let mut actual_neighbors = [nearest, farthest];
+
+            // When
+            insert_neighbor(
+                &mut actual_neighbors,
+                eligible_candidate.target,
+                eligible_candidate.distance,
+            );
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors);
+        }
+
+        #[test]
+        fn three_slot_insertion_places_the_nearest_candidate_first() {
+            // Given
+            let nearest = LeafNeighbor::new(1, 1.0);
+            let middle = LeafNeighbor::new(2, 2.0);
+            let farthest = LeafNeighbor::new(3, 4.0);
+            let nearer_candidate = LeafNeighbor::new(4, 0.5);
+            let expected_neighbors = [nearer_candidate, nearest, middle];
+            let mut actual_neighbors = [nearest, middle, farthest];
+
+            // When
+            insert_neighbor(
+                &mut actual_neighbors,
+                nearer_candidate.target,
+                nearer_candidate.distance,
+            );
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors);
+        }
+
+        #[test]
+        fn three_slot_insertion_places_a_middle_candidate_between_neighbors() {
+            // Given
+            let nearest = LeafNeighbor::new(1, 1.0);
+            let middle = LeafNeighbor::new(2, 2.0);
+            let farthest = LeafNeighbor::new(3, 4.0);
+            let middle_candidate = LeafNeighbor::new(4, 1.5);
+            let expected_neighbors = [nearest, middle_candidate, middle];
+            let mut actual_neighbors = [nearest, middle, farthest];
+
+            // When
+            insert_neighbor(
+                &mut actual_neighbors,
+                middle_candidate.target,
+                middle_candidate.distance,
+            );
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors);
+        }
+
+        #[test]
+        fn three_slot_insertion_replaces_the_farthest_neighbor() {
+            // Given
+            let nearest = LeafNeighbor::new(1, 1.0);
+            let middle = LeafNeighbor::new(2, 2.0);
+            let farthest = LeafNeighbor::new(3, 4.0);
+            let eligible_candidate = LeafNeighbor::new(4, 3.0);
+            let expected_neighbors = [nearest, middle, eligible_candidate];
+            let mut actual_neighbors = [nearest, middle, farthest];
+
+            // When
+            insert_neighbor(
+                &mut actual_neighbors,
+                eligible_candidate.target,
+                eligible_candidate.distance,
+            );
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors);
+        }
+
+        #[test]
+        fn runtime_width_insertion_shifts_only_the_later_neighbors() {
+            // Given
+            let first = LeafNeighbor::new(1, 1.0);
+            let second = LeafNeighbor::new(2, 2.0);
+            let third = LeafNeighbor::new(3, 3.0);
+            let fourth = LeafNeighbor::new(4, 5.0);
+            let candidate = LeafNeighbor::new(5, 2.5);
+            let expected_neighbors = [first, second, candidate, third];
+            let mut actual_neighbors = [first, second, third, fourth];
+
+            // When
+            insert_neighbor(
+                actual_neighbors.as_mut_slice(),
+                candidate.target,
+                candidate.distance,
+            );
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors);
+        }
+
+        #[test]
+        fn sorted_insertion_preserves_existing_order_for_equal_distances() {
+            // Given
+            let nearest = LeafNeighbor::new(1, 1.0);
+            let existing_tie = LeafNeighbor::new(2, 2.0);
+            let farthest = LeafNeighbor::new(3, 4.0);
+            let tied_candidate = LeafNeighbor::new(4, 2.0);
+            let expected_neighbors = [nearest, existing_tie, tied_candidate];
+            let mut actual_neighbors = [nearest, existing_tie, farthest];
+
+            // When
+            insert_neighbor(
+                &mut actual_neighbors,
+                tied_candidate.target,
+                tied_candidate.distance,
+            );
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors);
+        }
+    }
+
+    mod leaf_neighbor_count_tests {
+        use super::leaf_neighbor_count;
+
+        #[test]
+        fn empty_leaf_cannot_retain_neighbors() {
+            // Given
+            let point_count = 0;
+            let requested_k = 3;
+            let expected_neighbor_count = 0;
+
+            // When
+            let actual_neighbor_count = leaf_neighbor_count(point_count, requested_k);
+
+            // Then
+            assert_eq!(actual_neighbor_count, expected_neighbor_count);
+        }
+
+        #[test]
+        fn singleton_leaf_cannot_retain_its_source_point() {
+            // Given
+            let point_count = 1;
+            let requested_k = 3;
+            let expected_non_self_neighbor_count = 0;
+
+            // When
+            let actual_neighbor_count = leaf_neighbor_count(point_count, requested_k);
+
+            // Then
+            assert_eq!(actual_neighbor_count, expected_non_self_neighbor_count);
+        }
+
+        #[test]
+        fn requested_k_above_available_neighbors_is_clamped() {
+            // Given
+            let point_count = 4;
+            let requested_k = 4;
+            let expected_all_non_self_neighbors = point_count - 1;
+
+            // When
+            let actual_neighbor_count = leaf_neighbor_count(point_count, requested_k);
+
+            // Then
+            assert_eq!(actual_neighbor_count, expected_all_non_self_neighbors);
+        }
+
+        #[test]
+        fn requested_k_within_available_neighbors_is_unchanged() {
+            // Given
+            let point_count = 8;
+            let requested_k = 5;
+            let expected_requested_neighbor_count = requested_k;
+
+            // When
+            let actual_neighbor_count = leaf_neighbor_count(point_count, requested_k);
+
+            // Then
+            assert_eq!(actual_neighbor_count, expected_requested_neighbor_count);
+        }
+    }
+
+    mod select_leaf_neighbors_tests {
+        use super::*;
+
+        #[test]
+        fn l2_pipeline_orders_neighbors_by_squared_distance() {
+            // Given
+            let values = [0.0_f32, 1.0, 3.0, 10.0];
+            let points = MatrixView::try_from(&values[..], 4, 1).unwrap();
+            let expected_neighbors = [
+                LeafNeighbor::new(1, (values[0] - values[1]).powi(2)),
+                LeafNeighbor::new(2, (values[0] - values[2]).powi(2)),
+                LeafNeighbor::new(0, (values[1] - values[0]).powi(2)),
+                LeafNeighbor::new(2, (values[1] - values[2]).powi(2)),
+                LeafNeighbor::new(1, (values[2] - values[1]).powi(2)),
+                LeafNeighbor::new(0, (values[2] - values[0]).powi(2)),
+                LeafNeighbor::new(2, (values[3] - values[2]).powi(2)),
+                LeafNeighbor::new(1, (values[3] - values[1]).powi(2)),
+            ];
+            let mut actual_neighbors = [LeafNeighbor::default(); 8];
+
+            // When
+            select_leaf_neighbors::<_, L2>(
+                diskann_wide::ARCH,
+                points,
+                MutMatrixView::try_from(&mut actual_neighbors[..], 4, 2).unwrap(),
+                &mut LeafKernelWorkspace::default(),
+            )
+            .unwrap();
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors);
+        }
+
+        #[test]
+        fn reused_workspace_matches_fresh_neighbor_selection() {
+            // Given
+            let values = [0.0_f32, 1.0, 3.0, 10.0];
+            let smaller_points = MatrixView::try_from(&values[..3], 3, 1).unwrap();
+            let mut reused_workspace = LeafKernelWorkspace::default();
+            let mut discarded_large_output = [LeafNeighbor::default(); 8];
+            select_leaf_neighbors::<_, L2>(
+                diskann_wide::ARCH,
+                MatrixView::try_from(&values[..], 4, 1).unwrap(),
+                MutMatrixView::try_from(&mut discarded_large_output[..], 4, 2).unwrap(),
+                &mut reused_workspace,
+            )
+            .unwrap();
+            let mut expected_neighbors_from_fresh_workspace = [LeafNeighbor::default(); 6];
+            select_leaf_neighbors::<_, L2>(
+                diskann_wide::ARCH,
+                smaller_points,
+                MutMatrixView::try_from(&mut expected_neighbors_from_fresh_workspace[..], 3, 2)
+                    .unwrap(),
+                &mut LeafKernelWorkspace::default(),
+            )
+            .unwrap();
+
+            // When
+            let mut actual_neighbors_from_reused_workspace = [LeafNeighbor::default(); 6];
+            select_leaf_neighbors::<_, L2>(
+                diskann_wide::ARCH,
+                smaller_points,
+                MutMatrixView::try_from(&mut actual_neighbors_from_reused_workspace[..], 3, 2)
+                    .unwrap(),
+                &mut reused_workspace,
+            )
+            .unwrap();
+
+            // Then
+            assert_eq!(
+                actual_neighbors_from_reused_workspace,
+                expected_neighbors_from_fresh_workspace
+            );
+        }
+    }
+
+    mod rank_leaf_dots_tests {
+        use super::*;
+        use rstest::rstest;
+
+        #[rstest]
+        #[case::two_points_fixed_one(2, 1)]
+        #[case::scalar_fixed_two(7, 2)]
+        #[case::lane_minus_one_fixed_three(15, 3)]
+        #[case::one_complete_lane_fixed_three(16, 3)]
+        #[case::lane_plus_one_runtime_width(17, 4)]
+        #[case::two_lanes_minus_one_runtime_width(31, 7)]
+        #[case::two_complete_lanes_runtime_width(32, 7)]
+        #[case::two_lanes_plus_one_runtime_width(33, 7)]
+        #[case::four_complete_lanes_runtime_width(64, 7)]
+        #[case::sixteen_complete_lanes_runtime_width(256, 7)]
+        #[case::maximum_leaf_size_runtime_width(512, 7)]
+        #[trace]
+        fn dispatched_leaf_ranking_matches_scalar_reference_across_lane_boundaries(
+            #[values(
+                Metric::L2,
+                Metric::Cosine,
+                Metric::CosineNormalized,
+                Metric::InnerProduct
+            )]
+            metric: Metric,
+            #[case] point_count: usize,
+            #[case] requested_k: usize,
+        ) {
+            // Given
+            let dots = index_distance_lower_gram(point_count);
+            let expected_neighbors = reference_neighbors(metric, &dots, point_count, requested_k);
+
+            // When
+            let actual_neighbors = rank_neighbors(metric, &dots, point_count, requested_k).1;
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors);
+        }
+
+        const UNIT_SQUARED_NORM: f32 = 1.0;
+        const POINT_0_1_DOT: f32 = 0.0;
+        const POINT_0_2_DOT: f32 = -1.0;
+        const POINT_1_2_DOT: f32 = 0.5;
+
+        #[rustfmt::skip]
+    const THREE_POINT_LOWER_GRAM: [f32; 9] = [
+        UNIT_SQUARED_NORM, f32::NAN,          f32::NAN,
+        POINT_0_1_DOT,      UNIT_SQUARED_NORM, f32::NAN,
+        POINT_0_2_DOT,      POINT_1_2_DOT,     UNIT_SQUARED_NORM,
+    ];
+
+        fn rank_three_point_fixture(metric: Metric) -> Vec<LeafNeighbor> {
+            rank_neighbors(metric, &THREE_POINT_LOWER_GRAM, 3, 1).1
+        }
+
+        #[test]
+        fn l2_selects_the_nearest_target_from_the_lower_gram_triangle() {
+            // Given
+            let expected_neighbors = [
+                LeafNeighbor::new(1, 2.0 * UNIT_SQUARED_NORM - 2.0 * POINT_0_1_DOT),
+                LeafNeighbor::new(2, 2.0 * UNIT_SQUARED_NORM - 2.0 * POINT_1_2_DOT),
+                LeafNeighbor::new(1, 2.0 * UNIT_SQUARED_NORM - 2.0 * POINT_1_2_DOT),
+            ];
+
+            // When
+            let actual_neighbors = rank_three_point_fixture(Metric::L2);
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors);
+        }
+
+        #[test]
+        fn cosine_selects_the_nearest_target_from_the_lower_gram_triangle() {
+            // Given
+            let expected_neighbors = [
+                LeafNeighbor::new(1, 1.0 - POINT_0_1_DOT),
+                LeafNeighbor::new(2, 1.0 - POINT_1_2_DOT),
+                LeafNeighbor::new(1, 1.0 - POINT_1_2_DOT),
+            ];
+
+            // When
+            let actual_neighbors = rank_three_point_fixture(Metric::Cosine);
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors);
+        }
+
+        #[test]
+        fn normalized_cosine_selects_the_largest_dot_product() {
+            // Given
+            let expected_neighbors = [
+                LeafNeighbor::new(1, 1.0 - POINT_0_1_DOT),
+                LeafNeighbor::new(2, 1.0 - POINT_1_2_DOT),
+                LeafNeighbor::new(1, 1.0 - POINT_1_2_DOT),
+            ];
+
+            // When
+            let actual_neighbors = rank_three_point_fixture(Metric::CosineNormalized);
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors);
+        }
+
+        #[test]
+        fn inner_product_selects_the_largest_dot_product() {
+            // Given
+            let expected_neighbors = [
+                LeafNeighbor::new(1, -POINT_0_1_DOT),
+                LeafNeighbor::new(2, -POINT_1_2_DOT),
+                LeafNeighbor::new(1, -POINT_1_2_DOT),
+            ];
+
+            // When
+            let actual_neighbors = rank_three_point_fixture(Metric::InnerProduct);
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors);
+        }
+
+        #[test]
+        fn equal_l2_distances_keep_target_scan_order() {
+            // Given
+            let unit_squared_norm = 1.0;
+            let tied_dot_product = 0.0;
+            let expected_tied_distance = 2.0 * unit_squared_norm - 2.0 * tied_dot_product;
+            #[rustfmt::skip]
+        let dots = [
+            unit_squared_norm, f32::NAN,          f32::NAN,          f32::NAN,
+            tied_dot_product,  unit_squared_norm, f32::NAN,          f32::NAN,
+            tied_dot_product,  tied_dot_product,  unit_squared_norm, f32::NAN,
+            tied_dot_product,  tied_dot_product,  tied_dot_product,  unit_squared_norm,
+        ];
+            let expected_neighbors_in_scan_order = [
+                LeafNeighbor::new(1, expected_tied_distance),
+                LeafNeighbor::new(2, expected_tied_distance),
+                LeafNeighbor::new(0, expected_tied_distance),
+                LeafNeighbor::new(2, expected_tied_distance),
+                LeafNeighbor::new(0, expected_tied_distance),
+                LeafNeighbor::new(1, expected_tied_distance),
+                LeafNeighbor::new(0, expected_tied_distance),
+                LeafNeighbor::new(1, expected_tied_distance),
+            ];
+
+            // When
+            let actual_neighbors = rank_neighbors(Metric::L2, &dots, 4, 2).1;
+
+            // Then
+            assert_eq!(actual_neighbors, expected_neighbors_in_scan_order);
+        }
+
+        #[test]
+        fn scalar_l2_distance_stays_finite_when_twice_the_dot_product_overflows() {
+            // Given
+            let dot_product = f32::from_bits(f32::MAX.to_bits() - 1);
+            let unfused_twice_dot_product = 2.0 * dot_product;
+            let expected_fused_distance = (-2.0_f32).mul_add(dot_product, f32::MAX) + f32::MAX;
+            let dots = [f32::MAX, 0.0, dot_product, f32::MAX];
+
+            // When
+            let actual_neighbors = rank_neighbors(Metric::L2, &dots, 2, 1).1;
+
+            // Then
+            assert!(unfused_twice_dot_product.is_infinite());
+            assert!(expected_fused_distance.is_finite() && expected_fused_distance > 0.0);
+            assert_eq!(
+                actual_neighbors[0].distance.to_bits(),
+                expected_fused_distance.to_bits()
+            );
+        }
+
+        #[test]
+        fn simd_l2_distance_stays_finite_when_twice_the_dot_product_overflows() {
+            // Given
+            let dot_product = f32::from_bits(f32::MAX.to_bits() - 1);
+            let unfused_twice_dot_product = 2.0 * dot_product;
+            let expected_fused_distance = (-2.0_f32).mul_add(dot_product, f32::MAX) + f32::MAX;
+            let expected_simd_neighbor_target = 0;
+            let points = 17;
+            let mut dots = square_matrix_with_constant_diagonal(points, f32::MAX);
+            dots[16 * points] = dot_product;
+
+            // When
+            let actual_neighbors = rank_neighbors(Metric::L2, &dots, points, 1).1;
+
+            // Then
+            assert!(unfused_twice_dot_product.is_infinite());
+            assert_eq!(actual_neighbors[16].target, expected_simd_neighbor_target);
+            assert_eq!(
+                actual_neighbors[16].distance.to_bits(),
+                expected_fused_distance.to_bits()
+            );
+        }
+
+        #[test]
+        fn cosine_zero_norm_produces_unit_distance() {
+            // Given
+            #[rustfmt::skip]
+        let dots = [
+            0.0, 99.0, 99.0,
+            0.0,  1.0, 99.0,
+            0.0,  0.0,  1.0,
+        ];
+            let expected_zero_norm_neighbors =
+                [LeafNeighbor::new(1, 1.0), LeafNeighbor::new(2, 1.0)];
+
+            // When
+            let actual_neighbors = rank_neighbors(Metric::Cosine, &dots, 3, 2).1;
+
+            // Then
+            assert_eq!(&actual_neighbors[..2], &expected_zero_norm_neighbors);
+        }
+
+        #[test]
+        fn cosine_similarity_above_one_clamps_to_zero_distance() {
+            // Given
+            let dots = [1.0, 0.0, 2.0, 1.0];
+            let maximum_cosine_similarity = 1.0;
+            let expected_one_minus_maximum_similarity = 1.0 - maximum_cosine_similarity;
+
+            // When
+            let actual_neighbors = rank_neighbors(Metric::Cosine, &dots, 2, 1).1;
+
+            // Then
+            assert_eq!(
+                actual_neighbors[0].distance,
+                expected_one_minus_maximum_similarity
+            );
+        }
+
+        #[test]
+        fn cosine_similarity_below_negative_one_clamps_to_distance_two() {
+            // Given
+            let dots = [1.0, 0.0, -2.0, 1.0];
+            let minimum_cosine_similarity = -1.0;
+            let expected_one_minus_minimum_similarity = 1.0 - minimum_cosine_similarity;
+
+            // When
+            let actual_neighbors = rank_neighbors(Metric::Cosine, &dots, 2, 1).1;
+
+            // Then
+            assert_eq!(
+                actual_neighbors[0].distance,
+                expected_one_minus_minimum_similarity
+            );
+        }
+
+        #[test]
+        fn cosine_subnormal_norm_is_treated_as_zero() {
+            // Given
+            let dots = [f32::MIN_POSITIVE / 2.0, 0.0, 1.0, 1.0];
+            let zero_norm_similarity = 0.0;
+            let expected_one_minus_zero_similarity = 1.0 - zero_norm_similarity;
+
+            // When
+            let actual_neighbors = rank_neighbors(Metric::Cosine, &dots, 2, 1).1;
+
+            // Then
+            assert_eq!(
+                actual_neighbors[0].distance,
+                expected_one_minus_zero_similarity
+            );
+        }
+
+        #[test]
+        fn f32_max_distance_is_still_a_rankable_neighbor() {
+            // Given
+            let points = 4;
+            let expected_leaf_k = 3;
+            let expected_last_neighbor = LeafNeighbor::new(0, f32::MAX);
+            let mut dots = vec![0.0; points * points];
+            dots[3 * points] = -f32::MAX;
+
+            // When
+            let (actual_leaf_k, actual_neighbors) =
+                rank_neighbors(Metric::InnerProduct, &dots, points, expected_leaf_k);
+
+            // Then
+            assert_eq!(actual_leaf_k, expected_leaf_k);
+            assert_eq!(
+                actual_neighbors[3 * actual_leaf_k + actual_leaf_k - 1],
+                expected_last_neighbor
+            );
+        }
+
+        #[test]
+        fn scalar_nan_distance_leaves_the_neighbor_slot_unassigned() {
+            // Given
+            let dots = [1.0, 0.0, f32::NAN, 1.0];
+            let expected_unassigned_neighbors = [LeafNeighbor::default(), LeafNeighbor::default()];
+
+            // When
+            let actual_neighbors = rank_neighbors(Metric::CosineNormalized, &dots, 2, 1).1;
+
+            // Then
+            assert_eq!(actual_neighbors, expected_unassigned_neighbors);
+        }
+
+        #[test]
+        fn simd_nan_distance_cannot_replace_a_finite_neighbor() {
+            // Given
+            let points = 17;
+            let mut dots = square_matrix_with_constant_diagonal(points, 1.0);
+            dots[16 * points] = f32::NAN;
+            let expected_finite_neighbor = LeafNeighbor::new(1, 1.0);
+
+            // When
+            let actual_neighbors = rank_neighbors(Metric::CosineNormalized, &dots, points, 1).1;
+
+            // Then
+            assert_eq!(actual_neighbors[16], expected_finite_neighbor);
+        }
+
+        #[test]
+        fn empty_leaf_has_no_neighbors() {
+            // Given
+            let dots = [];
+            let expected_zero_neighbor_width = 0;
+            let expected_no_neighbors: [LeafNeighbor; 0] = [];
+
+            // When
+            let (actual_leaf_k, actual_neighbors) = rank_neighbors(Metric::L2, &dots, 0, 2);
+
+            // Then
+            assert_eq!(actual_leaf_k, expected_zero_neighbor_width);
+            assert_eq!(actual_neighbors, expected_no_neighbors);
+        }
+
+        #[test]
+        fn singleton_leaf_has_no_neighbors() {
+            // Given
+            let dots = [4.0];
+            let expected_zero_neighbor_width = 0;
+            let expected_no_neighbors: [LeafNeighbor; 0] = [];
+
+            // When
+            let (actual_leaf_k, actual_neighbors) = rank_neighbors(Metric::Cosine, &dots, 1, 2);
+
+            // Then
+            assert_eq!(actual_leaf_k, expected_zero_neighbor_width);
+            assert_eq!(actual_neighbors, expected_no_neighbors);
+        }
+
+        #[test]
+        fn zero_requested_k_has_no_neighbors() {
+            // Given
+            let dots = [1.0, 0.0, 0.0, 1.0];
+            let expected_zero_neighbor_width = 0;
+            let expected_no_neighbors: [LeafNeighbor; 0] = [];
+
+            // When
+            let (actual_leaf_k, actual_neighbors) =
+                rank_neighbors(Metric::InnerProduct, &dots, 2, 0);
+
+            // Then
+            assert_eq!(actual_leaf_k, expected_zero_neighbor_width);
+            assert_eq!(actual_neighbors, expected_no_neighbors);
+        }
+
+        #[test]
+        fn neighbor_width_equal_to_point_count_is_rejected() {
+            // Given
+            let dots = [0.0; 9];
+            let input = lower_gram_view(&dots, 3);
+            let norms = metric_norms(Metric::L2, input);
+            let expected_error = LeafKernelError::InvalidNeighborCount {
+                points: 3,
+                neighbors: 3,
+                maximum: 2,
+            };
+            let mut output = [LeafNeighbor::default(); 9];
+
+            // When
+            let actual_error = arch::dispatch1_no_features(
+                DispatchMetric(Metric::L2),
+                KernelCall {
+                    input,
+                    norms: &norms,
+                    output: MutMatrixView::try_from(&mut output[..], 3, 3).unwrap(),
+                    workspace: &mut LeafKernelWorkspace::default(),
+                },
+            )
+            .unwrap_err();
+
+            // Then
+            assert_eq!(actual_error, expected_error);
+        }
+    }
+}
diff --git a/diskann/src/graph/pipnn/mod.rs b/diskann/src/graph/pipnn/mod.rs
new file mode 100644
index 0000000000..092ee7549a
--- /dev/null
+++ b/diskann/src/graph/pipnn/mod.rs
@@ -0,0 +1,33 @@
+/*
+ * Copyright (c) Microsoft Corporation.
+ * Licensed under the MIT license.
+ */
+
+//! Numerical kernels for PiPNN graph construction.
+//!
+//! [`partition_kernel`] converts point-to-leader dot products into sorted leader
+//! positions. The output width sets the fanout. A scratch vector stores the
+//! ranked leaders and reuses its allocation for each point.
+//!
+//! [`leaf_kernel`] reads a lower-triangular Gram matrix. It evaluates each point
+//! pair once and updates both points. Each point retains at most three local
+//! neighbors.
+//!
+//! `kernel_metric` defines metric markers and shared math. Separate leaf and
+//! partition traits define norm preparation and ranking formulas.
+//!
+//! The graph builder selects architecture `A` and metric `M` once. It passes
+//! these concrete types to both kernels.
+//!
+//! Each kernel checks all view and norm relationships before unchecked SIMD
+//! access. The kernels borrow their matrices. They write only to caller-owned
+//! output and workspace.
+#[allow(dead_code)]
+mod kernel_metric;
+#[allow(dead_code)]
+mod simd;
+
+#[allow(dead_code)]
+mod leaf_kernel;
+#[allow(dead_code)]
+mod partition_kernel;
diff --git a/diskann/src/graph/pipnn/partition_kernel.rs b/diskann/src/graph/pipnn/partition_kernel.rs
new file mode 100644
index 0000000000..ba03067cb6
--- /dev/null
+++ b/diskann/src/graph/pipnn/partition_kernel.rs
@@ -0,0 +1,652 @@
+/*
+ * Copyright (c) Microsoft Corporation.
+ * Licensed under the MIT license.
+ */
+
+//! Select partition centers for PiPNN point assignment.
+//!
+//! A leader is a sampled dataset point that represents one child partition.
+//! The kernel prepares reusable leader norms, computes point-to-leader dot
+//! products, and returns nearest leader-column IDs for partition scatter.
+//!
+//! L2 omits the assigned point's norm because it is constant across all sampled
+//! leaders. Equal scores keep sampled-leader order. NaN is not rankable. An
+//! unfilled output slot contains [`UNASSIGNED_LEADER`].
+
+use std::marker::PhantomData;
+
+use crate::{ANNError, ANNResult};
+use diskann_linalg::Transpose;
+use diskann_utils::views::{MatrixView, MutMatrixView};
+use diskann_wide::{SIMDMask, SIMDVector};
+
+use super::{
+    kernel_metric::{PartitionMetric, PartitionNorms},
+    simd::{PiPNNSIMDSchema, PiPNNSIMDVector},
+};
+
+/// No sampled partition center was rankable for this output slot.
+pub(super) const UNASSIGNED_LEADER: u32 = u32::MAX;
+
+/// Sampled leader vectors with metric-specific reusable norms.
+pub(super) struct PreparedLeaders<'a, M> {
+    leader_values: MatrixView<'a, f32>,
+    leader_norms: Vec<f32>,
+    metric: PhantomData<M>,
+}
+
+impl<'a, M> PreparedLeaders<'a, M>
+where
+    M: PartitionMetric,
+{
+    /// Prepare leader state for all point stripes in one partition split.
+    pub(super) fn new(leader_values: MatrixView<'a, f32>) -> Self {
+        let mut leader_norms = Vec::new();
+        M::prepare_leader_norms(leader_values, &mut leader_norms);
+        Self {
+            leader_values,
+            leader_norms,
+            metric: PhantomData,
+        }
+    }
+
+    pub(super) fn len(&self) -> usize {
+        self.leader_values.nrows()
+    }
+}
+
+/// Reusable storage for one point-stripe numerical pipeline.
+#[derive(Default)]
+pub(super) struct PartitionKernelWorkspace {
+    dot_scratch: Vec<f32>,
+    point_norm_scratch: Vec<f32>,
+    ranked_leader_scratch: Vec<(u32, f32)>,
+}
+
+/// Dot products between assigned points and sampled partition centers.
+///
+/// Each row is one point being assigned. Each column is one sampled leader.
+/// [`Self::norms`] supplies the norm layout for metric `M`.
+#[derive(Clone, Copy, Debug)]
+struct PartitionInput<'a> {
+    dots: MatrixView<'a, f32>,
+    norms: PartitionNorms<'a>,
+}
+
+/// Assign one packed point stripe to prepared partition leaders.
+///
+/// A point can have fewer assignments than the output width. Each remaining
+/// slot contains [`UNASSIGNED_LEADER`].
+///
+/// # Errors
+///
+/// Returns an error for invalid GEMM input.
+pub(super) fn assign_leaders<A, M>(
+    arch: A,
+    points: MatrixView<'_, f32>,
+    leaders: &PreparedLeaders<'_, M>,
+    output: MutMatrixView<'_, u32>,
+    workspace: &mut PartitionKernelWorkspace,
+) -> ANNResult<()>
+where
+    A: PiPNNSIMDSchema,
+    M: PartitionMetric,
+{
+    let point_count = points.nrows();
+    let leader_count = leaders.len();
+    let dot_count = point_count * leader_count;
+    let PartitionKernelWorkspace {
+        dot_scratch,
+        point_norm_scratch,
+        ranked_leader_scratch,
+    } = workspace;
+    if dot_scratch.len() < dot_count {
+        dot_scratch.resize(dot_count, 0.0);
+    }
+    diskann_linalg::sgemm(
+        Transpose::None,
+        Transpose::Ordinary,
+        point_count,
+        leader_count,
+        points.ncols(),
+        1.0,
+        points.as_slice(),
+        leaders.leader_values.as_slice(),
+        None,
+        &mut dot_scratch[..dot_count],
+    )
+    .map_err(ANNError::new)?;
+    M::prepare_point_norms(points, point_norm_scratch);
+    let dots = MatrixView::try_from(&dot_scratch[..dot_count], point_count, leader_count)
+        .map_err(|error| ANNError::new(error.as_static()))?;
+    rank_leader_dots::<A, M>(
+        arch,
+        PartitionInput {
+            dots,
+            norms: PartitionNorms {
+                point_norms: point_norm_scratch,
+                leader_norms: &leaders.leader_norms,
+            },
+        },
+        output,
+        ranked_leader_scratch,
+    );
+    Ok(())
+}
+
+/// Rank prepared point-to-leader dot products.
+fn rank_leader_dots<A, M>(
+    arch: A,
+    input: PartitionInput<'_>,
+    output: MutMatrixView<'_, u32>,
+    ranked_leaders: &mut Vec<(u32, f32)>,
+) where
+    A: PiPNNSIMDSchema,
+    M: PartitionMetric,
+{
+    let fanout = output.ncols();
+    if fanout == 0 || input.dots.nrows() == 0 {
+        return;
+    }
+
+    ranked_leaders.resize(fanout, (UNASSIGNED_LEADER, f32::INFINITY));
+    select_point_leaders::<A, M>(arch, input.dots, input.norms, output, ranked_leaders);
+}
+
+/// Rank sampled partition centers for each assigned point.
+///
+/// The function keeps nearest-first order for every point. Full SIMD groups use
+/// metric-specific formulas. Remaining leaders use the matching single formula.
+fn select_point_leaders<A, M>(
+    arch: A,
+    dots: MatrixView<'_, f32>,
+    norms: PartitionNorms<'_>,
+    mut output: MutMatrixView<'_, u32>,
+    ranked_leaders: &mut [(u32, f32)],
+) where
+    A: PiPNNSIMDSchema,
+    M: PartitionMetric,
+{
+    let leader_count = dots.ncols();
+    let fanout = output.ncols();
+
+    for (point, (point_dots, point_output)) in dots
+        .row_iter()
+        .zip(output.as_mut_slice().chunks_exact_mut(fanout))
+        .enumerate()
+    {
+        ranked_leaders.fill((UNASSIGNED_LEADER, f32::INFINITY));
+        let point_simd = M::point_simd(arch, norms, point);
+        let point_single = M::point_single(norms, point);
+        let simd_prefix = leader_count - leader_count % M::Simd::<A>::LANES;
+
+        for first_leader in (0..simd_prefix).step_by(M::Simd::<A>::LANES) {
+            // SAFETY: This group is inside the point's leader row.
+            let dot_products =
+                unsafe { M::Simd::<A>::load_simd(arch, point_dots.as_ptr().add(first_leader)) };
+            let rankings = M::rankings_simd(arch, norms, point_simd, dot_products, first_leader);
+            insert_leader_lanes(rankings, first_leader, ranked_leaders);
+        }
+
+        for (leader, &dot_product) in point_dots.iter().enumerate().skip(simd_prefix) {
+            let ranking = M::ranking_single(norms, point_single, dot_product, leader);
+            insert_leader(ranked_leaders, leader as u32, ranking);
+        }
+        for (destination, &(leader, _)) in point_output.iter_mut().zip(ranked_leaders.iter()) {
+            *destination = leader;
+        }
+    }
+}
+
+/// Offer one SIMD group of sampled centers to the current point's ranked_leaders.
+///
+/// `first_leader` is the matrix-column ID of the first lane. Lanes enter in
+/// sampled-leader order, which preserves tie order.
+fn insert_leader_lanes<F>(scores: F, first_leader: usize, ranked_leaders: &mut [(u32, f32)])
+where
+    F: PiPNNSIMDVector,
+{
+    let threshold = F::splat(scores.arch(), ranked_leaders[ranked_leaders.len() - 1].1);
+    let eligible = scores.lt_simd(threshold);
+    if eligible.none() {
+        return;
+    }
+
+    let values = scores.to_array();
+    let values = values.as_ref();
+    let mut lanes = F::active_lanes(eligible);
+    while lanes != 0 {
+        let lane = lanes.trailing_zeros() as usize;
+        lanes &= lanes - 1;
+        insert_leader(ranked_leaders, (first_leader + lane) as u32, values[lane]);
+    }
+}
+
+/// Insert one sampled partition center into the current point's retained set.
+///
+/// `leader` is the center's column ID in the point-to-leader matrix. `ranked_leaders`
+/// stores retained centers in nearest-first order. Equal scores and NaN do not
+/// enter, so sampled-leader order resolves ties.
+#[inline(always)]
+fn insert_leader(ranked_leaders: &mut [(u32, f32)], leader: u32, score: f32) {
+    let threshold = ranked_leaders.len() - 1;
+    if score.partial_cmp(&ranked_leaders[threshold].1) != Some(std::cmp::Ordering::Less) {
+        return;
+    }
+
+    ranked_leaders[threshold] = (leader, score);
+    let mut slot = threshold;
+    while slot > 0 && ranked_leaders[slot].1 < ranked_leaders[slot - 1].1 {
+        ranked_leaders.swap(slot, slot - 1);
+        slot -= 1;
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::graph::pipnn::kernel_metric::{Cosine, CosineNormalized, InnerProduct, L2};
+    use diskann_utils::views::{Matrix, MatrixView, MutMatrixView};
+    use diskann_vector::distance::Metric;
+    use diskann_wide::arch::{self, Target1};
+
+    struct KernelCall<'a> {
+        input: PartitionInput<'a>,
+        output: MutMatrixView<'a, u32>,
+        ranked_leaders: &'a mut Vec<(u32, f32)>,
+    }
+
+    struct DispatchMetric(Metric);
+
+    impl<A> Target1<A, (), KernelCall<'_>> for DispatchMetric
+    where
+        A: PiPNNSIMDSchema,
+    {
+        fn run(self, arch: A, call: KernelCall<'_>) {
+            match self.0 {
+                Metric::L2 => {
+                    rank_leader_dots::<A, L2>(arch, call.input, call.output, call.ranked_leaders)
+                }
+                Metric::Cosine => rank_leader_dots::<A, Cosine>(
+                    arch,
+                    call.input,
+                    call.output,
+                    call.ranked_leaders,
+                ),
+                Metric::CosineNormalized => rank_leader_dots::<A, CosineNormalized>(
+                    arch,
+                    call.input,
+                    call.output,
+                    call.ranked_leaders,
+                ),
+                Metric::InnerProduct => rank_leader_dots::<A, InnerProduct>(
+                    arch,
+                    call.input,
+                    call.output,
+                    call.ranked_leaders,
+                ),
+            }
+        }
+    }
+
+    fn partition_input<'a>(
+        dots: &'a [f32],
+        point_count: usize,
+        leader_count: usize,
+        point_norms: &'a [f32],
+        leader_norms: &'a [f32],
+    ) -> PartitionInput<'a> {
+        PartitionInput {
+            dots: MatrixView::try_from(dots, point_count, leader_count).unwrap(),
+            norms: PartitionNorms {
+                point_norms,
+                leader_norms,
+            },
+        }
+    }
+
+    fn rank_partition_leaders(
+        metric: Metric,
+        input: PartitionInput<'_>,
+        fanout: usize,
+    ) -> Vec<u32> {
+        let mut output = Matrix::new(u32::MAX, input.dots.nrows(), fanout);
+        arch::dispatch1_no_features(
+            DispatchMetric(metric),
+            KernelCall {
+                input,
+                output: output.as_mut_view(),
+                ranked_leaders: &mut Vec::new(),
+            },
+        );
+        output.into_inner().into_vec()
+    }
+
+    fn reference_score(metric: Metric, dot: f32, point_norm: f32, leader_norm: f32) -> f32 {
+        match metric {
+            Metric::L2 => (-2.0_f32).mul_add(dot, leader_norm),
+            Metric::CosineNormalized => 1.0 - dot,
+            Metric::InnerProduct => -dot,
+            Metric::Cosine => {
+                if point_norm < f32::MIN_POSITIVE.sqrt() || leader_norm < f32::MIN_POSITIVE.sqrt() {
+                    1.0
+                } else {
+                    1.0 - (dot / (point_norm * leader_norm)).clamp(-1.0, 1.0)
+                }
+            }
+        }
+    }
+
+    fn reference_assignments(metric: Metric, input: PartitionInput<'_>, fanout: usize) -> Vec<u32> {
+        let mut output = vec![UNASSIGNED_LEADER; input.dots.nrows() * fanout];
+        for (point, (dots, assignments)) in input
+            .dots
+            .row_iter()
+            .zip(output.chunks_exact_mut(fanout))
+            .enumerate()
+        {
+            let point_norm = input.norms.point_norms.get(point).copied().unwrap_or(0.0);
+            let mut candidates: Vec<_> = dots
+                .iter()
+                .enumerate()
+                .filter_map(|(leader, &dot)| {
+                    let leader_norm = input.norms.leader_norms.get(leader).copied().unwrap_or(0.0);
+                    let score = reference_score(metric, dot, point_norm, leader_norm);
+                    (score.partial_cmp(&f32::INFINITY) == Some(std::cmp::Ordering::Less))
+                        .then_some((leader as u32, score))
+                })
+                .collect();
+            candidates.sort_by(|left, right| left.1.total_cmp(&right.1));
+            for (destination, (leader, _)) in assignments.iter_mut().zip(candidates) {
+                *destination = leader;
+            }
+        }
+        output
+    }
+
+    /// Build two score rows whose preferred leader direction is opposite.
+    /// The fractional scores cross complete SIMD groups and scalar tails.
+    fn lane_boundary_fixture(
+        metric: Metric,
+        leader_count: usize,
+    ) -> (Vec<f32>, Vec<f32>, Vec<f32>) {
+        let mut dots = Vec::with_capacity(2 * leader_count);
+        for point in 0..2 {
+            for leader in 0..leader_count {
+                let fraction = leader as f32 / leader_count as f32;
+                dots.push(if point == 0 { fraction } else { 1.0 - fraction });
+            }
+        }
+        let point_norms = if metric == Metric::Cosine {
+            vec![1.0; 2]
+        } else {
+            Vec::new()
+        };
+        let leader_norms = match metric {
+            Metric::L2 | Metric::Cosine => vec![1.0; leader_count],
+            Metric::CosineNormalized | Metric::InnerProduct => Vec::new(),
+        };
+        (dots, point_norms, leader_norms)
+    }
+
+    mod insert_leader_tests {
+        use super::*;
+
+        #[test]
+        fn topk_keeps_nearest_first_order_and_scan_order_ties() {
+            // Given
+            let expected_ranked_leaders = [(1, 1.0), (4, 1.0), (3, 2.0), (2, 3.0)];
+            let mut ranked_leaders = vec![(UNASSIGNED_LEADER, f32::INFINITY); 4];
+
+            // When
+            insert_leader(&mut ranked_leaders, 0, 4.0);
+            insert_leader(&mut ranked_leaders, 1, 1.0);
+            insert_leader(&mut ranked_leaders, 2, 3.0);
+            insert_leader(&mut ranked_leaders, 3, 2.0);
+            insert_leader(&mut ranked_leaders, 4, 1.0);
+
+            // Then
+            assert_eq!(ranked_leaders, expected_ranked_leaders);
+        }
+
+        #[test]
+        fn nan_score_does_not_enter_the_topk() {
+            // Given
+            let expected_ranked_leaders = [(0, 0.25), (UNASSIGNED_LEADER, f32::INFINITY)];
+            let mut ranked_leaders = vec![(UNASSIGNED_LEADER, f32::INFINITY); 2];
+
+            // When
+            insert_leader(&mut ranked_leaders, 0, 0.25);
+            insert_leader(&mut ranked_leaders, 1, f32::NAN);
+
+            // Then
+            assert_eq!(ranked_leaders, expected_ranked_leaders);
+        }
+    }
+
+    mod assign_leaders_tests {
+        use super::*;
+
+        #[test]
+        fn cosine_pipeline_assigns_each_point_to_its_nearest_leaders() {
+            // Given
+            let leader_values = [1.0, 0.0, 0.0, 1.0, -1.0, 0.0];
+            let leaders = PreparedLeaders::<Cosine>::new(
+                MatrixView::try_from(&leader_values[..], 3, 2).unwrap(),
+            );
+            let point_values = [0.9, 0.1, -0.8, 0.2];
+            let points = MatrixView::try_from(&point_values[..], 2, 2).unwrap();
+            let expected_leaders_by_descending_cosine_similarity = [0, 1, 2, 1];
+            let mut actual_assignments = [UNASSIGNED_LEADER; 4];
+
+            // When
+            assign_leaders::<_, Cosine>(
+                diskann_wide::ARCH,
+                points,
+                &leaders,
+                MutMatrixView::try_from(&mut actual_assignments[..], 2, 2).unwrap(),
+                &mut PartitionKernelWorkspace::default(),
+            )
+            .unwrap();
+
+            // Then
+            assert_eq!(
+                actual_assignments,
+                expected_leaders_by_descending_cosine_similarity
+            );
+        }
+
+        #[test]
+        fn reused_workspace_matches_fresh_leader_assignment() {
+            // Given
+            let leader_values = [1.0, 0.0, 0.0, 1.0, -1.0, 0.0];
+            let leaders = PreparedLeaders::<Cosine>::new(
+                MatrixView::try_from(&leader_values[..], 3, 2).unwrap(),
+            );
+            let point_values = [0.9, 0.1, -0.8, 0.2];
+            let smaller_points = MatrixView::try_from(&point_values[..2], 1, 2).unwrap();
+            let mut reused_workspace = PartitionKernelWorkspace::default();
+            let mut discarded_large_output = [UNASSIGNED_LEADER; 4];
+            assign_leaders::<_, Cosine>(
+                diskann_wide::ARCH,
+                MatrixView::try_from(&point_values[..], 2, 2).unwrap(),
+                &leaders,
+                MutMatrixView::try_from(&mut discarded_large_output[..], 2, 2).unwrap(),
+                &mut reused_workspace,
+            )
+            .unwrap();
+            let mut expected_assignments_from_fresh_workspace = [UNASSIGNED_LEADER; 2];
+            assign_leaders::<_, Cosine>(
+                diskann_wide::ARCH,
+                smaller_points,
+                &leaders,
+                MutMatrixView::try_from(&mut expected_assignments_from_fresh_workspace[..], 1, 2)
+                    .unwrap(),
+                &mut PartitionKernelWorkspace::default(),
+            )
+            .unwrap();
+
+            // When
+            let mut actual_assignments_from_reused_workspace = [UNASSIGNED_LEADER; 2];
+            assign_leaders::<_, Cosine>(
+                diskann_wide::ARCH,
+                smaller_points,
+                &leaders,
+                MutMatrixView::try_from(&mut actual_assignments_from_reused_workspace[..], 1, 2)
+                    .unwrap(),
+                &mut reused_workspace,
+            )
+            .unwrap();
+
+            // Then
+            assert_eq!(
+                actual_assignments_from_reused_workspace,
+                expected_assignments_from_fresh_workspace
+            );
+        }
+    }
+
+    mod rank_leader_dots_tests {
+        use super::*;
+        use rstest::rstest;
+
+        #[rstest]
+        #[case::two_leaders_fanout_one(2, 1)]
+        #[case::scalar_fanout_two(7, 2)]
+        #[case::lane_minus_one(15, 3)]
+        #[case::one_complete_lane(16, 3)]
+        #[case::lane_plus_one(17, 4)]
+        #[case::two_lanes_minus_one(31, 7)]
+        #[case::two_complete_lanes(32, 7)]
+        #[case::two_lanes_plus_one(33, 7)]
+        #[trace]
+        fn dispatched_partition_ranking_matches_scalar_reference_across_lane_boundaries(
+            #[values(
+                Metric::L2,
+                Metric::Cosine,
+                Metric::CosineNormalized,
+                Metric::InnerProduct
+            )]
+            metric: Metric,
+            #[case] leader_count: usize,
+            #[case] fanout: usize,
+        ) {
+            // Given
+            let (dots, point_norms, leader_norms) = lane_boundary_fixture(metric, leader_count);
+            let input = partition_input(&dots, 2, leader_count, &point_norms, &leader_norms);
+            let expected_assignments = reference_assignments(metric, input, fanout);
+
+            // When
+            let actual_assignments = rank_partition_leaders(metric, input, fanout);
+
+            // Then
+            assert_eq!(actual_assignments, expected_assignments);
+        }
+
+        #[test]
+        fn equal_l2_scores_keep_sampled_leader_order() {
+            // Given
+            let dots = [0.0, 0.0, 0.0, 0.0];
+            let leader_squared_norms = [1.0, 1.0, 1.0, 1.0];
+            let expected_sampled_leader_order = [0, 1];
+
+            // When
+            let actual_assignments = rank_partition_leaders(
+                Metric::L2,
+                partition_input(&dots, 1, 4, &[], &leader_squared_norms),
+                2,
+            );
+
+            // Then
+            assert_eq!(actual_assignments, expected_sampled_leader_order);
+        }
+
+        #[test]
+        fn cosine_zero_norm_keeps_sampled_leader_order() {
+            // Given
+            let dots = [100.0, -100.0];
+            let point_norms = [0.0];
+            let leader_norms = [1.0, 1.0];
+            let expected_sampled_leader_order = [0, 1];
+
+            // When
+            let actual_assignments = rank_partition_leaders(
+                Metric::Cosine,
+                partition_input(&dots, 1, 2, &point_norms, &leader_norms),
+                2,
+            );
+
+            // Then
+            assert_eq!(actual_assignments, expected_sampled_leader_order);
+        }
+
+        #[test]
+        fn f32_max_score_is_still_a_rankable_leader() {
+            // Given
+            let mut dots = [0.0; 8];
+            dots[7] = -f32::MAX;
+            let expected_all_leaders_in_scan_order = [0, 1, 2, 3, 4, 5, 6, 7];
+
+            // When
+            let actual_assignments = rank_partition_leaders(
+                Metric::InnerProduct,
+                partition_input(&dots, 1, 8, &[], &[]),
+                8,
+            );
+
+            // Then
+            assert_eq!(actual_assignments, expected_all_leaders_in_scan_order);
+        }
+
+        #[test]
+        fn nan_leader_does_not_displace_finite_leaders() {
+            // Given
+            let dots = [f32::NAN, 3.0, 2.0];
+            let expected_finite_leaders = [1, 2];
+
+            // When
+            let actual_assignments = rank_partition_leaders(
+                Metric::InnerProduct,
+                partition_input(&dots, 1, 3, &[], &[]),
+                2,
+            );
+
+            // Then
+            assert_eq!(actual_assignments, expected_finite_leaders);
+        }
+
+        #[test]
+        fn empty_point_matrix_produces_no_assignments() {
+            // Given
+            let dots = [];
+            let expected_no_assignments: [u32; 0] = [];
+
+            // When
+            let actual_assignments = rank_partition_leaders(
+                Metric::InnerProduct,
+                partition_input(&dots, 0, 3, &[], &[]),
+                2,
+            );
+
+            // Then
+            assert_eq!(actual_assignments, expected_no_assignments);
+        }
+
+        #[test]
+        fn zero_fanout_produces_no_assignments() {
+            // Given
+            let dots = [1.0, 2.0, 3.0];
+            let expected_no_assignments: [u32; 0] = [];
+
+            // When
+            let actual_assignments = rank_partition_leaders(
+                Metric::InnerProduct,
+                partition_input(&dots, 1, 3, &[], &[]),
+                0,
+            );
+
+            // Then
+            assert_eq!(actual_assignments, expected_no_assignments);
+        }
+    }
+}
diff --git a/diskann/src/graph/pipnn/simd.rs b/diskann/src/graph/pipnn/simd.rs
new file mode 100644
index 0000000000..42f6efe057
--- /dev/null
+++ b/diskann/src/graph/pipnn/simd.rs
@@ -0,0 +1,61 @@
+/*
+ * Copyright (c) Microsoft Corporation.
+ * Licensed under the MIT license.
+ */
+
+//! SIMD schema for PiPNN numerical kernels.
+
+use diskann_wide::{Architecture, SIMDFloat, SIMDMask, SIMDSelect, SIMDVector};
+
+/// Default SIMD representation used by every PiPNN numerical stage.
+///
+/// This alias is the single build-time width selection.
+type DefaultVector<A> = <A as Architecture>::f32x16;
+
+/// Operations required by PiPNN SIMD vectors.
+pub(super) trait PiPNNSIMDVector:
+    SIMDVector<Scalar = f32> + SIMDFloat + std::ops::Div<Output = Self>
+{
+    /// Return one bit for each selected lane.
+    fn active_lanes(mask: Self::Mask) -> u64;
+
+    /// Select one value from each pair of lanes.
+    fn select(mask: Self::Mask, if_true: Self, if_false: Self) -> Self;
+}
+
+impl<F> PiPNNSIMDVector for F
+where
+    F: SIMDVector<Scalar = f32> + SIMDFloat + std::ops::Div<Output = F>,
+    F::Mask: SIMDSelect<F>,
+    u64: From<<<F::Mask as SIMDMask>::BitMask as SIMDMask>::Underlying>,
+{
+    #[inline(always)]
+    fn active_lanes(mask: Self::Mask) -> u64 {
+        u64::from(mask.bitmask().to_underlying())
+    }
+
+    #[inline(always)]
+    fn select(mask: Self::Mask, if_true: Self, if_false: Self) -> Self {
+        mask.select(if_true, if_false)
+    }
+}
+
+/// Stage-specific SIMD representations for one architecture.
+pub(super) trait PiPNNSIMDSchema: Architecture {
+    /// SIMD representation for leaf distance scores.
+    type LeafScore: PiPNNSIMDVector<Arch = Self>;
+    /// SIMD representation for partition ranking scores.
+    type PartitionScore: PiPNNSIMDVector<Arch = Self>;
+    /// SIMD representation for relative-hash sketch comparisons.
+    type HashScore: PiPNNSIMDVector<Arch = Self>;
+}
+
+impl<A> PiPNNSIMDSchema for A
+where
+    A: Architecture,
+    DefaultVector<A>: PiPNNSIMDVector<Arch = A>,
+{
+    type LeafScore = DefaultVector<A>;
+    type PartitionScore = DefaultVector<A>;
+    type HashScore = DefaultVector<A>;
+}