diff --git a/crates/clawhdf5-accel/src/avx2.rs b/crates/clawhdf5-accel/src/avx2.rs index dcdbf7e..f2da7e3 100644 --- a/crates/clawhdf5-accel/src/avx2.rs +++ b/crates/clawhdf5-accel/src/avx2.rs @@ -25,6 +25,55 @@ unsafe fn hsum_256(v: __m256) -> f32 { _mm_cvtss_f32(result) } +/// AVX2 dot product of two `i8` slices, widened to `i32`. +/// +/// Each 16-byte half is sign-extended to sixteen `i16` lanes and multiplied +/// pairwise with `madd_epi16`, which sums adjacent products straight into +/// eight `i32` lanes — the widening that an autovectorised scalar loop does +/// in several shuffles is one instruction here. A pair sum is at most +/// `2 * 127 * 127`, far inside `i32`. +/// +/// # Safety +/// Caller must verify is_x86_feature_detected!("avx2"). +// SAFETY: Caller must have verified AVX2 via is_x86_feature_detected!. +#[target_feature(enable = "avx2")] +pub unsafe fn dot_i8(a: &[i8], b: &[i8]) -> i32 { + // SAFETY: Caller guarantees AVX2 is available per the # Safety contract; + // every load reads 32 bytes at an index checked against `len` first. + unsafe { + assert_eq!(a.len(), b.len()); + let len = a.len(); + let mut i = 0; + let mut acc0 = _mm256_setzero_si256(); + let mut acc1 = _mm256_setzero_si256(); + + while i + 32 <= len { + let va = _mm256_loadu_si256(a.as_ptr().add(i).cast()); + let vb = _mm256_loadu_si256(b.as_ptr().add(i).cast()); + let a_lo = _mm256_cvtepi8_epi16(_mm256_castsi256_si128(va)); + let b_lo = _mm256_cvtepi8_epi16(_mm256_castsi256_si128(vb)); + let a_hi = _mm256_cvtepi8_epi16(_mm256_extracti128_si256(va, 1)); + let b_hi = _mm256_cvtepi8_epi16(_mm256_extracti128_si256(vb, 1)); + acc0 = _mm256_add_epi32(acc0, _mm256_madd_epi16(a_lo, b_lo)); + acc1 = _mm256_add_epi32(acc1, _mm256_madd_epi16(a_hi, b_hi)); + i += 32; + } + + // Horizontal sum of the eight i32 lanes. + let v = _mm256_add_epi32(acc0, acc1); + let s128 = _mm_add_epi32(_mm256_castsi256_si128(v), _mm256_extracti128_si256(v, 1)); + let s64 = _mm_add_epi32(s128, _mm_unpackhi_epi64(s128, s128)); + let s32 = _mm_add_epi32(s64, _mm_shuffle_epi32(s64, 0b01)); + let mut sum = _mm_cvtsi128_si32(s32); + + while i < len { + sum += i32::from(a[i]) * i32::from(b[i]); + i += 1; + } + sum + } +} + /// AVX2 dot product for f32 slices. /// /// # Safety diff --git a/crates/clawhdf5-accel/src/lib.rs b/crates/clawhdf5-accel/src/lib.rs index 95987ef..6fcd9b3 100644 --- a/crates/clawhdf5-accel/src/lib.rs +++ b/crates/clawhdf5-accel/src/lib.rs @@ -122,6 +122,23 @@ pub fn dot_product(a: &[f32], b: &[f32]) -> f32 { } } +/// Dot product of two `i8` slices, widened to `i32`. +/// +/// The kernel behind int8-quantised vector search. Uses the AVX2 path +/// whenever AVX2 is present — including on AVX-512 machines, where it is +/// what the f32 kernels use too on a default build. +pub fn dot_i8(a: &[i8], b: &[i8]) -> i32 { + match detect_backend() { + #[cfg(target_arch = "x86_64")] + // SAFETY: both variants imply AVX2 was detected at runtime (the + // AVX-512 backend is only selected on CPUs that also have AVX2). + Backend::Avx2 | Backend::Avx512 if is_x86_feature_detected!("avx2") => unsafe { + avx2::dot_i8(a, b) + }, + _ => scalar::dot_i8(a, b), + } +} + /// Compute the L2 norm (magnitude) of a vector. pub fn vector_norm(v: &[f32]) -> f32 { dot_product(v, v).sqrt() @@ -713,3 +730,42 @@ mod tests { } } } + +#[cfg(test)] +mod dot_i8_tests { + use super::*; + + fn codes(n: usize, seed: u64) -> Vec { + let mut state = seed; + (0..n) + .map(|_| { + state = state + .wrapping_mul(6_364_136_223_846_793_005) + .wrapping_add(1_442_695_040_888_963_407); + // Full range, including the extremes. + ((state >> 56) as u8) as i8 + }) + .collect() + } + + #[test] + fn dispatched_kernel_matches_scalar_exactly() { + // Integer arithmetic: the SIMD path must agree bit for bit, at every + // length — including ones that are not multiples of the 32-byte block, + // which exercise the tail. + for len in [0, 1, 7, 31, 32, 33, 63, 64, 100, 384, 385, 1536] { + let a = codes(len, 1 + len as u64); + let b = codes(len, 1000 + len as u64); + assert_eq!(dot_i8(&a, &b), scalar::dot_i8(&a, &b), "len {len}"); + } + } + + #[test] + fn extremes_do_not_overflow() { + // -128 * -128 is the largest product; a long run of it must still fit. + let a = vec![-128i8; 4096]; + assert_eq!(dot_i8(&a, &a), 4096 * 128 * 128); + let b = vec![127i8; 4096]; + assert_eq!(dot_i8(&a, &b), -4096 * 128 * 127); + } +} diff --git a/crates/clawhdf5-accel/src/scalar.rs b/crates/clawhdf5-accel/src/scalar.rs index 3c9a507..e8dd5b4 100644 --- a/crates/clawhdf5-accel/src/scalar.rs +++ b/crates/clawhdf5-accel/src/scalar.rs @@ -140,3 +140,33 @@ fn f16_to_f32_soft(h: u16) -> f32 { f32::from_bits(f32_bits) } + +/// Dot product of two `i8` slices, widened to `i32`. +/// +/// `dim` terms of at most `127 * 127` fit an `i32` for any realistic +/// dimension (over 130 000 terms before overflow is possible). +pub fn dot_i8(a: &[i8], b: &[i8]) -> i32 { + assert_eq!(a.len(), b.len()); + // Four independent accumulators over 32-lane blocks: the widening product + // has to sit in a fixed-length chunk for the vectoriser to see it, and the + // separate accumulators keep it off one dependency chain. + const LANE: usize = 8; + let (a_blocks, a_tail) = a.as_chunks::<{ LANE * 4 }>(); + let (b_blocks, b_tail) = b.as_chunks::<{ LANE * 4 }>(); + let mut acc = [0i32; 4]; + for (x, y) in a_blocks.iter().zip(b_blocks) { + for (lane, slot) in acc.iter_mut().enumerate() { + let mut sum = 0i32; + for k in 0..LANE { + sum += i32::from(x[lane * LANE + k]) * i32::from(y[lane * LANE + k]); + } + *slot += sum; + } + } + let tail: i32 = a_tail + .iter() + .zip(b_tail) + .map(|(&x, &y)| i32::from(x) * i32::from(y)) + .sum(); + acc[0] + acc[1] + acc[2] + acc[3] + tail +} diff --git a/crates/clawhdf5-ann/src/hnsw.rs b/crates/clawhdf5-ann/src/hnsw.rs index a4104e4..a526bee 100644 --- a/crates/clawhdf5-ann/src/hnsw.rs +++ b/crates/clawhdf5-ann/src/hnsw.rs @@ -358,31 +358,12 @@ enum Query { Int8(Vec, f32), } -/// Sum of products, widened so it cannot overflow: `dim` terms of at most -/// `127 * 127`, so `i32` suffices for any realistic dimension. +/// Sum of products, widened so it cannot overflow. Runtime-dispatched to the +/// same SIMD backend as the f32 kernels, so the two storages are compared on +/// equal terms. +#[inline] fn dot_i8(a: &[i8], b: &[i8]) -> i32 { - // Four independent accumulators over 32-lane blocks: the widening product - // has to sit in a fixed-length chunk for the vectoriser to see it, and the - // separate accumulators keep it off one dependency chain. - const LANE: usize = 8; - let (a_blocks, a_tail) = a.as_chunks::<{ LANE * 4 }>(); - let (b_blocks, b_tail) = b.as_chunks::<{ LANE * 4 }>(); - let mut acc = [0i32; 4]; - for (x, y) in a_blocks.iter().zip(b_blocks) { - for (lane, slot) in acc.iter_mut().enumerate() { - let mut sum = 0i32; - for k in 0..LANE { - sum += i32::from(x[lane * LANE + k]) * i32::from(y[lane * LANE + k]); - } - *slot += sum; - } - } - let tail: i32 = a_tail - .iter() - .zip(b_tail) - .map(|(&x, &y)| i32::from(x) * i32::from(y)) - .sum(); - acc[0] + acc[1] + acc[2] + acc[3] + tail + clawhdf5_accel::dot_i8(a, b) } /// Magic for [`HnswIndex::graph_to_bytes`].