From 1ac49a57448e4594c0a0878e5a6f6c5ca9f261e0 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Tue, 1 Sep 2026 15:45:43 +0300 Subject: [PATCH 01/10] [MOD-17526] Compute SQ8-FP32 L2 via direct residual accumulation, not cancellation SQ8_FP32_L2Sqr and its SIMD variants computed L2^2 via the algebraic identity ||x||^2 + ||y||^2 - 2*IP(x, y), reading precomputed sums from blob metadata. This catastrophically cancels in FP32 when x and y share a large common offset relative to their spread (e.g. x=[100000,100008], y=[100000,100000]: true L2^2=64, the old kernel could return 0/negative/garbage). Replace it with direct residual accumulation: dequantize each stored byte, subtract the query value, square, and accumulate. The scalar and SSE4/AVX2 paths no longer read y_sum/x_sum_sq/y_sum_sq. The FMA-capable paths (AVX2_FMA, AVX512, NEON, SVE) fuse the dequantize-and-subtract into a single FMA (diff = fma(delta, q, min - y)) rather than computing the subtract separately, since that fusion is the difference between a 10-23% and a 60-80% slowdown on ARM. SQ8-FP16 L2, SQ8-SQ8 L2, and all IP/Cosine kernels are unchanged. --- src/VecSim/spaces/L2/L2.cpp | 32 +++-- src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h | 115 +++++++++++++--- src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h | 113 ++++++++++++--- .../L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h | 110 ++++++++++++--- src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h | 122 +++++++++++++--- src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h | 121 +++++++++++++--- src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h | 130 +++++++++++++++--- 7 files changed, 613 insertions(+), 130 deletions(-) diff --git a/src/VecSim/spaces/L2/L2.cpp b/src/VecSim/spaces/L2/L2.cpp index 253041fb6..971da2ec9 100644 --- a/src/VecSim/spaces/L2/L2.cpp +++ b/src/VecSim/spaces/L2/L2.cpp @@ -19,29 +19,33 @@ using float16 = vecsim_types::float16; using sq8 = vecsim_types::sq8; /* - * Optimized asymmetric SQ8-FP32 L2 squared distance using algebraic identity: - * ||x - y||² = Σx_i² - 2*IP(x, y) + Σy_i² - * = x_sum_squares - 2 * IP(x, y) + y_sum_squares - * where IP(x, y) = min * y_sum + delta * Σ(q_i * y_i) + * Asymmetric SQ8-FP32 L2 squared distance computed via direct residual accumulation: + * ||x - y||² = Σ(dequant(x_i) - y_i)² + * where dequant(x_i) = min_val + delta * q_i + * + * This avoids the algebraic-identity/cancellation approach (||x||² + ||y||² - 2*IP(x, y)), + * which catastrophically cancels in FP32 when x and y share a large common offset relative to + * their spread. * * pVect1 is storage (SQ8): [uint8_t values (dim)] [min_val] [delta] [x_sum] [x_sum_squares] * pVect2 is query (FP32): [float values (dim)] [y_sum] [y_sum_squares] */ float SQ8_FP32_L2Sqr(const void *pVect1v, const void *pVect2v, size_t dimension) { - // Get the raw inner product using the common implementation - const float ip = SQ8_FP32_InnerProduct_Impl(pVect1v, pVect2v, dimension); - // Storage metadata follows a byte payload and is not necessarily float-aligned. const auto *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); - - // Get precomputed sum of squares from query blob (pVect2 is FP32) const auto *pVect2 = static_cast(pVect2v); - const float y_sum_sq = pVect2[dimension + sq8::SUM_SQUARES_QUERY]; - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + const auto *params1 = pVect1 + dimension; + const float min_val = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); + const float delta = load_unaligned(params1 + sq8::DELTA * sizeof(float)); + + float res = 0; + for (size_t i = 0; i < dimension; i++) { + // diff = dequant(x_i) - y_i = min_val + delta * q_i - y_i + float diff = min_val + delta * static_cast(pVect1[i]) - pVect2[i]; + res += diff * diff; + } + return res; } /* diff --git a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h index 0f9d5bde9..4249ea07d 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h @@ -9,38 +9,117 @@ #pragma once #include "VecSim/spaces/space_includes.h" #include "VecSim/spaces/AVX_utils.h" -#include "VecSim/spaces/IP/IP_AVX2_FMA_SQ8_FP32.h" #include "VecSim/types/sq8.h" using sq8 = vecsim_types::sq8; /* - * Optimized asymmetric SQ8 L2 squared distance using algebraic identity: + * Asymmetric SQ8 L2 squared distance computed via direct residual accumulation: * - * ||x - y||² = Σx_i² - 2*IP(x, y) + Σy_i² - * = x_sum_squares - 2 * IP(x, y) + y_sum_squares + * ||x - y||² = Σ(dequant(x_i) - y_i)² + * where dequant(x_i) = min_val + delta * q_i * - * where: - * - IP(x, y) = min * y_sum + delta * Σ(q_i * y_i) (computed via SQ8_FP32_InnerProductImp_FMA) - * - x_sum_squares and y_sum_squares are precomputed + * This avoids the algebraic-identity/cancellation approach, which catastrophically cancels in + * FP32 when x and y share a large common offset relative to their spread. * - * This avoids dequantization in the hot loop. + * This version uses FMA instructions. Critically, the subtract is fused into the FMA + * (diff = fma(delta, q, min - y)) rather than computed separately, which matters for + * performance. */ +// Helper: compute Σ(diff_i²) for 8 elements, where diff_i = dequant(x_i) - y_i. +// pVect1 = SQ8 storage (quantized values), pVect2 = FP32 query. +// min_val_vec/delta_vec are broadcast scalars from the stored vector's metadata. +static inline void L2StepSQ8_FP32_FMA(const uint8_t *&pVect1, const float *&pVect2, __m256 &sum, + __m256 min_val_vec, __m256 delta_vec) { + // Load 8 uint8 elements and convert to float + __m128i v1_128 = _mm_loadl_epi64(reinterpret_cast(pVect1)); + pVect1 += 8; + __m256i v1_256 = _mm256_cvtepu8_epi32(v1_128); + __m256 v1_f = _mm256_cvtepi32_ps(v1_256); + + // Load 8 float elements from query + __m256 v2 = _mm256_loadu_ps(pVect2); + pVect2 += 8; + + // min - y computed once per lane, then fuse the dequantize-and-subtract into a single FMA: + // diff = delta*q + (min - y). + __m256 min_minus_y = _mm256_sub_ps(min_val_vec, v2); + __m256 diff = _mm256_fmadd_ps(delta_vec, v1_f, min_minus_y); + + sum = _mm256_fmadd_ps(diff, diff, sum); +} + // pVect1v = SQ8 storage, pVect2v = FP32 query template // 0..31 float SQ8_FP32_L2SqrSIMD16_AVX2_FMA(const void *pVect1v, const void *pVect2v, size_t dimension) { - // Get the raw inner product using the common SIMD implementation - const float ip = SQ8_FP32_InnerProductImp_FMA(pVect1v, pVect2v, dimension); + const uint8_t *pVect1 = static_cast(pVect1v); // SQ8 storage + const float *pVect2 = static_cast(pVect2v); // FP32 query + const uint8_t *pEnd1 = pVect1 + dimension; + + // Get quantization parameters from stored vector (after quantized data) + const uint8_t *pVect1Base = static_cast(pVect1v); + const auto *params1 = pVect1Base + dimension; + const float min_val_scalar = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); + const float delta_scalar = load_unaligned(params1 + sq8::DELTA * sizeof(float)); + const __m256 min_val_vec = _mm256_set1_ps(min_val_scalar); + const __m256 delta_vec = _mm256_set1_ps(delta_scalar); + + // Initialize sum accumulators. Four accumulators break the FMA dependency chain, letting + // more FMAs be in flight at once. + __m256 sum0 = _mm256_setzero_ps(); + __m256 sum1 = _mm256_setzero_ps(); + __m256 sum2 = _mm256_setzero_ps(); + __m256 sum3 = _mm256_setzero_ps(); + + // Handle residual elements first (0-7 elements). The full-width query load is safe because + // `dim` is at least 8, so the query spans at least 8 floats. + if constexpr (residual % 8) { + __mmask8 constexpr mask = (1 << (residual % 8)) - 1; + + // Load uint8 elements and convert to float + __m128i v1_128 = _mm_loadl_epi64(reinterpret_cast(pVect1)); + pVect1 += residual % 8; + + __m256i v1_256 = _mm256_cvtepu8_epi32(v1_128); + __m256 v1_f = _mm256_cvtepi32_ps(v1_256); + + // Load masked float elements from query + __m256 v2 = my_mm256_maskz_loadu_ps(pVect2); + pVect2 += residual % 8; + + // min - y, then dequantize-and-subtract + __m256 min_minus_y = _mm256_sub_ps(min_val_vec, v2); + __m256 diff = _mm256_fmadd_ps(delta_vec, v1_f, min_minus_y); + + // Masked-out lanes carry garbage (v2 was zeroed, not set to min_val), so blend the + // squared diff with zero for those lanes before accumulating. + __m256 diff_sq = _mm256_mul_ps(diff, diff); + sum0 = _mm256_blend_ps(_mm256_setzero_ps(), diff_sq, mask); + } - // Get precomputed sum of squares from storage blob (pVect1v is SQ8 storage) - const uint8_t *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + // Handle the remaining full 8-element blocks of the residual (compile-time resolved). + if constexpr (residual >= 8) { + L2StepSQ8_FP32_FMA(pVect1, pVect2, sum1, min_val_vec, delta_vec); + } + if constexpr (residual >= 16) { + L2StepSQ8_FP32_FMA(pVect1, pVect2, sum2, min_val_vec, delta_vec); + } + if constexpr (residual >= 24) { + L2StepSQ8_FP32_FMA(pVect1, pVect2, sum3, min_val_vec, delta_vec); + } - // Get precomputed sum of squares from query blob (pVect2v is FP32 query) - const float y_sum_sq = static_cast(pVect2v)[dimension + sq8::SUM_SQUARES_QUERY]; + // We dealt with the residual part. We are left with some multiple of 32 elements. + // In each iteration we calculate 32 elements = 4 chunks of 8. The loop may run zero times + // (dim can be as small as 8). + while (pVect1 < pEnd1) { + L2StepSQ8_FP32_FMA(pVect1, pVect2, sum0, min_val_vec, delta_vec); + L2StepSQ8_FP32_FMA(pVect1, pVect2, sum1, min_val_vec, delta_vec); + L2StepSQ8_FP32_FMA(pVect1, pVect2, sum2, min_val_vec, delta_vec); + L2StepSQ8_FP32_FMA(pVect1, pVect2, sum3, min_val_vec, delta_vec); + } - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + // Reduce to get Σ(diff_i²) + return my_mm256_reduce_add_ps( + _mm256_add_ps(_mm256_add_ps(sum0, sum1), _mm256_add_ps(sum2, sum3))); } diff --git a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h index d08474f71..68bf78bfa 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h @@ -9,38 +9,113 @@ #pragma once #include "VecSim/spaces/space_includes.h" #include "VecSim/spaces/AVX_utils.h" -#include "VecSim/spaces/IP/IP_AVX2_SQ8_FP32.h" #include "VecSim/types/sq8.h" using sq8 = vecsim_types::sq8; /* - * Optimized asymmetric SQ8 L2 squared distance using algebraic identity: + * Asymmetric SQ8 L2 squared distance computed via direct residual accumulation: * - * ||x - y||² = Σx_i² - 2*IP(x, y) + Σy_i² - * = x_sum_squares - 2 * IP(x, y) + y_sum_squares + * ||x - y||² = Σ(dequant(x_i) - y_i)² + * where dequant(x_i) = min_val + delta * q_i * - * where: - * - IP(x, y) = min * y_sum + delta * Σ(q_i * y_i) (computed via SQ8_FP32_InnerProductImp_AVX2) - * - x_sum_squares and y_sum_squares are precomputed - * - * This avoids dequantization in the hot loop. + * This avoids the algebraic-identity/cancellation approach, which catastrophically cancels in + * FP32 when x and y share a large common offset relative to their spread. */ +// Helper: compute Σ(diff_i²) for 8 elements, where diff_i = dequant(x_i) - y_i. +// pVect1 = SQ8 storage (quantized values), pVect2 = FP32 query. +// min_val_vec/delta_vec are broadcast scalars from the stored vector's metadata. +static inline void L2StepSQ8_FP32_AVX2(const uint8_t *&pVect1, const float *&pVect2, __m256 &sum, + __m256 min_val_vec, __m256 delta_vec) { + // Load 8 uint8 elements and convert to float + __m128i v1_128 = _mm_loadl_epi64(reinterpret_cast(pVect1)); + pVect1 += 8; + __m256i v1_256 = _mm256_cvtepu8_epi32(v1_128); + __m256 v1_f = _mm256_cvtepi32_ps(v1_256); + + // Load 8 float elements from query + __m256 v2 = _mm256_loadu_ps(pVect2); + pVect2 += 8; + + // min - y computed once per lane, then dequantize-and-subtract: diff = delta*q + (min - y). + // No FMA in this variant, so mul + add. + __m256 min_minus_y = _mm256_sub_ps(min_val_vec, v2); + __m256 diff = _mm256_add_ps(_mm256_mul_ps(delta_vec, v1_f), min_minus_y); + + sum = _mm256_add_ps(sum, _mm256_mul_ps(diff, diff)); +} + // pVect1v = SQ8 storage, pVect2v = FP32 query template // 0..31 float SQ8_FP32_L2SqrSIMD16_AVX2(const void *pVect1v, const void *pVect2v, size_t dimension) { - // Get the raw inner product using the common SIMD implementation - const float ip = SQ8_FP32_InnerProductImp_AVX2(pVect1v, pVect2v, dimension); + const uint8_t *pVect1 = static_cast(pVect1v); // SQ8 storage + const float *pVect2 = static_cast(pVect2v); // FP32 query + const uint8_t *pEnd1 = pVect1 + dimension; + + // Get quantization parameters from stored vector (after quantized data) + const uint8_t *pVect1Base = static_cast(pVect1v); + const auto *params1 = pVect1Base + dimension; + const float min_val_scalar = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); + const float delta_scalar = load_unaligned(params1 + sq8::DELTA * sizeof(float)); + const __m256 min_val_vec = _mm256_set1_ps(min_val_scalar); + const __m256 delta_vec = _mm256_set1_ps(delta_scalar); + + // Initialize sum accumulators. Four accumulators break the dependency chain, letting more + // ops be in flight at once. + __m256 sum0 = _mm256_setzero_ps(); + __m256 sum1 = _mm256_setzero_ps(); + __m256 sum2 = _mm256_setzero_ps(); + __m256 sum3 = _mm256_setzero_ps(); + + // Handle residual elements first (0-7 elements). The full-width query load is safe because + // `dim` is at least 8, so the query spans at least 8 floats. + if constexpr (residual % 8) { + __mmask8 constexpr mask = (1 << (residual % 8)) - 1; + + // Load uint8 elements and convert to float + __m128i v1_128 = _mm_loadl_epi64(reinterpret_cast(pVect1)); + pVect1 += residual % 8; + + __m256i v1_256 = _mm256_cvtepu8_epi32(v1_128); + __m256 v1_f = _mm256_cvtepi32_ps(v1_256); + + // Load masked float elements from query + __m256 v2 = my_mm256_maskz_loadu_ps(pVect2); + pVect2 += residual % 8; + + // min - y, then dequantize-and-subtract + __m256 min_minus_y = _mm256_sub_ps(min_val_vec, v2); + __m256 diff = _mm256_add_ps(_mm256_mul_ps(delta_vec, v1_f), min_minus_y); + + // Masked-out lanes carry garbage (v2 was zeroed, not set to min_val), so blend the + // squared diff with zero for those lanes before accumulating. + __m256 diff_sq = _mm256_mul_ps(diff, diff); + sum0 = _mm256_blend_ps(_mm256_setzero_ps(), diff_sq, mask); + } - // Get precomputed sum of squares from storage blob (pVect1v is SQ8 storage) - const uint8_t *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + // Handle the remaining full 8-element blocks of the residual (compile-time resolved). + if constexpr (residual >= 8) { + L2StepSQ8_FP32_AVX2(pVect1, pVect2, sum1, min_val_vec, delta_vec); + } + if constexpr (residual >= 16) { + L2StepSQ8_FP32_AVX2(pVect1, pVect2, sum2, min_val_vec, delta_vec); + } + if constexpr (residual >= 24) { + L2StepSQ8_FP32_AVX2(pVect1, pVect2, sum3, min_val_vec, delta_vec); + } - // Get precomputed sum of squares from query blob (pVect2v is FP32 query) - const float y_sum_sq = static_cast(pVect2v)[dimension + sq8::SUM_SQUARES_QUERY]; + // We dealt with the residual part. We are left with some multiple of 32 elements. + // In each iteration we calculate 32 elements = 4 chunks of 8. The loop may run zero times + // (dim can be as small as 8). + while (pVect1 < pEnd1) { + L2StepSQ8_FP32_AVX2(pVect1, pVect2, sum0, min_val_vec, delta_vec); + L2StepSQ8_FP32_AVX2(pVect1, pVect2, sum1, min_val_vec, delta_vec); + L2StepSQ8_FP32_AVX2(pVect1, pVect2, sum2, min_val_vec, delta_vec); + L2StepSQ8_FP32_AVX2(pVect1, pVect2, sum3, min_val_vec, delta_vec); + } - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + // Reduce to get Σ(diff_i²) + return my_mm256_reduce_add_ps( + _mm256_add_ps(_mm256_add_ps(sum0, sum1), _mm256_add_ps(sum2, sum3))); } diff --git a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h index f1c86d689..0c321a228 100644 --- a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h @@ -8,39 +8,109 @@ */ #pragma once #include "VecSim/spaces/space_includes.h" -#include "VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_SQ8_FP32.h" #include "VecSim/types/sq8.h" +#include using sq8 = vecsim_types::sq8; /* - * Optimized asymmetric SQ8 L2 squared distance using algebraic identity: + * Asymmetric SQ8 L2 squared distance computed via direct residual accumulation: * - * ||x - y||² = Σx_i² - 2*IP(x, y) + Σy_i² - * = x_sum_squares - 2 * IP(x, y) + y_sum_squares + * ||x - y||² = Σ(dequant(x_i) - y_i)² + * where dequant(x_i) = min_val + delta * q_i * - * where: - * - IP(x, y) = min * y_sum + delta * Σ(q_i * y_i) (computed via SQ8_FP32_InnerProductImp_AVX512) - * - x_sum_squares and y_sum_squares are precomputed + * This avoids the algebraic-identity/cancellation approach, which catastrophically cancels in + * FP32 when x and y share a large common offset relative to their spread. * - * This avoids dequantization in the hot loop. + * The subtract is fused into the FMA (diff = fma(delta, q, min - y)) rather than computed + * separately, which matters for performance. */ -// pVect1v = SQ8 storage, pVect2v = FP32 query +// Helper: compute Σ(diff_i²) for 16 elements, where diff_i = dequant(x_i) - y_i. +// pVec1 = SQ8 storage (quantized values), pVec2 = FP32 query. +// min_val_vec/delta_vec are broadcast scalars from the stored vector's metadata. +static inline void L2StepSQ8_FP32_AVX512(const uint8_t *&pVec1, const float *&pVec2, __m512 &sum, + __m512 min_val_vec, __m512 delta_vec) { + // Load 16 uint8 elements from quantized vector and convert to float + __m128i v1_128 = _mm_loadu_si128(reinterpret_cast(pVec1)); + __m512i v1_512 = _mm512_cvtepu8_epi32(v1_128); + __m512 v1_f = _mm512_cvtepi32_ps(v1_512); + + // Load 16 float elements from query (pVec2) + __m512 v2 = _mm512_loadu_ps(pVec2); + + // min - y computed once per lane, then fuse the dequantize-and-subtract into a single FMA: + // diff = delta*q + (min - y). + __m512 min_minus_y = _mm512_sub_ps(min_val_vec, v2); + __m512 diff = _mm512_fmadd_ps(delta_vec, v1_f, min_minus_y); + + sum = _mm512_fmadd_ps(diff, diff, sum); + + pVec1 += 16; + pVec2 += 16; +} + +// pVec1v = SQ8 storage, pVec2v = FP32 query template // 0..31 -float SQ8_FP32_L2SqrSIMD16_AVX512F_BW_VL_VNNI(const void *pVect1v, const void *pVect2v, +float SQ8_FP32_L2SqrSIMD16_AVX512F_BW_VL_VNNI(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Get the raw inner product using the common SIMD implementation - const float ip = SQ8_FP32_InnerProductImp_AVX512(pVect1v, pVect2v, dimension); + const uint8_t *pVec1 = static_cast(pVec1v); // SQ8 storage + const float *pVec2 = static_cast(pVec2v); // FP32 query + const uint8_t *pEnd1 = pVec1 + dimension; + + // Get quantization parameters from stored vector (after quantized data) + const uint8_t *pVec1Base = static_cast(pVec1v); + const auto *params1 = pVec1Base + dimension; + const float min_val_scalar = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); + const float delta_scalar = load_unaligned(params1 + sq8::DELTA * sizeof(float)); + const __m512 min_val_vec = _mm512_set1_ps(min_val_scalar); + const __m512 delta_vec = _mm512_set1_ps(delta_scalar); + + // Initialize sum accumulators for Σ(diff_i²). Two accumulators break the FMA dependency + // chain, letting more FMAs be in flight at once. + __m512 sum0 = _mm512_setzero_ps(); + __m512 sum1 = _mm512_setzero_ps(); + + // Handle the sub-16 residual elements first + if constexpr (residual % 16) { + __mmask16 constexpr mask = (1U << (residual % 16)) - 1; + + // Load uint8 elements (safe to load 16 bytes due to the metadata padding after the + // quantized values). The query load is masked, which suppresses faults on masked-out + // lanes, so both loads are safe for any dimension. + __m128i v1_128 = _mm_loadu_si128(reinterpret_cast(pVec1)); + __m512i v1_512 = _mm512_cvtepu8_epi32(v1_128); + __m512 v1_f = _mm512_cvtepi32_ps(v1_512); + + // Load masked float elements from query + __m512 v2 = _mm512_maskz_loadu_ps(mask, pVec2); + + // min - y, then dequantize-and-subtract + __m512 min_minus_y = _mm512_sub_ps(min_val_vec, v2); + __m512 diff = _mm512_fmadd_ps(delta_vec, v1_f, min_minus_y); + + // Masked-out lanes carry garbage (v2 was zeroed, not set to min_val), so mask the + // squared diff to zero for those lanes before accumulating. + sum0 = _mm512_maskz_mul_ps(mask, diff, diff); + + pVec1 += residual % 16; + pVec2 += residual % 16; + } - // Get precomputed sum of squares from storage blob (pVect1v is SQ8 storage) - const uint8_t *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + // Handle the remaining full 16-element block of the residual (compile-time resolved). + if constexpr (residual >= 16) { + L2StepSQ8_FP32_AVX512(pVec1, pVec2, sum1, min_val_vec, delta_vec); + } - // Get precomputed sum of squares from query blob (pVect2v is FP32 query) - const float y_sum_sq = static_cast(pVect2v)[dimension + sq8::SUM_SQUARES_QUERY]; + // We dealt with the residual part. We are left with some multiple of 32 elements. + // In each iteration we calculate 32 elements = 2 chunks of 16. The loop may run zero times + // (dim can be as small as 8). + while (pVec1 < pEnd1) { + L2StepSQ8_FP32_AVX512(pVec1, pVec2, sum0, min_val_vec, delta_vec); + L2StepSQ8_FP32_AVX512(pVec1, pVec2, sum1, min_val_vec, delta_vec); + } - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + // Reduce to get Σ(diff_i²) + __m512 sum = _mm512_add_ps(sum0, sum1); + return _mm512_reduce_add_ps(sum); } diff --git a/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h index f6f7a6bc0..30ece5808 100644 --- a/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h @@ -8,40 +8,124 @@ */ #pragma once #include "VecSim/spaces/space_includes.h" -#include "VecSim/spaces/IP/IP_NEON_SQ8_FP32.h" #include "VecSim/types/sq8.h" #include using sq8 = vecsim_types::sq8; /* - * Optimized asymmetric SQ8 L2 squared distance using algebraic identity: + * Asymmetric SQ8 L2 squared distance computed via direct residual accumulation: * - * ||x - y||² = Σx_i² - 2*IP(x, y) + Σy_i² - * = x_sum_squares - 2 * IP(x, y) + y_sum_squares + * ||x - y||² = Σ(dequant(x_i) - y_i)² + * where dequant(x_i) = min_val + delta * q_i * - * where: - * - IP(x, y) = min * y_sum + delta * Σ(q_i * y_i) (computed via - * SQ8_FP32_InnerProductSIMD16_NEON_IMP) - * - x_sum_squares and y_sum_squares are precomputed + * This avoids the algebraic-identity/cancellation approach, which catastrophically cancels in + * FP32 when x and y share a large common offset relative to their spread. * - * This avoids dequantization in the hot loop. + * The subtract is fused into the FMA (diff = fma(delta, q, min - y)) using the explicit + * vfmaq_f32 intrinsic rather than computed separately (and rather than relying on + * autovectorizing vmlaq_f32-style code), which matters for performance on ARM. */ +// Helper: compute Σ(diff_i²) for 4 elements, where diff_i = dequant(x_i) - y_i. +// pVect1 = SQ8 storage (quantized values), pVect2 = FP32 query. +// min_val_vec/delta_vec are broadcast scalars from the stored vector's metadata. +static inline void L2StepSQ8_FP32_NEON(const uint8_t *&pVect1, const float *&pVect2, + float32x4_t &sum, float32x4_t min_val_vec, + float32x4_t delta_vec) { + // Load 4 uint8 elements and convert to float + uint8x8_t v1_u8 = vld1_u8(pVect1); + pVect1 += 4; + + uint32x4_t v1_u32 = vmovl_u16(vget_low_u16(vmovl_u8(v1_u8))); + float32x4_t v1_f = vcvtq_f32_u32(v1_u32); + + // Load 4 float elements from query + float32x4_t v2 = vld1q_f32(pVect2); + pVect2 += 4; + + // min - y computed once per lane, then fuse the dequantize-and-subtract into a single FMA: + // diff = fma(delta, q, min - y). Uses the explicit vfmaq_f32 intrinsic, not vmlaq_f32. + float32x4_t min_minus_y = vsubq_f32(min_val_vec, v2); + float32x4_t diff = vfmaq_f32(min_minus_y, delta_vec, v1_f); + + sum = vfmaq_f32(sum, diff, diff); +} + // pVect1v = SQ8 storage, pVect2v = FP32 query template // 0..15 float SQ8_FP32_L2SqrSIMD16_NEON(const void *pVect1v, const void *pVect2v, size_t dimension) { - // Get the raw inner product using the common SIMD implementation - const float ip = SQ8_FP32_InnerProductSIMD16_NEON_IMP(pVect1v, pVect2v, dimension); + const uint8_t *pVect1 = static_cast(pVect1v); // SQ8 storage + const float *pVect2 = static_cast(pVect2v); // FP32 query + + // Get quantization parameters from stored vector (after quantized data) + const uint8_t *pVect1Base = static_cast(pVect1v); + const auto *params1 = pVect1Base + dimension; + const float min_val_scalar = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); + const float delta_scalar = load_unaligned(params1 + sq8::DELTA * sizeof(float)); + const float32x4_t min_val_vec = vdupq_n_f32(min_val_scalar); + const float32x4_t delta_vec = vdupq_n_f32(delta_scalar); + + // Multiple accumulators for ILP + float32x4_t sum0 = vdupq_n_f32(0.0f); + float32x4_t sum1 = vdupq_n_f32(0.0f); + float32x4_t sum2 = vdupq_n_f32(0.0f); + float32x4_t sum3 = vdupq_n_f32(0.0f); + + const size_t num_of_chunks = dimension / 16; + + // Process 16 elements at a time in the main loop + for (size_t i = 0; i < num_of_chunks; i++) { + L2StepSQ8_FP32_NEON(pVect1, pVect2, sum0, min_val_vec, delta_vec); + L2StepSQ8_FP32_NEON(pVect1, pVect2, sum1, min_val_vec, delta_vec); + L2StepSQ8_FP32_NEON(pVect1, pVect2, sum2, min_val_vec, delta_vec); + L2StepSQ8_FP32_NEON(pVect1, pVect2, sum3, min_val_vec, delta_vec); + } + + // Handle remaining complete 4-element blocks within residual + if constexpr (residual >= 4) { + L2StepSQ8_FP32_NEON(pVect1, pVect2, sum0, min_val_vec, delta_vec); + } + if constexpr (residual >= 8) { + L2StepSQ8_FP32_NEON(pVect1, pVect2, sum1, min_val_vec, delta_vec); + } + if constexpr (residual >= 12) { + L2StepSQ8_FP32_NEON(pVect1, pVect2, sum2, min_val_vec, delta_vec); + } + + // Handle final residual elements (0-3 elements) + constexpr size_t final_residual = residual % 4; + if constexpr (final_residual > 0) { + // Padding lanes get q=0, y=min_val_scalar, so diff = delta*0 + (min - min) = 0. + float32x4_t v1_f = vdupq_n_f32(0.0f); + float32x4_t v2 = vdupq_n_f32(min_val_scalar); + + if constexpr (final_residual >= 1) { + float q0 = static_cast(pVect1[0]); + v1_f = vld1q_lane_f32(&q0, v1_f, 0); + v2 = vld1q_lane_f32(pVect2, v2, 0); + } + if constexpr (final_residual >= 2) { + float q1 = static_cast(pVect1[1]); + v1_f = vld1q_lane_f32(&q1, v1_f, 1); + v2 = vld1q_lane_f32(pVect2 + 1, v2, 1); + } + if constexpr (final_residual >= 3) { + float q2 = static_cast(pVect1[2]); + v1_f = vld1q_lane_f32(&q2, v1_f, 2); + v2 = vld1q_lane_f32(pVect2 + 2, v2, 2); + } - // Get precomputed sum of squares from storage blob (pVect1v is SQ8 storage) - const uint8_t *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + float32x4_t min_minus_y = vsubq_f32(min_val_vec, v2); + float32x4_t diff = vfmaq_f32(min_minus_y, delta_vec, v1_f); + sum3 = vfmaq_f32(sum3, diff, diff); + } - // Get precomputed sum of squares from query blob (pVect2v is FP32 query) - const float y_sum_sq = static_cast(pVect2v)[dimension + sq8::SUM_SQUARES_QUERY]; + // Combine all four sum accumulators + float32x4_t sum_combined = vaddq_f32(vaddq_f32(sum0, sum1), vaddq_f32(sum2, sum3)); - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + // Horizontal sum to get Σ(diff_i²) + float32x2_t sum_halves = vadd_f32(vget_low_f32(sum_combined), vget_high_f32(sum_combined)); + float32x2_t summed = vpadd_f32(sum_halves, sum_halves); + return vget_lane_f32(summed, 0); } diff --git a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h index 3a3a4d12d..7e56bb141 100644 --- a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h @@ -8,39 +8,120 @@ */ #pragma once #include "VecSim/spaces/space_includes.h" -#include "VecSim/spaces/IP/IP_SSE4_SQ8_FP32.h" #include "VecSim/types/sq8.h" using sq8 = vecsim_types::sq8; /* - * Optimized asymmetric SQ8 L2 squared distance using algebraic identity: + * Asymmetric SQ8 L2 squared distance computed via direct residual accumulation: * - * ||x - y||² = Σx_i² - 2*IP(x, y) + Σy_i² - * = x_sum_squares - 2 * IP(x, y) + y_sum_squares + * ||x - y||² = Σ(dequant(x_i) - y_i)² + * where dequant(x_i) = min_val + delta * q_i * - * where: - * - IP(x, y) = min * y_sum + delta * Σ(q_i * y_i) (computed via - * SQ8_FP32_InnerProductSIMD16_SSE4_IMP) - * - x_sum_squares and y_sum_squares are precomputed - * - * This avoids dequantization in the hot loop. + * This avoids the algebraic-identity/cancellation approach, which catastrophically cancels in + * FP32 when x and y share a large common offset relative to their spread. */ +// Helper: compute Σ(diff_i²) for 4 elements, where diff_i = dequant(x_i) - y_i. +// pVect1 = SQ8 storage (quantized values), pVect2 = FP32 query. +// min_val/delta are broadcast scalars from the stored vector's metadata. +static inline void L2StepSQ8_FP32_SSE4(const uint8_t *&pVect1, const float *&pVect2, + __m128 &sum, __m128 min_val, __m128 delta) { + // Load 4 uint8 elements and convert to float + __m128i v1_i = _mm_cvtepu8_epi32(_mm_cvtsi32_si128(load_unaligned(pVect1))); + pVect1 += 4; + __m128 v1_f = _mm_cvtepi32_ps(v1_i); + + // Load 4 float elements from query + __m128 v2 = _mm_loadu_ps(pVect2); + pVect2 += 4; + + // min - y computed once per lane, then fuse the dequantize-and-subtract: diff = delta*q + + // (min - y). SSE has no FMA, so this is mul + add. + __m128 min_minus_y = _mm_sub_ps(min_val, v2); + __m128 diff = _mm_add_ps(_mm_mul_ps(delta, v1_f), min_minus_y); + + sum = _mm_add_ps(sum, _mm_mul_ps(diff, diff)); +} + // pVect1v = SQ8 storage, pVect2v = FP32 query template // 0..15 float SQ8_FP32_L2SqrSIMD16_SSE4(const void *pVect1v, const void *pVect2v, size_t dimension) { - // Get the raw inner product using the common SIMD implementation - const float ip = SQ8_FP32_InnerProductSIMD16_SSE4_IMP(pVect1v, pVect2v, dimension); + const uint8_t *pVect1 = static_cast(pVect1v); // SQ8 storage + const float *pVect2 = static_cast(pVect2v); // FP32 query + const uint8_t *pEnd1 = pVect1 + dimension; + + // Get quantization parameters from stored vector (after quantized data) + const uint8_t *pVect1Base = static_cast(pVect1v); + const auto *params1 = pVect1Base + dimension; + const float min_val_scalar = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); + const float delta_scalar = load_unaligned(params1 + sq8::DELTA * sizeof(float)); + const __m128 min_val = _mm_set1_ps(min_val_scalar); + const __m128 delta = _mm_set1_ps(delta_scalar); + + // Initialize sum accumulators. Four accumulators break the dependency chain, letting more + // ops be in flight at once. + __m128 sum0 = _mm_setzero_ps(); + __m128 sum1 = _mm_setzero_ps(); + __m128 sum2 = _mm_setzero_ps(); + __m128 sum3 = _mm_setzero_ps(); + + // Process residual elements first (1-3 elements). Loads touch only the residual elements, + // so they are safe for any dimension. + if constexpr (residual % 4) { + float PORTABLE_ALIGN16 q_arr[4] = {0, 0, 0, 0}; + float PORTABLE_ALIGN16 y_arr[4] = {0, 0, 0, 0}; + + if constexpr (residual % 4 >= 1) { + q_arr[0] = static_cast(pVect1[0]); + y_arr[0] = pVect2[0]; + } + if constexpr (residual % 4 >= 2) { + q_arr[1] = static_cast(pVect1[1]); + y_arr[1] = pVect2[1]; + } + if constexpr (residual % 4 >= 3) { + q_arr[2] = static_cast(pVect1[2]); + y_arr[2] = pVect2[2]; + } + // Padding lanes get q=0, y=min_val_scalar, so diff = delta*0 + (min - min) = 0. + for (size_t i = residual % 4; i < 4; i++) { + y_arr[i] = min_val_scalar; + } + + pVect1 += residual % 4; + pVect2 += residual % 4; + + __m128 v1_f = _mm_load_ps(q_arr); + __m128 v2 = _mm_load_ps(y_arr); + __m128 min_minus_y = _mm_sub_ps(min_val, v2); + __m128 diff = _mm_add_ps(_mm_mul_ps(delta, v1_f), min_minus_y); + sum0 = _mm_mul_ps(diff, diff); + } - // Get precomputed sum of squares from storage blob (pVect1v is SQ8 storage) - const uint8_t *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + // Handle remaining residual in chunks of 4 (for residual 4-15) + if constexpr (residual >= 4) { + L2StepSQ8_FP32_SSE4(pVect1, pVect2, sum1, min_val, delta); + } + if constexpr (residual >= 8) { + L2StepSQ8_FP32_SSE4(pVect1, pVect2, sum2, min_val, delta); + } + if constexpr (residual >= 12) { + L2StepSQ8_FP32_SSE4(pVect1, pVect2, sum3, min_val, delta); + } - // Get precomputed sum of squares from query blob (pVect2v is FP32 query) - const float y_sum_sq = static_cast(pVect2v)[dimension + sq8::SUM_SQUARES_QUERY]; + // Process remaining full chunks of 16 elements (4x4). The loop may run zero times + // (dim can be as small as 8). + while (pVect1 < pEnd1) { + L2StepSQ8_FP32_SSE4(pVect1, pVect2, sum0, min_val, delta); + L2StepSQ8_FP32_SSE4(pVect1, pVect2, sum1, min_val, delta); + L2StepSQ8_FP32_SSE4(pVect1, pVect2, sum2, min_val, delta); + L2StepSQ8_FP32_SSE4(pVect1, pVect2, sum3, min_val, delta); + } - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + // Horizontal sum to get Σ(diff_i²) + __m128 sum = _mm_add_ps(_mm_add_ps(sum0, sum1), _mm_add_ps(sum2, sum3)); + float PORTABLE_ALIGN16 TmpRes[4]; + _mm_store_ps(TmpRes, sum); + return TmpRes[0] + TmpRes[1] + TmpRes[2] + TmpRes[3]; } diff --git a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h index a1a800727..da7a4d0d6 100644 --- a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h @@ -8,7 +8,6 @@ */ #pragma once #include "VecSim/spaces/space_includes.h" -#include "VecSim/spaces/IP/IP_SVE_SQ8_FP32.h" #include "VecSim/types/sq8.h" #include @@ -21,35 +20,126 @@ namespace { using sq8 = vecsim_types::sq8; /* - * Optimized asymmetric SQ8-FP32 L2 squared distance using algebraic identity: + * Asymmetric SQ8-FP32 L2 squared distance computed via direct residual accumulation: * - * ||x - y||² = Σx_i² - 2*IP(x, y) + Σy_i² - * = x_sum_squares - 2 * IP(x, y) + y_sum_squares + * ||x - y||² = Σ(dequant(x_i) - y_i)² + * where dequant(x_i) = min_val + delta * q_i * - * where: - * - IP(x, y) = min * y_sum + delta * Σ(q_i * y_i) (computed via - * SQ8_FP32_InnerProductSIMD_SVE_IMP) - * - x_sum_squares and y_sum_squares are precomputed + * This avoids the algebraic-identity/cancellation approach, which catastrophically cancels in + * FP32 when x and y share a large common offset relative to their spread. * - * This avoids dequantization in the hot loop. + * The subtract is fused into the multiply-add (diff = min_minus_y + delta*q) via svmla_f32_x, + * which matters for performance. */ +// Helper: compute Σ(diff_i²) for one SVE vector width, where diff_i = dequant(x_i) - y_i. +// pVect1 = SQ8 storage (quantized values), pVect2 = FP32 query. +// min_val_vec/delta_vec are broadcast scalars from the stored vector's metadata. +static inline void L2StepSQ8_FP32_SVE(const uint8_t *pVect1, const float *pVect2, size_t &offset, + svfloat32_t &sum, const size_t chunk, + svfloat32_t min_val_vec, svfloat32_t delta_vec) { + svbool_t pg = svptrue_b32(); + + // Load uint8 elements and zero-extend to uint32 + svuint32_t v1_u32 = svld1ub_u32(pg, pVect1 + offset); + + // Convert uint32 to float32 + svfloat32_t v1_f = svcvt_f32_u32_x(pg, v1_u32); + + // Load float elements from query + svfloat32_t v2 = svld1_f32(pg, pVect2 + offset); + + // min - y computed once per lane, then fuse the dequantize-and-subtract: + // diff = min_minus_y + delta*q, via a single fused multiply-add. + svfloat32_t min_minus_y = svsub_f32_x(pg, min_val_vec, v2); + svfloat32_t diff = svmla_f32_x(pg, min_minus_y, delta_vec, v1_f); + + sum = svmla_f32_x(pg, sum, diff, diff); + + offset += chunk; +} + // pVect1v = SQ8 storage, pVect2v = FP32 query template float SQ8_FP32_L2SqrSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t dimension) { - // Get the raw inner product using the common SIMD implementation - const float ip = SQ8_FP32_InnerProductSIMD_SVE_IMP( - pVect1v, pVect2v, dimension); + const uint8_t *pVect1 = static_cast(pVect1v); // SQ8 storage + const float *pVect2 = static_cast(pVect2v); // FP32 query + size_t offset = 0; + + svbool_t pg = svptrue_b32(); + + // Get the number of 32-bit elements per vector at runtime + uint64_t chunk = svcntw(); + + // Get quantization parameters from stored vector (after quantized data) + const auto *params1 = pVect1 + dimension; + const float min_val_scalar = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); + const float delta_scalar = load_unaligned(params1 + sq8::DELTA * sizeof(float)); + const svfloat32_t min_val_vec = svdup_f32(min_val_scalar); + const svfloat32_t delta_vec = svdup_f32(delta_scalar); + + // Multiple accumulators for ILP + svfloat32_t sum0 = svdup_f32(0.0f); + svfloat32_t sum1 = svdup_f32(0.0f); + svfloat32_t sum2 = svdup_f32(0.0f); + svfloat32_t sum3 = svdup_f32(0.0f); + + // Handle partial chunk if needed + if constexpr (partial_chunk) { + size_t remaining = dimension % chunk; + if (remaining > 0) { + // Create predicate for the remaining elements + svbool_t pg_partial = + svwhilelt_b32(static_cast(0), static_cast(remaining)); + + // Load uint8 elements and zero-extend to uint32 + svuint32_t v1_u32 = svld1ub_u32(pg_partial, pVect1 + offset); + + // Convert uint32 to float32 + svfloat32_t v1_f = svcvt_f32_u32_z(pg_partial, v1_u32); + + // Load float elements from query with predicate + svfloat32_t v2 = svld1_f32(pg_partial, pVect2); + + // min - y, then dequantize-and-subtract. Inactive lanes of a `_z` (zeroing) + // predicated op become zero, which is exactly what we want when accumulating. + svfloat32_t min_minus_y = svsub_f32_z(pg_partial, min_val_vec, v2); + svfloat32_t diff = svmla_f32_z(pg_partial, min_minus_y, delta_vec, v1_f); + sum0 = svmla_f32_z(pg_partial, sum0, diff, diff); + + offset += remaining; + } + } + + // Process 4 chunks at a time in the main loop + auto chunk_size = 4 * chunk; + const size_t number_of_chunks = + (dimension - (partial_chunk ? dimension % chunk : 0)) / chunk_size; + + for (size_t i = 0; i < number_of_chunks; i++) { + L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum0, chunk, min_val_vec, delta_vec); + L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum1, chunk, min_val_vec, delta_vec); + L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum2, chunk, min_val_vec, delta_vec); + L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum3, chunk, min_val_vec, delta_vec); + } - // Get precomputed sum of squares from storage blob (pVect1v is SQ8 storage) - const uint8_t *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + // Handle remaining steps (0-3) + if constexpr (additional_steps > 0) { + L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum0, chunk, min_val_vec, delta_vec); + } + if constexpr (additional_steps > 1) { + L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum1, chunk, min_val_vec, delta_vec); + } + if constexpr (additional_steps > 2) { + L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum2, chunk, min_val_vec, delta_vec); + } - // Get precomputed sum of squares from query blob (pVect2v is FP32 query) - const float y_sum_sq = static_cast(pVect2v)[dimension + sq8::SUM_SQUARES_QUERY]; + // Combine the accumulators + svfloat32_t sum = svadd_f32_z(pg, sum0, sum1); + sum = svadd_f32_z(pg, sum, sum2); + sum = svadd_f32_z(pg, sum, sum3); - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + // Horizontal sum to get Σ(diff_i²) + return svaddv_f32(pg, sum); } } // namespace From 0e1a4d516b134b78ce038e954bff8b248b0710c9 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Tue, 1 Sep 2026 17:52:21 +0300 Subject: [PATCH 02/10] test(sq8): add MOD-17526 regression tests for SQ8-FP32 L2 Covers the gap that let the original bug ship: translation-invariance (L2(x,y) == L2(x+C,y+C)) at realistic dims (128-1024), the ticket's literal large-offset repro, and a regular-vector sanity check -- all checked against an independent double-precision reference rather than another kernel variant sharing the same formula. --- tests/unit/test_spaces.cpp | 214 +++++++++++++++++++++++++++++++++++++ 1 file changed, 214 insertions(+) diff --git a/tests/unit/test_spaces.cpp b/tests/unit/test_spaces.cpp index fe0138246..9a577f1fd 100644 --- a/tests/unit/test_spaces.cpp +++ b/tests/unit/test_spaces.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include "gtest/gtest.h" #include "VecSim/spaces/space_includes.h" @@ -409,6 +410,219 @@ TEST_F(SpacesTest, SQ8_FP32_odd_dim_unaligned_metadata_test) { } } +/* ==================== MOD-17526 regression tests ==================== + * SQ8_FP32_L2Sqr and its SIMD variants used to compute L2^2 via the algebraic identity + * ||x||^2 + ||y||^2 - 2*IP(x, y), reading precomputed sums from the SQ8 blob metadata. That + * identity catastrophically cancels in FP32 when the two vectors share a large common offset + * relative to their spread (e.g. x=[100000,100008], y=[100000,100000]: true L2^2=64, the old + * kernel could return 0, a negative value, or otherwise unbounded garbage because it computed a + * huge positive value minus a nearly-equal huge positive value in single precision). + * + * The tests below use an independent double-precision reference computed by manually + * dequantizing the SQ8 storage and accumulating squared residuals -- NOT the algebraic identity, + * and NOT test_utils::SQ8_FP32_NotOptimized_L2Sqr or any other kernel variant -- so they would + * have failed against the old kernel and remain meaningful now that both are fixed. They also + * use *relative* tolerances (scaled to the true distance) rather than the ~0.01 absolute + * tolerance used elsewhere in this file, since an absolute tolerance is exactly what let the old + * bug hide: the large-offset case produces errors many orders of magnitude larger than 0.01. + */ + +namespace { + +// Manually dequantizes the SQ8 storage and accumulates the squared residual against the FP32 +// query in double precision. Deliberately independent of both the production kernels and of +// test_utils::SQ8_FP32_NotOptimized_L2Sqr (which accumulates in float) so it serves as ground +// truth rather than another instance of the same computation. +double SQ8_FP32_L2Sqr_DoubleReference(const uint8_t *storage, const float *query, size_t dim) { + const float min_val = load_unaligned(storage + dim + sq8::MIN_VAL * sizeof(float)); + const float delta = load_unaligned(storage + dim + sq8::DELTA * sizeof(float)); + double res = 0.0; + for (size_t i = 0; i < dim; i++) { + const double dequantized = + static_cast(min_val) + static_cast(delta) * storage[i]; + const double diff = static_cast(query[i]) - dequantized; + res += diff * diff; + } + return res; +} + +// Runs the scalar kernel, the runtime-dispatched kernel, and every SIMD variant compiled into +// this binary against `storage`/`query`, asserting each is within `rel_tol` of `expected` +// (relative to max(expected, 1.0), so a tiny true distance isn't swamped by an absolute +// epsilon and a huge one isn't judged by an absolute epsilon that's too tight). +void ExpectSQ8_FP32_L2SqrNear(const uint8_t *storage, const float *query, size_t dim, + double expected, double rel_tol, const std::string &context) { + const double scale = std::max(expected, 1.0); + auto check = [&](double actual, const char *label) { + EXPECT_NEAR(actual, expected, rel_tol * scale) + << context << " [" << label << "]: actual=" << actual << ", expected=" << expected; + }; + + check(SQ8_FP32_L2Sqr(storage, query, dim), "scalar"); + check(L2_SQ8_FP32_GetDistFunc(dim, nullptr)(storage, query, dim), "dispatched"); + + auto optimization = getCpuOptimizationFeatures(); +#ifdef OPT_AVX512_F_BW_VL_VNNI + if (optimization.avx512f && optimization.avx512bw && optimization.avx512vl && + optimization.avx512vnni) { + check(Choose_SQ8_FP32_L2_implementation_AVX512F_BW_VL_VNNI(dim)(storage, query, dim), + "AVX512F_BW_VL_VNNI"); + } +#endif +#ifdef OPT_AVX2_FMA + if (optimization.avx2 && optimization.fma3) { + check(Choose_SQ8_FP32_L2_implementation_AVX2_FMA(dim)(storage, query, dim), "AVX2_FMA"); + } +#endif +#ifdef OPT_AVX2 + if (optimization.avx2) { + check(Choose_SQ8_FP32_L2_implementation_AVX2(dim)(storage, query, dim), "AVX2"); + } +#endif +#ifdef OPT_SSE4 + if (optimization.sse4_1) { + check(Choose_SQ8_FP32_L2_implementation_SSE4(dim)(storage, query, dim), "SSE4"); + } +#endif +#ifdef OPT_SVE2 + if (optimization.sve2) { + check(Choose_SQ8_FP32_L2_implementation_SVE2(dim)(storage, query, dim), "SVE2"); + } +#endif +#ifdef OPT_SVE + if (optimization.sve) { + check(Choose_SQ8_FP32_L2_implementation_SVE(dim)(storage, query, dim), "SVE"); + } +#endif +#ifdef OPT_NEON + if (optimization.asimd) { + check(Choose_SQ8_FP32_L2_implementation_NEON(dim)(storage, query, dim), "NEON"); + } +#endif +} + +// Builds an SQ8 storage blob and an FP32 query for L2 from an explicit pair of float vectors +// (query_values, storage_values). Kept separate from generation so callers can build one base +// pair and derive a shifted pair by literally adding a constant to the same values -- re-sampling +// std::uniform_real_distribution with a shifted [min,max] range and the same seed does NOT +// reproduce "the same vector plus a constant" (it's an independently-shaped draw over the new +// range), which is not what a translation-invariance test needs. +struct SQ8_FP32_L2_TestVectors { + std::vector query; // [float values (dim)] [sum] [sum_squares] + std::vector storage; // [uint8_t values (dim)] [min] [delta] [sum] [sum_squares] +}; + +SQ8_FP32_L2_TestVectors BuildSQ8_FP32_L2_TestVectorsFromValues(const std::vector &query_values, + const std::vector &storage_values) { + const size_t dim = query_values.size(); + SQ8_FP32_L2_TestVectors v; + v.query.resize(dim + sq8::query_metadata_count()); + std::copy(query_values.begin(), query_values.end(), v.query.begin()); + test_utils::preprocess_sq8_fp32_query(v.query.data(), dim); + + v.storage.resize(dim * sizeof(uint8_t) + + sq8::storage_metadata_count() * sizeof(float)); + test_utils::quantize_float_vec_to_sq8_with_metadata(storage_values.data(), dim, v.storage.data()); + return v; +} + +// Draws one base (query, storage) float pair from N(0, spread) and returns it alongside the +// same pair with `offset` added to every element -- a real shift of the identical values, not +// two independently-seeded draws over different ranges. +std::pair +BuildSQ8_FP32_L2_ShiftedPair(size_t dim, float spread, float offset, int seed) { + std::vector query_values(dim), storage_values(dim); + test_utils::populate_float_vec(query_values.data(), dim, seed, -spread, spread); + test_utils::populate_float_vec(storage_values.data(), dim, seed + 1, -spread, spread); + + std::vector shifted_query(dim), shifted_storage(dim); + for (size_t i = 0; i < dim; i++) { + shifted_query[i] = query_values[i] + offset; + shifted_storage[i] = storage_values[i] + offset; + } + + return {BuildSQ8_FP32_L2_TestVectorsFromValues(query_values, storage_values), + BuildSQ8_FP32_L2_TestVectorsFromValues(shifted_query, shifted_storage)}; +} + +} // namespace + +TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_TranslationInvariance) { + // L2 must satisfy L2(x, y) == L2(x + C, y + C) for any constant offset C. The old + // identity-based kernel violated this for large C relative to the vectors' spread; the + // direct-residual kernel must satisfy it (up to FP32 rounding, hence the relative tolerance). + const float spread = 1.0f; + const float large_offset = 100000.0f; // offset/spread == 100000, far past the ~4000 threshold + // where the old identity started to catastrophically + // cancel in FP32. + for (const size_t dim : {128UL, 512UL, 768UL}) { + auto [no_offset, with_offset] = BuildSQ8_FP32_L2_ShiftedPair(dim, spread, large_offset, 1234); + + const double expected_no_offset = SQ8_FP32_L2Sqr_DoubleReference( + no_offset.storage.data(), no_offset.query.data(), dim); + const double expected_with_offset = SQ8_FP32_L2Sqr_DoubleReference( + with_offset.storage.data(), with_offset.query.data(), dim); + + // Sanity check on the reference itself: the two builds only differ by a shared additive + // constant on both vectors, so their true L2^2 must match closely. A small tolerance + // (rather than exact equality) accounts for the shifted float32 inputs themselves losing + // a bit of precision at the larger magnitude (float32 spacing near 1e5 is coarser than + // near 0), which perturbs which byte a value quantizes to right at a code boundary. + ASSERT_NEAR(expected_no_offset, expected_with_offset, + 1e-2 * std::max(expected_no_offset, 1.0)) + << "dim " << dim << ": independent reference is not translation-invariant"; + + ExpectSQ8_FP32_L2SqrNear(no_offset.storage.data(), no_offset.query.data(), dim, + expected_no_offset, 1e-3, "dim " + std::to_string(dim) + ", C=0"); + + // The direct-diff kernel computes `diff = fma(delta, q, min - y)`. `min` and `y` both + // carry the large shared offset, so `min - y` is itself a large-minus-large FP32 + // subtraction (magnitude ~1e5, absolute rounding ~ULP(1e5)/2 ~ 0.008) even though the + // *result* is a small O(1..10) value -- far milder than the old identity's ~1e10-scale + // cancellation, but not zero. This residual is a known, accepted precision limit of the + // direct-diff fix at large offsets (per MOD-17526 benchmarking); removing it entirely + // would require accumulating in double, which was scoped out as a separate product + // decision. Hence a looser tolerance here than the no-offset case above. + ExpectSQ8_FP32_L2SqrNear(with_offset.storage.data(), with_offset.query.data(), dim, + expected_with_offset, 5e-3, + "dim " + std::to_string(dim) + ", C=" + + std::to_string(large_offset)); + } +} + +TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_TicketReproShape) { + // The ticket's literal repro: x=[100000,100008], y=[100000,100000]. True L2^2 = (100000 - + // 100000)^2 + (100008 - 100000)^2 = 64. The old ||x||^2+||y||^2-2*IP identity computed this + // as a huge value minus a nearly-equal huge value in FP32 and could return 0, a negative + // number, or other unbounded garbage instead of 64. + const size_t dim = 2; + const float query_values[dim] = {100000.0f, 100008.0f}; + const float storage_values[dim] = {100000.0f, 100000.0f}; + + std::vector query(dim + sq8::query_metadata_count()); + std::copy(query_values, query_values + dim, query.begin()); + test_utils::preprocess_sq8_fp32_query(query.data(), dim); + + std::vector storage(dim * sizeof(uint8_t) + + sq8::storage_metadata_count() * sizeof(float)); + test_utils::quantize_float_vec_to_sq8_with_metadata(storage_values, dim, storage.data()); + + ExpectSQ8_FP32_L2SqrNear(storage.data(), query.data(), dim, /*expected=*/64.0, 1e-4, + "MOD-17526 ticket repro shape"); +} + +TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_RegularVectorSanity) { + // Sanity check that the fix didn't regress the common case: ordinary small-magnitude + // vectors, no large shared offset, realistic embedding-sized dims. + for (const size_t dim : {128UL, 768UL}) { + auto [v, _unused] = BuildSQ8_FP32_L2_ShiftedPair(dim, /*spread=*/1.0f, /*offset=*/0.0f, 4242); + const double expected = + SQ8_FP32_L2Sqr_DoubleReference(v.storage.data(), v.query.data(), dim); + ExpectSQ8_FP32_L2SqrNear(v.storage.data(), v.query.data(), dim, expected, 1e-3, + "dim " + std::to_string(dim) + ", regular vectors"); + } +} + /* ======================== Tests SQ8-FP16 ========================= */ TEST_F(SpacesTest, SQ8_FP16_ip_no_optimization_norm_func_test) { From a1f7cebbeebf8e1273be9e8074bc75345a9462c9 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Tue, 1 Sep 2026 18:03:18 +0300 Subject: [PATCH 03/10] fix(sq8): reorder scalar SQ8-FP32 L2 to avoid dequant-then-subtract rounding Scalar diff was computed as (min_val + delta*q) - y, left-to-right: the dequantized value gets rounded to FP32 at the large offset's precision before y is subtracted, discarding the residual and reintroducing the cancellation this kernel exists to avoid (only the SIMD kernels used the correct delta*q + (min_val - y) order). Affects x86 dims <8 and any scalar-fallback path. Also fixes the ticket regression test's operands (storage and query were swapped, so the quantized byte was always 0 and delta*q was never exercised), tightens the translation-invariance tolerance now that the dominant error source is gone (min-y is Sterbenz-exact, not a lossy large-offset subtraction as the prior comment claimed), and adds a dedicated regression pinning the exact associativity bug. --- src/VecSim/spaces/L2/L2.cpp | 9 ++++-- tests/unit/test_spaces.cpp | 59 ++++++++++++++++++++++++++++++------- 2 files changed, 55 insertions(+), 13 deletions(-) diff --git a/src/VecSim/spaces/L2/L2.cpp b/src/VecSim/spaces/L2/L2.cpp index 971da2ec9..7ea2ec78a 100644 --- a/src/VecSim/spaces/L2/L2.cpp +++ b/src/VecSim/spaces/L2/L2.cpp @@ -41,8 +41,13 @@ float SQ8_FP32_L2Sqr(const void *pVect1v, const void *pVect2v, size_t dimension) float res = 0; for (size_t i = 0; i < dimension; i++) { - // diff = dequant(x_i) - y_i = min_val + delta * q_i - y_i - float diff = min_val + delta * static_cast(pVect1[i]) - pVect2[i]; + // diff = dequant(x_i) - y_i = delta * q_i + (min_val - y_i). Order matters: min_val and + // y_i are both large and close in magnitude (Sterbenz's lemma makes their subtraction + // exact), so computing that first and adding the small delta*q_i correction preserves the + // residual. Computing (min_val + delta*q_i) - y_i first rounds the dequantized value to + // FP32 at the large offset's precision, discarding the residual before the subtraction + // ever happens -- silently reintroducing the cancellation this kernel exists to avoid. + float diff = delta * static_cast(pVect1[i]) + (min_val - pVect2[i]); res += diff * diff; } return res; diff --git a/tests/unit/test_spaces.cpp b/tests/unit/test_spaces.cpp index 9a577f1fd..e46ca1f0b 100644 --- a/tests/unit/test_spaces.cpp +++ b/tests/unit/test_spaces.cpp @@ -575,16 +575,13 @@ TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_TranslationInvariance) { ExpectSQ8_FP32_L2SqrNear(no_offset.storage.data(), no_offset.query.data(), dim, expected_no_offset, 1e-3, "dim " + std::to_string(dim) + ", C=0"); - // The direct-diff kernel computes `diff = fma(delta, q, min - y)`. `min` and `y` both - // carry the large shared offset, so `min - y` is itself a large-minus-large FP32 - // subtraction (magnitude ~1e5, absolute rounding ~ULP(1e5)/2 ~ 0.008) even though the - // *result* is a small O(1..10) value -- far milder than the old identity's ~1e10-scale - // cancellation, but not zero. This residual is a known, accepted precision limit of the - // direct-diff fix at large offsets (per MOD-17526 benchmarking); removing it entirely - // would require accumulating in double, which was scoped out as a separate product - // decision. Hence a looser tolerance here than the no-offset case above. + // The direct-diff kernel computes `diff = delta*q + (min - y)`. `min` and `y` both carry + // the large shared offset and are within a factor of 2 of each other, so by Sterbenz's + // lemma `min - y` is computed *exactly* in FP32 (no rounding at all) -- it's the small + // `delta*q` correction, not the offset, that's added afterward. So there's no large-scale + // cancellation left to tolerate here; the tolerance stays as tight as the no-offset case. ExpectSQ8_FP32_L2SqrNear(with_offset.storage.data(), with_offset.query.data(), dim, - expected_with_offset, 5e-3, + expected_with_offset, 1e-3, "dim " + std::to_string(dim) + ", C=" + std::to_string(large_offset)); } @@ -596,8 +593,8 @@ TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_TicketReproShape) { // as a huge value minus a nearly-equal huge value in FP32 and could return 0, a negative // number, or other unbounded garbage instead of 64. const size_t dim = 2; - const float query_values[dim] = {100000.0f, 100008.0f}; - const float storage_values[dim] = {100000.0f, 100000.0f}; + const float storage_values[dim] = {100000.0f, 100008.0f}; + const float query_values[dim] = {100000.0f, 100000.0f}; std::vector query(dim + sq8::query_metadata_count()); std::copy(query_values, query_values + dim, query.begin()); @@ -611,6 +608,46 @@ TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_TicketReproShape) { "MOD-17526 ticket repro shape"); } +TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_ScalarAssociativityRegression) { + // Pins a scalar-only regression found in review: `diff = (min_val + delta*q) - y` and + // `diff = delta*q + (min_val - y)` are the same expression algebraically, but not in FP32 -- + // the first rounds the dequantized value to the large offset's precision *before* subtracting + // y, silently discarding the residual the direct-diff kernel exists to preserve; the second + // computes the (Sterbenz-exact) min-y cancellation first and adds the small correction after. + // + // Concrete numbers where this bites: min=100000, delta=8/255, q=1 dequantizes to + // 100000.031372549..., which FP32 rounds to 100000.03125 (float32 ULP near 1e5 is 1/64). + // With y=100000.03125 chosen to land on exactly that rounded value, the buggy order gives + // diff=0 (100000.03125 - 100000.03125), while the correct order gives the true residual + // delta*1 + (100000 - 100000.03125) ~= 0.0001225, matching the double reference. + const size_t dim = 8; + // storage: min=100000 (index 0), max=100008 (index 7) => delta = 8/255, so index 1 + // quantizes to q=1 (dequantizes to 100000 + 8/255 = 100000.031372549...). + const float storage_values[dim] = {100000.0f, 100000.0f + 8.0f / 255.0f, 100000.0f, + 100000.0f, 100000.0f, 100000.0f, + 100000.0f, 100008.0f}; + // query matches storage exactly everywhere except index 1, so every other element + // contributes exactly 0 and the whole result isolates the associativity bug at index 1. + const float query_values[dim] = {100000.0f, 100000.03125f, 100000.0f, 100000.0f, + 100000.0f, 100000.0f, 100000.0f, 100008.0f}; + + std::vector query(dim + sq8::query_metadata_count()); + std::copy(query_values, query_values + dim, query.begin()); + test_utils::preprocess_sq8_fp32_query(query.data(), dim); + + std::vector storage(dim * sizeof(uint8_t) + + sq8::storage_metadata_count() * sizeof(float)); + test_utils::quantize_float_vec_to_sq8_with_metadata(storage_values, dim, storage.data()); + + const double expected = SQ8_FP32_L2Sqr_DoubleReference(storage.data(), query.data(), dim); + ASSERT_GT(expected, 0.0) << "test construction should produce a nonzero true distance"; + // Absolute tolerance, not relative: the true value is tiny (~1.5e-8), and the bug this + // pins is "kernel returns exactly 0 instead of a small nonzero value", not a rounding-scale + // discrepancy -- a relative check against a near-zero expected value isn't meaningful here. + EXPECT_NEAR(SQ8_FP32_L2Sqr(storage.data(), query.data(), dim), expected, 1e-9) + << "scalar kernel: expected=" << expected; +} + TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_RegularVectorSanity) { // Sanity check that the fix didn't regress the common case: ordinary small-magnitude // vectors, no large shared offset, realistic embedding-sized dims. From 07a2d9eb7491c79f452948fb4befbb002d705a77 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 2 Sep 2026 10:10:21 +0300 Subject: [PATCH 04/10] fix(sq8): stop the L2 tests calling x86 kernels below their dim >= 8 contract The test helper calls every ISA kernel directly, bypassing the chooser, and one case runs at dim=2. On x86 the chooser floors SIMD at dim >= 8 ("Optimizations assume at least 8 elements"), and the AVX2 residual path reads a full 32-byte vector through my_mm256_maskz_loadu_ps, which is an unconditional _mm256_loadu_ps -- a 16-byte overread of a dim=2 query blob that ASAN reports. Gate the x86 tiers on the same floor production uses; aarch64 has no such floor (its tail uses per-lane loads) so NEON/SVE stay ungated to keep matching production. Also: spell out the missing `+ offset` on the SVE partial-chunk query load (zero today only because that block runs first), extend the translation-invariance dims to 129/145/513/527 so the large-offset case exercises the residual tails and not just aligned main loops, note that -ffast-math would reassociate the fix away, and apply clang-format. --- src/VecSim/spaces/L2/L2.cpp | 4 ++ src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h | 4 +- src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h | 10 ++-- tests/unit/test_spaces.cpp | 63 +++++++++++++++---------- 4 files changed, 51 insertions(+), 30 deletions(-) diff --git a/src/VecSim/spaces/L2/L2.cpp b/src/VecSim/spaces/L2/L2.cpp index 7ea2ec78a..23744016a 100644 --- a/src/VecSim/spaces/L2/L2.cpp +++ b/src/VecSim/spaces/L2/L2.cpp @@ -27,6 +27,10 @@ using sq8 = vecsim_types::sq8; * which catastrophically cancels in FP32 when x and y share a large common offset relative to * their spread. * + * The operand order in the loop below is load-bearing, not stylistic: it relies on FP addition + * NOT being reassociated. `-ffast-math` / `-Ofast` permit exactly that reassociation and would + * reinstate the bug this kernel fixes. The repo's -O3 builds are safe; keep it that way. + * * pVect1 is storage (SQ8): [uint8_t values (dim)] [min_val] [delta] [x_sum] [x_sum_squares] * pVect2 is query (FP32): [float values (dim)] [y_sum] [y_sum_squares] */ diff --git a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h index 7e56bb141..217ec584e 100644 --- a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h @@ -25,8 +25,8 @@ using sq8 = vecsim_types::sq8; // Helper: compute Σ(diff_i²) for 4 elements, where diff_i = dequant(x_i) - y_i. // pVect1 = SQ8 storage (quantized values), pVect2 = FP32 query. // min_val/delta are broadcast scalars from the stored vector's metadata. -static inline void L2StepSQ8_FP32_SSE4(const uint8_t *&pVect1, const float *&pVect2, - __m128 &sum, __m128 min_val, __m128 delta) { +static inline void L2StepSQ8_FP32_SSE4(const uint8_t *&pVect1, const float *&pVect2, __m128 &sum, + __m128 min_val, __m128 delta) { // Load 4 uint8 elements and convert to float __m128i v1_i = _mm_cvtepu8_epi32(_mm_cvtsi32_si128(load_unaligned(pVect1))); pVect1 += 4; diff --git a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h index da7a4d0d6..a8bcac2e8 100644 --- a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h @@ -36,8 +36,8 @@ using sq8 = vecsim_types::sq8; // pVect1 = SQ8 storage (quantized values), pVect2 = FP32 query. // min_val_vec/delta_vec are broadcast scalars from the stored vector's metadata. static inline void L2StepSQ8_FP32_SVE(const uint8_t *pVect1, const float *pVect2, size_t &offset, - svfloat32_t &sum, const size_t chunk, - svfloat32_t min_val_vec, svfloat32_t delta_vec) { + svfloat32_t &sum, const size_t chunk, svfloat32_t min_val_vec, + svfloat32_t delta_vec) { svbool_t pg = svptrue_b32(); // Load uint8 elements and zero-extend to uint32 @@ -98,8 +98,10 @@ float SQ8_FP32_L2SqrSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t di // Convert uint32 to float32 svfloat32_t v1_f = svcvt_f32_u32_z(pg_partial, v1_u32); - // Load float elements from query with predicate - svfloat32_t v2 = svld1_f32(pg_partial, pVect2); + // Load float elements from query with predicate. `+ offset` is 0 here (this block + // runs before any full-width step), but spell it out to match the load above and + // so the two stay correct if the block order ever changes. + svfloat32_t v2 = svld1_f32(pg_partial, pVect2 + offset); // min - y, then dequantize-and-subtract. Inactive lanes of a `_z` (zeroing) // predicated op become zero, which is exactly what we want when accumulating. diff --git a/tests/unit/test_spaces.cpp b/tests/unit/test_spaces.cpp index e46ca1f0b..3bdfcf155 100644 --- a/tests/unit/test_spaces.cpp +++ b/tests/unit/test_spaces.cpp @@ -462,25 +462,34 @@ void ExpectSQ8_FP32_L2SqrNear(const uint8_t *storage, const float *query, size_t check(L2_SQ8_FP32_GetDistFunc(dim, nullptr)(storage, query, dim), "dispatched"); auto optimization = getCpuOptimizationFeatures(); + + // The x86 kernels are only ever reached through the chooser at dim >= 8 (L2_space.cpp: + // "Optimizations assume at least 8 elements (see the residual handling in the kernels)"). + // Calling them directly below that floor tests a configuration production never produces, + // and the AVX2 residual path reads a full 32-byte vector via my_mm256_maskz_loadu_ps + // (AVX_utils.h) -- an unconditional _mm256_loadu_ps -- which overreads a dim<8 query blob + // and trips ASAN. aarch64 has no such floor (the chooser returns NEON/SVE at any dim, and + // the NEON tail uses per-lane loads), so those tiers stay ungated here to match production. + const bool x86_simd_reachable = dim >= 8; #ifdef OPT_AVX512_F_BW_VL_VNNI - if (optimization.avx512f && optimization.avx512bw && optimization.avx512vl && - optimization.avx512vnni) { + if (x86_simd_reachable && optimization.avx512f && optimization.avx512bw && + optimization.avx512vl && optimization.avx512vnni) { check(Choose_SQ8_FP32_L2_implementation_AVX512F_BW_VL_VNNI(dim)(storage, query, dim), "AVX512F_BW_VL_VNNI"); } #endif #ifdef OPT_AVX2_FMA - if (optimization.avx2 && optimization.fma3) { + if (x86_simd_reachable && optimization.avx2 && optimization.fma3) { check(Choose_SQ8_FP32_L2_implementation_AVX2_FMA(dim)(storage, query, dim), "AVX2_FMA"); } #endif #ifdef OPT_AVX2 - if (optimization.avx2) { + if (x86_simd_reachable && optimization.avx2) { check(Choose_SQ8_FP32_L2_implementation_AVX2(dim)(storage, query, dim), "AVX2"); } #endif #ifdef OPT_SSE4 - if (optimization.sse4_1) { + if (x86_simd_reachable && optimization.sse4_1) { check(Choose_SQ8_FP32_L2_implementation_SSE4(dim)(storage, query, dim), "SSE4"); } #endif @@ -508,12 +517,13 @@ void ExpectSQ8_FP32_L2SqrNear(const uint8_t *storage, const float *query, size_t // reproduce "the same vector plus a constant" (it's an independently-shaped draw over the new // range), which is not what a translation-invariance test needs. struct SQ8_FP32_L2_TestVectors { - std::vector query; // [float values (dim)] [sum] [sum_squares] + std::vector query; // [float values (dim)] [sum] [sum_squares] std::vector storage; // [uint8_t values (dim)] [min] [delta] [sum] [sum_squares] }; -SQ8_FP32_L2_TestVectors BuildSQ8_FP32_L2_TestVectorsFromValues(const std::vector &query_values, - const std::vector &storage_values) { +SQ8_FP32_L2_TestVectors +BuildSQ8_FP32_L2_TestVectorsFromValues(const std::vector &query_values, + const std::vector &storage_values) { const size_t dim = query_values.size(); SQ8_FP32_L2_TestVectors v; v.query.resize(dim + sq8::query_metadata_count()); @@ -522,7 +532,8 @@ SQ8_FP32_L2_TestVectors BuildSQ8_FP32_L2_TestVectorsFromValues(const std::vector v.storage.resize(dim * sizeof(uint8_t) + sq8::storage_metadata_count() * sizeof(float)); - test_utils::quantize_float_vec_to_sq8_with_metadata(storage_values.data(), dim, v.storage.data()); + test_utils::quantize_float_vec_to_sq8_with_metadata(storage_values.data(), dim, + v.storage.data()); return v; } @@ -542,7 +553,7 @@ BuildSQ8_FP32_L2_ShiftedPair(size_t dim, float spread, float offset, int seed) { } return {BuildSQ8_FP32_L2_TestVectorsFromValues(query_values, storage_values), - BuildSQ8_FP32_L2_TestVectorsFromValues(shifted_query, shifted_storage)}; + BuildSQ8_FP32_L2_TestVectorsFromValues(shifted_query, shifted_storage)}; } } // namespace @@ -555,11 +566,15 @@ TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_TranslationInvariance) { const float large_offset = 100000.0f; // offset/spread == 100000, far past the ~4000 threshold // where the old identity started to catastrophically // cancel in FP32. - for (const size_t dim : {128UL, 512UL, 768UL}) { - auto [no_offset, with_offset] = BuildSQ8_FP32_L2_ShiftedPair(dim, spread, large_offset, 1234); - - const double expected_no_offset = SQ8_FP32_L2Sqr_DoubleReference( - no_offset.storage.data(), no_offset.query.data(), dim); + // Dims deliberately mix exact multiples of the SIMD chunk widths with awkward remainders + // (129, 145, 513, 527), so the large-offset cancellation case also runs through the + // residual/masked-lane tails, not just the aligned main loops. + for (const size_t dim : {128UL, 129UL, 145UL, 512UL, 513UL, 527UL, 768UL}) { + auto [no_offset, with_offset] = + BuildSQ8_FP32_L2_ShiftedPair(dim, spread, large_offset, 1234); + + const double expected_no_offset = + SQ8_FP32_L2Sqr_DoubleReference(no_offset.storage.data(), no_offset.query.data(), dim); const double expected_with_offset = SQ8_FP32_L2Sqr_DoubleReference( with_offset.storage.data(), with_offset.query.data(), dim); @@ -580,10 +595,9 @@ TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_TranslationInvariance) { // lemma `min - y` is computed *exactly* in FP32 (no rounding at all) -- it's the small // `delta*q` correction, not the offset, that's added afterward. So there's no large-scale // cancellation left to tolerate here; the tolerance stays as tight as the no-offset case. - ExpectSQ8_FP32_L2SqrNear(with_offset.storage.data(), with_offset.query.data(), dim, - expected_with_offset, 1e-3, - "dim " + std::to_string(dim) + ", C=" + - std::to_string(large_offset)); + ExpectSQ8_FP32_L2SqrNear( + with_offset.storage.data(), with_offset.query.data(), dim, expected_with_offset, 1e-3, + "dim " + std::to_string(dim) + ", C=" + std::to_string(large_offset)); } } @@ -623,13 +637,13 @@ TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_ScalarAssociativityRegression) { const size_t dim = 8; // storage: min=100000 (index 0), max=100008 (index 7) => delta = 8/255, so index 1 // quantizes to q=1 (dequantizes to 100000 + 8/255 = 100000.031372549...). - const float storage_values[dim] = {100000.0f, 100000.0f + 8.0f / 255.0f, 100000.0f, - 100000.0f, 100000.0f, 100000.0f, - 100000.0f, 100008.0f}; + const float storage_values[dim] = { + 100000.0f, 100000.0f + 8.0f / 255.0f, 100000.0f, 100000.0f, 100000.0f, 100000.0f, 100000.0f, + 100008.0f}; // query matches storage exactly everywhere except index 1, so every other element // contributes exactly 0 and the whole result isolates the associativity bug at index 1. const float query_values[dim] = {100000.0f, 100000.03125f, 100000.0f, 100000.0f, - 100000.0f, 100000.0f, 100000.0f, 100008.0f}; + 100000.0f, 100000.0f, 100000.0f, 100008.0f}; std::vector query(dim + sq8::query_metadata_count()); std::copy(query_values, query_values + dim, query.begin()); @@ -652,7 +666,8 @@ TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_RegularVectorSanity) { // Sanity check that the fix didn't regress the common case: ordinary small-magnitude // vectors, no large shared offset, realistic embedding-sized dims. for (const size_t dim : {128UL, 768UL}) { - auto [v, _unused] = BuildSQ8_FP32_L2_ShiftedPair(dim, /*spread=*/1.0f, /*offset=*/0.0f, 4242); + auto [v, _unused] = + BuildSQ8_FP32_L2_ShiftedPair(dim, /*spread=*/1.0f, /*offset=*/0.0f, 4242); const double expected = SQ8_FP32_L2Sqr_DoubleReference(v.storage.data(), v.query.data(), dim); ExpectSQ8_FP32_L2SqrNear(v.storage.data(), v.query.data(), dim, expected, 1e-3, From 47c4f9d3de1040ad1f1703273f94b464989fba4c Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 2 Sep 2026 10:30:54 +0300 Subject: [PATCH 05/10] perf(sq8): stop the SQ8-FP32 L2 kernels paying for loads and folds they don't need Four independent costs, each measured with the repo's own bm_spaces_sq8_fp32 (median of 5 reps, cv <= 1%). NEON: the 4-element step read 8 bytes (vld1_u8) but consumed only the low 4 and advanced by 4, so consecutive steps re-read half of every load. Widening the main loop to one 16-byte load feeding all four accumulators, with the per-lane arithmetic and accumulator mapping unchanged. Neoverse-N1: dim16 10.9->8.59ns, dim64 34.5->25.8ns, dim256 128->95.2ns, dim1024 502->370ns (-21% to -26%). The 4-element helper stays for the tail. SVE: `chunk` is svcntw(), a runtime value, so `dimension % chunk` and `/ chunk_size` emitted real udiv (~12-20 cycles, unpipelined) on every call. They were also redundant -- CHOOSE_SVE_IMPLEMENTATION already divides once at chooser time and passes the answers down as template params. Restructured to full vectors first with a compared bound and an svwhilelt-predicated tail, which also keeps every unpredicated load at a multiple of the vector length instead of pushing them off a leading partial chunk. Isolated like-for-like compile of two instantiations: 3 udiv before, 0 after; 0 udiv across all 32 instantiations in both SVE.cpp.o and SVE2.cpp.o. Loop algebra verified exhaustively off-hardware for chunk in {4,8,16,32,64} x dim in [1,4096]: every element covered exactly once, no overread. SSE4: the residual path stored scalars into two 16-byte stack arrays and immediately reloaded them with _mm_load_ps, which cannot be store-to-load forwarded from narrower stores. It showed in the residual sweep as a fixed penalty on every dim where residual % 4 != 0 (125-127ns) versus the dims that skip the path entirely (116-117ns). Now both operands load at full width -- in bounds because the kernel is unreachable below dim 8 -- with a compile-time lane mask zeroing the lanes the main loop will handle. The sweep is flat at 116-119ns. AVX2 / AVX2+FMA: my_mm256_reduce_add_ps spills 8 floats and sums them with 7 dependent scalar adds, a fixed cost that dominates at small dims. Added my_mm256_reduce_add_ps_tree (three in-register steps) alongside it rather than changing it, since the original has callers across the repo that would each need their own benchmarking. AVX2+FMA dim16 4.42->3.10ns (-30%), dim64 9.05->7.65ns, dim1024 85.4->81.9ns (-4%); AVX2 dim16 5.25->3.91ns (-26%). AVX512 keeps _mm512_reduce_add_ps, which is already a tree. Verified: 595/595 SQ8 tests and 595/595 under ASAN on x86 (AVX512/AVX2+FMA/ AVX2/SSE4), 593/593 on aarch64 (NEON). Left alone deliberately, per measurements: AVX512 stays at 2 accumulators, the scalar path stays at 1 (it only ever sees dims 1-7), and the residual arithmetic itself is unchanged. --- src/VecSim/spaces/AVX_utils.h | 16 +++++ src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h | 4 +- src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h | 4 +- src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h | 46 ++++++++++++-- src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h | 54 ++++++++-------- src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h | 70 ++++++++++----------- 6 files changed, 124 insertions(+), 70 deletions(-) diff --git a/src/VecSim/spaces/AVX_utils.h b/src/VecSim/spaces/AVX_utils.h index 2fe0b904e..a87081487 100644 --- a/src/VecSim/spaces/AVX_utils.h +++ b/src/VecSim/spaces/AVX_utils.h @@ -35,3 +35,19 @@ static inline float my_mm256_reduce_add_ps(__m256 x) { return TmpRes[0] + TmpRes[1] + TmpRes[2] + TmpRes[3] + TmpRes[4] + TmpRes[5] + TmpRes[6] + TmpRes[7]; } + +// Same result as my_mm256_reduce_add_ps, folded in-register instead of through the stack. +// +// The version above spills 8 floats and sums them with 7 *dependent* scalar adds, so it pays a +// store-to-load stall plus a serial ~7-add latency chain every call. That is a fixed cost, so +// it dominates at small dimensions. This does it in three in-register steps. +// +// Reassociating changes which partial sums are formed, so the low bits of the result can +// differ. Added alongside the original rather than replacing it: the original has many callers +// across the repo that would each need their own benchmarking and tolerance review. +static inline float my_mm256_reduce_add_ps_tree(__m256 x) { + __m128 sum128 = _mm_add_ps(_mm256_castps256_ps128(x), _mm256_extractf128_ps(x, 1)); + sum128 = _mm_add_ps(sum128, _mm_movehl_ps(sum128, sum128)); + sum128 = _mm_add_ss(sum128, _mm_shuffle_ps(sum128, sum128, 0x1)); + return _mm_cvtss_f32(sum128); +} diff --git a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h index 4249ea07d..285bb2ab7 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h @@ -120,6 +120,8 @@ float SQ8_FP32_L2SqrSIMD16_AVX2_FMA(const void *pVect1v, const void *pVect2v, si } // Reduce to get Σ(diff_i²) - return my_mm256_reduce_add_ps( + // Tree reduction: the fixed cost of the stack-based fold is a large share of + // total time at small dimensions. See my_mm256_reduce_add_ps_tree in AVX_utils.h. + return my_mm256_reduce_add_ps_tree( _mm256_add_ps(_mm256_add_ps(sum0, sum1), _mm256_add_ps(sum2, sum3))); } diff --git a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h index 68bf78bfa..f3f1a31f8 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h @@ -116,6 +116,8 @@ float SQ8_FP32_L2SqrSIMD16_AVX2(const void *pVect1v, const void *pVect2v, size_t } // Reduce to get Σ(diff_i²) - return my_mm256_reduce_add_ps( + // Tree reduction: the fixed cost of the stack-based fold is a large share of + // total time at small dimensions. See my_mm256_reduce_add_ps_tree in AVX_utils.h. + return my_mm256_reduce_add_ps_tree( _mm256_add_ps(_mm256_add_ps(sum0, sum1), _mm256_add_ps(sum2, sum3))); } diff --git a/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h index 30ece5808..8c2287476 100644 --- a/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h @@ -52,6 +52,45 @@ static inline void L2StepSQ8_FP32_NEON(const uint8_t *&pVect1, const float *&pVe sum = vfmaq_f32(sum, diff, diff); } +// Helper: same arithmetic as above, but for 16 elements off a single 16-byte load. +// +// The 4-element helper reads 8 bytes (vld1_u8) and consumes only the low 4 (vget_low_u16), +// then advances the pointer by 4, so back-to-back calls re-read half of what they just loaded +// and defeat load pairing. Widening one load into all four float32x4_t groups removes that. +// Per-lane arithmetic and accumulator assignment are unchanged (group g still lands in sum, +// still fma(delta, q, min - y)), so results stay bit-identical to the 4-element path. +static inline void L2Step16SQ8_FP32_NEON(const uint8_t *&pVect1, const float *&pVect2, + float32x4_t &sum0, float32x4_t &sum1, float32x4_t &sum2, + float32x4_t &sum3, float32x4_t min_val_vec, + float32x4_t delta_vec) { + uint8x16_t v1_u8 = vld1q_u8(pVect1); + pVect1 += 16; + + const uint16x8_t wide_lo = vmovl_u8(vget_low_u8(v1_u8)); + const uint16x8_t wide_hi = vmovl_u8(vget_high_u8(v1_u8)); + + const float32x4_t q0 = vcvtq_f32_u32(vmovl_u16(vget_low_u16(wide_lo))); + const float32x4_t q1 = vcvtq_f32_u32(vmovl_u16(vget_high_u16(wide_lo))); + const float32x4_t q2 = vcvtq_f32_u32(vmovl_u16(vget_low_u16(wide_hi))); + const float32x4_t q3 = vcvtq_f32_u32(vmovl_u16(vget_high_u16(wide_hi))); + + const float32x4_t y0 = vld1q_f32(pVect2); + const float32x4_t y1 = vld1q_f32(pVect2 + 4); + const float32x4_t y2 = vld1q_f32(pVect2 + 8); + const float32x4_t y3 = vld1q_f32(pVect2 + 12); + pVect2 += 16; + + const float32x4_t d0 = vfmaq_f32(vsubq_f32(min_val_vec, y0), delta_vec, q0); + const float32x4_t d1 = vfmaq_f32(vsubq_f32(min_val_vec, y1), delta_vec, q1); + const float32x4_t d2 = vfmaq_f32(vsubq_f32(min_val_vec, y2), delta_vec, q2); + const float32x4_t d3 = vfmaq_f32(vsubq_f32(min_val_vec, y3), delta_vec, q3); + + sum0 = vfmaq_f32(sum0, d0, d0); + sum1 = vfmaq_f32(sum1, d1, d1); + sum2 = vfmaq_f32(sum2, d2, d2); + sum3 = vfmaq_f32(sum3, d3, d3); +} + // pVect1v = SQ8 storage, pVect2v = FP32 query template // 0..15 float SQ8_FP32_L2SqrSIMD16_NEON(const void *pVect1v, const void *pVect2v, size_t dimension) { @@ -74,12 +113,9 @@ float SQ8_FP32_L2SqrSIMD16_NEON(const void *pVect1v, const void *pVect2v, size_t const size_t num_of_chunks = dimension / 16; - // Process 16 elements at a time in the main loop + // Process 16 elements at a time in the main loop, one 16-byte load per iteration. for (size_t i = 0; i < num_of_chunks; i++) { - L2StepSQ8_FP32_NEON(pVect1, pVect2, sum0, min_val_vec, delta_vec); - L2StepSQ8_FP32_NEON(pVect1, pVect2, sum1, min_val_vec, delta_vec); - L2StepSQ8_FP32_NEON(pVect1, pVect2, sum2, min_val_vec, delta_vec); - L2StepSQ8_FP32_NEON(pVect1, pVect2, sum3, min_val_vec, delta_vec); + L2Step16SQ8_FP32_NEON(pVect1, pVect2, sum0, sum1, sum2, sum3, min_val_vec, delta_vec); } // Handle remaining complete 4-element blocks within residual diff --git a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h index 217ec584e..a461c79a6 100644 --- a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h @@ -66,36 +66,36 @@ float SQ8_FP32_L2SqrSIMD16_SSE4(const void *pVect1v, const void *pVect2v, size_t __m128 sum2 = _mm_setzero_ps(); __m128 sum3 = _mm_setzero_ps(); - // Process residual elements first (1-3 elements). Loads touch only the residual elements, - // so they are safe for any dimension. + // Process residual elements first (1-3 elements). + // + // Both operands are loaded at full width and the lanes past the residual are masked off, + // rather than staged through the stack. The previous form stored scalars into two 16-byte + // stack arrays and immediately reloaded them with _mm_load_ps; a 16-byte load cannot be + // store-to-load forwarded from four narrower stores, so it stalls. That showed up in the + // residual benchmark sweep as a fixed ~8-9 ns penalty on every dim where residual % 4 != 0. + // + // The wide loads are in bounds because this kernel is only reachable at dim >= 8 (the x86 + // chooser floors SIMD there), so both operands have at least 8 elements ahead of offset 0, + // and the metadata trailing each blob keeps even a 16-byte read inside the allocation. if constexpr (residual % 4) { - float PORTABLE_ALIGN16 q_arr[4] = {0, 0, 0, 0}; - float PORTABLE_ALIGN16 y_arr[4] = {0, 0, 0, 0}; - - if constexpr (residual % 4 >= 1) { - q_arr[0] = static_cast(pVect1[0]); - y_arr[0] = pVect2[0]; - } - if constexpr (residual % 4 >= 2) { - q_arr[1] = static_cast(pVect1[1]); - y_arr[1] = pVect2[1]; - } - if constexpr (residual % 4 >= 3) { - q_arr[2] = static_cast(pVect1[2]); - y_arr[2] = pVect2[2]; - } - // Padding lanes get q=0, y=min_val_scalar, so diff = delta*0 + (min - min) = 0. - for (size_t i = residual % 4; i < 4; i++) { - y_arr[i] = min_val_scalar; - } - - pVect1 += residual % 4; - pVect2 += residual % 4; - - __m128 v1_f = _mm_load_ps(q_arr); - __m128 v2 = _mm_load_ps(y_arr); + constexpr unsigned char r = residual % 4; + + __m128i v1_i = _mm_cvtepu8_epi32(_mm_cvtsi32_si128(load_unaligned(pVect1))); + __m128 v1_f = _mm_cvtepi32_ps(v1_i); + __m128 v2 = _mm_loadu_ps(pVect2); + __m128 min_minus_y = _mm_sub_ps(min_val, v2); __m128 diff = _mm_add_ps(_mm_mul_ps(delta, v1_f), min_minus_y); + + // Lanes >= r hold elements the main loop will process; zero them so this step adds + // nothing for them. r is a compile-time value, so the mask is a constant. + const __m128 lane_mask = _mm_castsi128_ps( + _mm_set_epi32(r > 3 ? -1 : 0, r > 2 ? -1 : 0, r > 1 ? -1 : 0, r > 0 ? -1 : 0)); + diff = _mm_and_ps(diff, lane_mask); + + pVect1 += r; + pVect2 += r; + sum0 = _mm_mul_ps(diff, diff); } diff --git a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h index a8bcac2e8..7c3a41ca9 100644 --- a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h @@ -84,48 +84,30 @@ float SQ8_FP32_L2SqrSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t di svfloat32_t sum2 = svdup_f32(0.0f); svfloat32_t sum3 = svdup_f32(0.0f); - // Handle partial chunk if needed - if constexpr (partial_chunk) { - size_t remaining = dimension % chunk; - if (remaining > 0) { - // Create predicate for the remaining elements - svbool_t pg_partial = - svwhilelt_b32(static_cast(0), static_cast(remaining)); - - // Load uint8 elements and zero-extend to uint32 - svuint32_t v1_u32 = svld1ub_u32(pg_partial, pVect1 + offset); - - // Convert uint32 to float32 - svfloat32_t v1_f = svcvt_f32_u32_z(pg_partial, v1_u32); - - // Load float elements from query with predicate. `+ offset` is 0 here (this block - // runs before any full-width step), but spell it out to match the load above and - // so the two stay correct if the block order ever changes. - svfloat32_t v2 = svld1_f32(pg_partial, pVect2 + offset); - - // min - y, then dequantize-and-subtract. Inactive lanes of a `_z` (zeroing) - // predicated op become zero, which is exactly what we want when accumulating. - svfloat32_t min_minus_y = svsub_f32_z(pg_partial, min_val_vec, v2); - svfloat32_t diff = svmla_f32_z(pg_partial, min_minus_y, delta_vec, v1_f); - sum0 = svmla_f32_z(pg_partial, sum0, diff, diff); - - offset += remaining; - } - } - - // Process 4 chunks at a time in the main loop - auto chunk_size = 4 * chunk; - const size_t number_of_chunks = - (dimension - (partial_chunk ? dimension % chunk : 0)) / chunk_size; - - for (size_t i = 0; i < number_of_chunks; i++) { + // Full-width groups first, predicated tail last, and every bound compared rather than + // divided. + // + // `chunk` is a runtime value (svcntw), so `dimension % chunk` and `/ chunk_size` compile to + // real `udiv` instructions -- ~12-20 cycles each and not pipelined, on a function called + // thousands of times per query. They are also redundant: CHOOSE_SVE_IMPLEMENTATION already + // divides once, at chooser time, and hands the results down as `partial_chunk` and + // `additional_steps`. Doing the full vectors first additionally keeps every unpredicated + // load at a multiple of the vector length; a leading partial chunk would push all of them + // to a non-VL-multiple offset. + // + // This is why the shape here differs from the prefix-first sibling SQ8 SVE kernels. Given + // dimension = k*chunk + r (r = dimension % chunk, so r > 0 exactly when partial_chunk), + // the loop below runs floor(k/4) times, leaving (k % 4) == additional_steps full vectors + // plus r tail elements. + const size_t chunk_size = 4 * chunk; + while (offset + chunk_size <= dimension) { L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum0, chunk, min_val_vec, delta_vec); L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum1, chunk, min_val_vec, delta_vec); L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum2, chunk, min_val_vec, delta_vec); L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum3, chunk, min_val_vec, delta_vec); } - // Handle remaining steps (0-3) + // Handle remaining full-width steps (0-3), resolved at compile time. if constexpr (additional_steps > 0) { L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum0, chunk, min_val_vec, delta_vec); } @@ -136,6 +118,22 @@ float SQ8_FP32_L2SqrSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t di L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum2, chunk, min_val_vec, delta_vec); } + // Predicated tail for the final partial vector. svwhilelt derives the predicate from the + // current offset against `dimension`, so no residual arithmetic is needed, and inactive + // lanes of the `_z` (zeroing) forms contribute 0 to the accumulator. + if constexpr (partial_chunk) { + svbool_t pg_tail = + svwhilelt_b32(static_cast(offset), static_cast(dimension)); + + svuint32_t v1_u32 = svld1ub_u32(pg_tail, pVect1 + offset); + svfloat32_t v1_f = svcvt_f32_u32_z(pg_tail, v1_u32); + svfloat32_t v2 = svld1_f32(pg_tail, pVect2 + offset); + + svfloat32_t min_minus_y = svsub_f32_z(pg_tail, min_val_vec, v2); + svfloat32_t diff = svmla_f32_z(pg_tail, min_minus_y, delta_vec, v1_f); + sum3 = svmla_f32_z(pg_tail, sum3, diff, diff); + } + // Combine the accumulators svfloat32_t sum = svadd_f32_z(pg, sum0, sum1); sum = svadd_f32_z(pg, sum, sum2); From fc7995d4f0ae02a75cb24a5e7e9c25ca2b2c77ed Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 2 Sep 2026 12:27:31 +0300 Subject: [PATCH 06/10] fix(sq8): keep the SVE L2 tail from zeroing live accumulator lanes The restructured tail accumulated with svmla_f32_z, whose zeroing form clears every lane outside the predicate. sum3 already holds full-vector partial sums by then -- the main loop's fourth step writes it -- so any dimension with a partial chunk and at least one main-loop iteration (dim >= 4*chunk && dim % chunk != 0) silently dropped part of the distance. The pre-restructure kernel used _z safely here only because the partial block ran first, against an all-zero sum0; moving it after the main loop removed that precondition. Use the merging form instead. The temporaries above stay _z since they are fresh values rather than accumulators. Caught by Cursor Bugbot and then reproduced on Graviton4 (Neoverse-V2, 128-bit VL so chunk=4): the parameterized SQ8_FP32_L2Sqr sweep failed on exactly the dims where dim % 4 != 0 and passed on the rest. With the fix, 593/593 SQ8 tests pass on that host -- the first run of these kernels on hardware that actually has SVE. My off-hardware coverage check did not catch this: it verified which elements each step touches, not what the predicate does to inactive lanes. --- src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h index 7c3a41ca9..908b23130 100644 --- a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h @@ -131,7 +131,13 @@ float SQ8_FP32_L2SqrSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t di svfloat32_t min_minus_y = svsub_f32_z(pg_tail, min_val_vec, v2); svfloat32_t diff = svmla_f32_z(pg_tail, min_minus_y, delta_vec, v1_f); - sum3 = svmla_f32_z(pg_tail, sum3, diff, diff); + + // Merging (_m), not zeroing (_z): sum3 already holds full-vector partial sums from the + // main loop, and _z would zero every lane outside pg_tail, silently dropping them. _z is + // only safe on a freshly zeroed accumulator, which is why the pre-restructure kernel + // could use it here -- back then this block ran first, against an all-zero sum0. + // The temporaries above stay _z because they are fresh values, not accumulators. + sum3 = svmla_f32_m(pg_tail, sum3, diff, diff); } // Combine the accumulators From 891c8762d10564d9694d74e8c55220c9e22a128c Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 2 Sep 2026 12:39:34 +0300 Subject: [PATCH 07/10] docs(sq8): record the measured SVE loop-shape trade-off The restructure is not a uniform win and the numbers should be in the file rather than only in the PR: on Graviton4 (128-bit VL, median of 9, cv <= 0.33%) SVE2 improves everywhere while plain SVE improves only where the dimension divides evenly by the vector length and costs ~8% where it does not. Kept because SVE2 is preferred when both tiers are present and real embedding dims land on the faster path. Also records the tail variant that measured worse, and why a bitmask cannot replace the modulo (SVE vector length is a multiple of 128 bits, not necessarily a power of two). --- src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h index 908b23130..2e6950363 100644 --- a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h @@ -99,6 +99,18 @@ float SQ8_FP32_L2SqrSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t di // dimension = k*chunk + r (r = dimension % chunk, so r > 0 exactly when partial_chunk), // the loop below runs floor(k/4) times, leaving (k % 4) == additional_steps full vectors // plus r tail elements. + // + // Measured on Graviton4 (Neoverse-V2, 128-bit VL), median of 9, cv <= 0.33%, against the + // prefix-first shape: SVE2 is faster everywhere (dim 1024 160 -> 141ns, dim 513 81.9 -> + // 79.4ns), and SVE is faster where dimension is a multiple of the vector length (dim 1024 + // 164 -> 153ns) but ~8% slower at dims that leave a tail (dim 513 78.3 -> 84.6ns). Kept + // because SVE2 is the tier the chooser prefers when both are present, and because real + // embedding dims (128/256/384/512/768/1024/1536) all divide evenly at 128- and 256-bit VL, + // landing on the faster path. Feeding the tail through its own accumulator instead of + // svmla_f32_m was also tried and measured slightly worse (dim 513 85.7ns). + // + // Note `dimension & (chunk - 1)` is NOT a valid substitute for the modulo: SVE permits any + // vector length that is a multiple of 128 bits, not only powers of two. const size_t chunk_size = 4 * chunk; while (offset + chunk_size <= dimension) { L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum0, chunk, min_val_vec, delta_vec); From 766aabcd0dbbcfb062eaf7831a33c62bf9b1d406 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 2 Sep 2026 12:50:22 +0300 Subject: [PATCH 08/10] style(sq8): cut the comment volume down and match the for-loop convention Comments on the changed lines drop from 226 to 128. Removed the narration that restates the intrinsic on the next line, and moved the benchmark tables and the rejected-alternative notes out of the sources -- they are already in the commit messages and the PR, which is where they belong. What stays is the reasoning a reader cannot recover from the code: why the operand order is load-bearing, why the SVE tail merges instead of zeroing, why the wide residual load is in bounds, and why the SVE loop shape differs from its siblings. Also switch the SVE main loop from while to for, matching NEON and the sibling SVE kernels. The step helper no longer mutates the caller's offset: it takes one by value, so the cursor advances in the loop header and each unrolled step names the sub-offset it reads. Same instruction sequence -- 0 udiv in both SVE.cpp.o and SVE2.cpp.o, and Graviton4 medians are unchanged (SVE 1024 153ns / 513 84.6ns, SVE2 1024 141ns / 513 79.4ns). 595/595 SQ8 on x86, 593/593 on Graviton4. --- src/VecSim/spaces/AVX_utils.h | 13 ++-- src/VecSim/spaces/L2/L2.cpp | 14 ++-- src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h | 32 +++------ src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h | 17 ++--- src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h | 94 ++++++++---------------- tests/unit/test_spaces.cpp | 95 ++++++++----------------- 6 files changed, 82 insertions(+), 183 deletions(-) diff --git a/src/VecSim/spaces/AVX_utils.h b/src/VecSim/spaces/AVX_utils.h index a87081487..b631ee5d6 100644 --- a/src/VecSim/spaces/AVX_utils.h +++ b/src/VecSim/spaces/AVX_utils.h @@ -36,15 +36,10 @@ static inline float my_mm256_reduce_add_ps(__m256 x) { TmpRes[7]; } -// Same result as my_mm256_reduce_add_ps, folded in-register instead of through the stack. -// -// The version above spills 8 floats and sums them with 7 *dependent* scalar adds, so it pays a -// store-to-load stall plus a serial ~7-add latency chain every call. That is a fixed cost, so -// it dominates at small dimensions. This does it in three in-register steps. -// -// Reassociating changes which partial sums are formed, so the low bits of the result can -// differ. Added alongside the original rather than replacing it: the original has many callers -// across the repo that would each need their own benchmarking and tolerance review. +// As my_mm256_reduce_add_ps, but folded in-register in three steps instead of spilling 8 floats +// and summing them with 7 dependent scalar adds. That fixed cost dominates at small dimensions. +// Reassociating changes the low bits, so this is added alongside the original rather than +// replacing it -- the original's other callers would each need their own benchmarking. static inline float my_mm256_reduce_add_ps_tree(__m256 x) { __m128 sum128 = _mm_add_ps(_mm256_castps256_ps128(x), _mm256_extractf128_ps(x, 1)); sum128 = _mm_add_ps(sum128, _mm_movehl_ps(sum128, sum128)); diff --git a/src/VecSim/spaces/L2/L2.cpp b/src/VecSim/spaces/L2/L2.cpp index 23744016a..3ba75fb4c 100644 --- a/src/VecSim/spaces/L2/L2.cpp +++ b/src/VecSim/spaces/L2/L2.cpp @@ -27,9 +27,8 @@ using sq8 = vecsim_types::sq8; * which catastrophically cancels in FP32 when x and y share a large common offset relative to * their spread. * - * The operand order in the loop below is load-bearing, not stylistic: it relies on FP addition - * NOT being reassociated. `-ffast-math` / `-Ofast` permit exactly that reassociation and would - * reinstate the bug this kernel fixes. The repo's -O3 builds are safe; keep it that way. + * The operand order below relies on FP addition NOT being reassociated, so `-ffast-math` / + * `-Ofast` would reinstate the bug. The repo's -O3 builds are safe. * * pVect1 is storage (SQ8): [uint8_t values (dim)] [min_val] [delta] [x_sum] [x_sum_squares] * pVect2 is query (FP32): [float values (dim)] [y_sum] [y_sum_squares] @@ -45,12 +44,9 @@ float SQ8_FP32_L2Sqr(const void *pVect1v, const void *pVect2v, size_t dimension) float res = 0; for (size_t i = 0; i < dimension; i++) { - // diff = dequant(x_i) - y_i = delta * q_i + (min_val - y_i). Order matters: min_val and - // y_i are both large and close in magnitude (Sterbenz's lemma makes their subtraction - // exact), so computing that first and adding the small delta*q_i correction preserves the - // residual. Computing (min_val + delta*q_i) - y_i first rounds the dequantized value to - // FP32 at the large offset's precision, discarding the residual before the subtraction - // ever happens -- silently reintroducing the cancellation this kernel exists to avoid. + // Order matters: min_val - y_i is exact (Sterbenz) since both are large and close, so + // adding the small delta*q_i correction afterward preserves the residual. Computing + // (min_val + delta*q_i) - y_i instead rounds it away at the large offset's precision. float diff = delta * static_cast(pVect1[i]) + (min_val - pVect2[i]); res += diff * diff; } diff --git a/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h index 8c2287476..5de7f6a34 100644 --- a/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h @@ -14,51 +14,37 @@ using sq8 = vecsim_types::sq8; /* - * Asymmetric SQ8 L2 squared distance computed via direct residual accumulation: + * Asymmetric SQ8-FP32 L2 squared distance via direct residual accumulation: * - * ||x - y||² = Σ(dequant(x_i) - y_i)² - * where dequant(x_i) = min_val + delta * q_i + * ||x - y||² = Σ(dequant(x_i) - y_i)², where dequant(x_i) = min_val + delta * q_i * - * This avoids the algebraic-identity/cancellation approach, which catastrophically cancels in - * FP32 when x and y share a large common offset relative to their spread. - * - * The subtract is fused into the FMA (diff = fma(delta, q, min - y)) using the explicit - * vfmaq_f32 intrinsic rather than computed separately (and rather than relying on - * autovectorizing vmlaq_f32-style code), which matters for performance on ARM. + * Not the ||x||² + ||y||² - 2*IP identity, which cancels catastrophically in FP32 when x and y + * share a large common offset relative to their spread (MOD-17526). */ -// Helper: compute Σ(diff_i²) for 4 elements, where diff_i = dequant(x_i) - y_i. -// pVect1 = SQ8 storage (quantized values), pVect2 = FP32 query. -// min_val_vec/delta_vec are broadcast scalars from the stored vector's metadata. +// 4 elements of Σ(diff_i²). Used for the tail; the main loop uses the 16-element form below. static inline void L2StepSQ8_FP32_NEON(const uint8_t *&pVect1, const float *&pVect2, float32x4_t &sum, float32x4_t min_val_vec, float32x4_t delta_vec) { - // Load 4 uint8 elements and convert to float uint8x8_t v1_u8 = vld1_u8(pVect1); pVect1 += 4; uint32x4_t v1_u32 = vmovl_u16(vget_low_u16(vmovl_u8(v1_u8))); float32x4_t v1_f = vcvtq_f32_u32(v1_u32); - // Load 4 float elements from query float32x4_t v2 = vld1q_f32(pVect2); pVect2 += 4; - // min - y computed once per lane, then fuse the dequantize-and-subtract into a single FMA: - // diff = fma(delta, q, min - y). Uses the explicit vfmaq_f32 intrinsic, not vmlaq_f32. + // Explicit vfmaq_f32, not vmlaq_f32: keeping min - y first is what preserves the residual. float32x4_t min_minus_y = vsubq_f32(min_val_vec, v2); float32x4_t diff = vfmaq_f32(min_minus_y, delta_vec, v1_f); sum = vfmaq_f32(sum, diff, diff); } -// Helper: same arithmetic as above, but for 16 elements off a single 16-byte load. -// -// The 4-element helper reads 8 bytes (vld1_u8) and consumes only the low 4 (vget_low_u16), -// then advances the pointer by 4, so back-to-back calls re-read half of what they just loaded -// and defeat load pairing. Widening one load into all four float32x4_t groups removes that. -// Per-lane arithmetic and accumulator assignment are unchanged (group g still lands in sum, -// still fma(delta, q, min - y)), so results stay bit-identical to the 4-element path. +// 16 elements off a single 16-byte load. The 4-element form reads 8 bytes and uses only 4, so +// back-to-back calls re-read half of every load. Per-lane arithmetic and accumulator mapping are +// unchanged, so results match the 4-element path exactly. static inline void L2Step16SQ8_FP32_NEON(const uint8_t *&pVect1, const float *&pVect2, float32x4_t &sum0, float32x4_t &sum1, float32x4_t &sum2, float32x4_t &sum3, float32x4_t min_val_vec, diff --git a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h index a461c79a6..2cf5cec6c 100644 --- a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h @@ -66,17 +66,9 @@ float SQ8_FP32_L2SqrSIMD16_SSE4(const void *pVect1v, const void *pVect2v, size_t __m128 sum2 = _mm_setzero_ps(); __m128 sum3 = _mm_setzero_ps(); - // Process residual elements first (1-3 elements). - // - // Both operands are loaded at full width and the lanes past the residual are masked off, - // rather than staged through the stack. The previous form stored scalars into two 16-byte - // stack arrays and immediately reloaded them with _mm_load_ps; a 16-byte load cannot be - // store-to-load forwarded from four narrower stores, so it stalls. That showed up in the - // residual benchmark sweep as a fixed ~8-9 ns penalty on every dim where residual % 4 != 0. - // - // The wide loads are in bounds because this kernel is only reachable at dim >= 8 (the x86 - // chooser floors SIMD there), so both operands have at least 8 elements ahead of offset 0, - // and the metadata trailing each blob keeps even a 16-byte read inside the allocation. + // Residual elements (1-3), loaded at full width with the lanes past the residual masked + // off. Staging them through stack arrays instead costs a store-to-load stall. In bounds + // because the x86 chooser never reaches this kernel below dim 8. if constexpr (residual % 4) { constexpr unsigned char r = residual % 4; @@ -87,8 +79,7 @@ float SQ8_FP32_L2SqrSIMD16_SSE4(const void *pVect1v, const void *pVect2v, size_t __m128 min_minus_y = _mm_sub_ps(min_val, v2); __m128 diff = _mm_add_ps(_mm_mul_ps(delta, v1_f), min_minus_y); - // Lanes >= r hold elements the main loop will process; zero them so this step adds - // nothing for them. r is a compile-time value, so the mask is a constant. + // Lanes >= r hold elements the main loop will process; zero their contribution. const __m128 lane_mask = _mm_castsi128_ps( _mm_set_epi32(r > 3 ? -1 : 0, r > 2 ? -1 : 0, r > 1 ? -1 : 0, r > 0 ? -1 : 0)); diff = _mm_and_ps(diff, lane_mask); diff --git a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h index 2e6950363..2b891adae 100644 --- a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h @@ -20,43 +20,29 @@ namespace { using sq8 = vecsim_types::sq8; /* - * Asymmetric SQ8-FP32 L2 squared distance computed via direct residual accumulation: + * Asymmetric SQ8-FP32 L2 squared distance via direct residual accumulation: * - * ||x - y||² = Σ(dequant(x_i) - y_i)² - * where dequant(x_i) = min_val + delta * q_i + * ||x - y||² = Σ(dequant(x_i) - y_i)², where dequant(x_i) = min_val + delta * q_i * - * This avoids the algebraic-identity/cancellation approach, which catastrophically cancels in - * FP32 when x and y share a large common offset relative to their spread. - * - * The subtract is fused into the multiply-add (diff = min_minus_y + delta*q) via svmla_f32_x, - * which matters for performance. + * Not the ||x||² + ||y||² - 2*IP identity, which cancels catastrophically in FP32 when x and y + * share a large common offset relative to their spread (MOD-17526). */ -// Helper: compute Σ(diff_i²) for one SVE vector width, where diff_i = dequant(x_i) - y_i. -// pVect1 = SQ8 storage (quantized values), pVect2 = FP32 query. -// min_val_vec/delta_vec are broadcast scalars from the stored vector's metadata. -static inline void L2StepSQ8_FP32_SVE(const uint8_t *pVect1, const float *pVect2, size_t &offset, - svfloat32_t &sum, const size_t chunk, svfloat32_t min_val_vec, +// One SVE vector width of Σ(diff_i²). pVect1 = SQ8 storage, pVect2 = FP32 query. +static inline void L2StepSQ8_FP32_SVE(const uint8_t *pVect1, const float *pVect2, size_t offset, + svfloat32_t &sum, svfloat32_t min_val_vec, svfloat32_t delta_vec) { svbool_t pg = svptrue_b32(); - // Load uint8 elements and zero-extend to uint32 svuint32_t v1_u32 = svld1ub_u32(pg, pVect1 + offset); - - // Convert uint32 to float32 svfloat32_t v1_f = svcvt_f32_u32_x(pg, v1_u32); - - // Load float elements from query svfloat32_t v2 = svld1_f32(pg, pVect2 + offset); - // min - y computed once per lane, then fuse the dequantize-and-subtract: - // diff = min_minus_y + delta*q, via a single fused multiply-add. + // Subtract fused into the multiply-add; keeping min - y first is what preserves the residual. svfloat32_t min_minus_y = svsub_f32_x(pg, min_val_vec, v2); svfloat32_t diff = svmla_f32_x(pg, min_minus_y, delta_vec, v1_f); sum = svmla_f32_x(pg, sum, diff, diff); - - offset += chunk; } // pVect1v = SQ8 storage, pVect2v = FP32 query @@ -84,55 +70,35 @@ float SQ8_FP32_L2SqrSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t di svfloat32_t sum2 = svdup_f32(0.0f); svfloat32_t sum3 = svdup_f32(0.0f); - // Full-width groups first, predicated tail last, and every bound compared rather than - // divided. - // - // `chunk` is a runtime value (svcntw), so `dimension % chunk` and `/ chunk_size` compile to - // real `udiv` instructions -- ~12-20 cycles each and not pipelined, on a function called - // thousands of times per query. They are also redundant: CHOOSE_SVE_IMPLEMENTATION already - // divides once, at chooser time, and hands the results down as `partial_chunk` and - // `additional_steps`. Doing the full vectors first additionally keeps every unpredicated - // load at a multiple of the vector length; a leading partial chunk would push all of them - // to a non-VL-multiple offset. - // - // This is why the shape here differs from the prefix-first sibling SQ8 SVE kernels. Given - // dimension = k*chunk + r (r = dimension % chunk, so r > 0 exactly when partial_chunk), - // the loop below runs floor(k/4) times, leaving (k % 4) == additional_steps full vectors - // plus r tail elements. - // - // Measured on Graviton4 (Neoverse-V2, 128-bit VL), median of 9, cv <= 0.33%, against the - // prefix-first shape: SVE2 is faster everywhere (dim 1024 160 -> 141ns, dim 513 81.9 -> - // 79.4ns), and SVE is faster where dimension is a multiple of the vector length (dim 1024 - // 164 -> 153ns) but ~8% slower at dims that leave a tail (dim 513 78.3 -> 84.6ns). Kept - // because SVE2 is the tier the chooser prefers when both are present, and because real - // embedding dims (128/256/384/512/768/1024/1536) all divide evenly at 128- and 256-bit VL, - // landing on the faster path. Feeding the tail through its own accumulator instead of - // svmla_f32_m was also tried and measured slightly worse (dim 513 85.7ns). - // - // Note `dimension & (chunk - 1)` is NOT a valid substitute for the modulo: SVE permits any - // vector length that is a multiple of 128 bits, not only powers of two. + // Full vectors first, predicated tail last, bounds compared rather than divided: `chunk` is + // runtime (svcntw), so `dimension % chunk` would emit a real udiv, and the chooser already + // divided once to produce `partial_chunk`/`additional_steps`. Deliberately unlike the + // prefix-first sibling SQ8 SVE kernels; see MOD-17526 for the measured trade-off. + // A bitmask cannot replace the modulo here: SVE vector length is a multiple of 128 bits, + // not necessarily a power of two. const size_t chunk_size = 4 * chunk; - while (offset + chunk_size <= dimension) { - L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum0, chunk, min_val_vec, delta_vec); - L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum1, chunk, min_val_vec, delta_vec); - L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum2, chunk, min_val_vec, delta_vec); - L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum3, chunk, min_val_vec, delta_vec); + for (; offset + chunk_size <= dimension; offset += chunk_size) { + L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum0, min_val_vec, delta_vec); + L2StepSQ8_FP32_SVE(pVect1, pVect2, offset + chunk, sum1, min_val_vec, delta_vec); + L2StepSQ8_FP32_SVE(pVect1, pVect2, offset + 2 * chunk, sum2, min_val_vec, delta_vec); + L2StepSQ8_FP32_SVE(pVect1, pVect2, offset + 3 * chunk, sum3, min_val_vec, delta_vec); } // Handle remaining full-width steps (0-3), resolved at compile time. if constexpr (additional_steps > 0) { - L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum0, chunk, min_val_vec, delta_vec); + L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum0, min_val_vec, delta_vec); + offset += chunk; } if constexpr (additional_steps > 1) { - L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum1, chunk, min_val_vec, delta_vec); + L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum1, min_val_vec, delta_vec); + offset += chunk; } if constexpr (additional_steps > 2) { - L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum2, chunk, min_val_vec, delta_vec); + L2StepSQ8_FP32_SVE(pVect1, pVect2, offset, sum2, min_val_vec, delta_vec); + offset += chunk; } - // Predicated tail for the final partial vector. svwhilelt derives the predicate from the - // current offset against `dimension`, so no residual arithmetic is needed, and inactive - // lanes of the `_z` (zeroing) forms contribute 0 to the accumulator. + // Predicated tail for the final partial vector. if constexpr (partial_chunk) { svbool_t pg_tail = svwhilelt_b32(static_cast(offset), static_cast(dimension)); @@ -144,11 +110,9 @@ float SQ8_FP32_L2SqrSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t di svfloat32_t min_minus_y = svsub_f32_z(pg_tail, min_val_vec, v2); svfloat32_t diff = svmla_f32_z(pg_tail, min_minus_y, delta_vec, v1_f); - // Merging (_m), not zeroing (_z): sum3 already holds full-vector partial sums from the - // main loop, and _z would zero every lane outside pg_tail, silently dropping them. _z is - // only safe on a freshly zeroed accumulator, which is why the pre-restructure kernel - // could use it here -- back then this block ran first, against an all-zero sum0. - // The temporaries above stay _z because they are fresh values, not accumulators. + // Merging (_m), not zeroing: sum3 already holds partial sums from the main loop, and _z + // would zero every lane outside pg_tail. Only the accumulator needs this; the + // temporaries above are fresh values. sum3 = svmla_f32_m(pg_tail, sum3, diff, diff); } diff --git a/tests/unit/test_spaces.cpp b/tests/unit/test_spaces.cpp index 3bdfcf155..66d82fe17 100644 --- a/tests/unit/test_spaces.cpp +++ b/tests/unit/test_spaces.cpp @@ -411,28 +411,18 @@ TEST_F(SpacesTest, SQ8_FP32_odd_dim_unaligned_metadata_test) { } /* ==================== MOD-17526 regression tests ==================== - * SQ8_FP32_L2Sqr and its SIMD variants used to compute L2^2 via the algebraic identity - * ||x||^2 + ||y||^2 - 2*IP(x, y), reading precomputed sums from the SQ8 blob metadata. That - * identity catastrophically cancels in FP32 when the two vectors share a large common offset - * relative to their spread (e.g. x=[100000,100008], y=[100000,100000]: true L2^2=64, the old - * kernel could return 0, a negative value, or otherwise unbounded garbage because it computed a - * huge positive value minus a nearly-equal huge positive value in single precision). + * The SQ8-FP32 L2 kernels used to compute L2^2 via the ||x||^2 + ||y||^2 - 2*IP identity, which + * cancels catastrophically in FP32 when both vectors share a large common offset. * - * The tests below use an independent double-precision reference computed by manually - * dequantizing the SQ8 storage and accumulating squared residuals -- NOT the algebraic identity, - * and NOT test_utils::SQ8_FP32_NotOptimized_L2Sqr or any other kernel variant -- so they would - * have failed against the old kernel and remain meaningful now that both are fixed. They also - * use *relative* tolerances (scaled to the true distance) rather than the ~0.01 absolute - * tolerance used elsewhere in this file, since an absolute tolerance is exactly what let the old - * bug hide: the large-offset case produces errors many orders of magnitude larger than 0.01. + * These check against an independent double-precision reference rather than another kernel + * (pre-fix, scalar and SIMD were wrong together), and use relative rather than absolute + * tolerances, since a ~0.01 absolute bound is what let the bug hide. */ namespace { -// Manually dequantizes the SQ8 storage and accumulates the squared residual against the FP32 -// query in double precision. Deliberately independent of both the production kernels and of -// test_utils::SQ8_FP32_NotOptimized_L2Sqr (which accumulates in float) so it serves as ground -// truth rather than another instance of the same computation. +// Ground truth: dequantizes the storage by hand and accumulates in double. Independent of the +// production kernels and of test_utils::SQ8_FP32_NotOptimized_L2Sqr, which accumulates in float. double SQ8_FP32_L2Sqr_DoubleReference(const uint8_t *storage, const float *query, size_t dim) { const float min_val = load_unaligned(storage + dim + sq8::MIN_VAL * sizeof(float)); const float delta = load_unaligned(storage + dim + sq8::DELTA * sizeof(float)); @@ -446,10 +436,8 @@ double SQ8_FP32_L2Sqr_DoubleReference(const uint8_t *storage, const float *query return res; } -// Runs the scalar kernel, the runtime-dispatched kernel, and every SIMD variant compiled into -// this binary against `storage`/`query`, asserting each is within `rel_tol` of `expected` -// (relative to max(expected, 1.0), so a tiny true distance isn't swamped by an absolute -// epsilon and a huge one isn't judged by an absolute epsilon that's too tight). +// Checks the scalar kernel, the dispatched kernel, and every SIMD variant in this binary against +// `expected`, within `rel_tol` scaled by max(expected, 1.0). void ExpectSQ8_FP32_L2SqrNear(const uint8_t *storage, const float *query, size_t dim, double expected, double rel_tol, const std::string &context) { const double scale = std::max(expected, 1.0); @@ -463,13 +451,10 @@ void ExpectSQ8_FP32_L2SqrNear(const uint8_t *storage, const float *query, size_t auto optimization = getCpuOptimizationFeatures(); - // The x86 kernels are only ever reached through the chooser at dim >= 8 (L2_space.cpp: - // "Optimizations assume at least 8 elements (see the residual handling in the kernels)"). - // Calling them directly below that floor tests a configuration production never produces, - // and the AVX2 residual path reads a full 32-byte vector via my_mm256_maskz_loadu_ps - // (AVX_utils.h) -- an unconditional _mm256_loadu_ps -- which overreads a dim<8 query blob - // and trips ASAN. aarch64 has no such floor (the chooser returns NEON/SVE at any dim, and - // the NEON tail uses per-lane loads), so those tiers stay ungated here to match production. + // The chooser only reaches the x86 kernels at dim >= 8, and the AVX2 residual path reads a + // full 32-byte vector (my_mm256_maskz_loadu_ps is an unconditional load), which overreads a + // smaller query blob. aarch64 has no such floor, so those tiers stay ungated to match + // production. const bool x86_simd_reachable = dim >= 8; #ifdef OPT_AVX512_F_BW_VL_VNNI if (x86_simd_reachable && optimization.avx512f && optimization.avx512bw && @@ -510,12 +495,10 @@ void ExpectSQ8_FP32_L2SqrNear(const uint8_t *storage, const float *query, size_t #endif } -// Builds an SQ8 storage blob and an FP32 query for L2 from an explicit pair of float vectors -// (query_values, storage_values). Kept separate from generation so callers can build one base -// pair and derive a shifted pair by literally adding a constant to the same values -- re-sampling -// std::uniform_real_distribution with a shifted [min,max] range and the same seed does NOT -// reproduce "the same vector plus a constant" (it's an independently-shaped draw over the new -// range), which is not what a translation-invariance test needs. +// Builds the blobs from explicit float vectors, kept separate from generation so a caller can +// derive a shifted pair by adding a constant to the same values. Re-sampling +// std::uniform_real_distribution over a shifted [min,max] with the same seed does NOT give "the +// same vector plus a constant", which is what a translation-invariance test needs. struct SQ8_FP32_L2_TestVectors { std::vector query; // [float values (dim)] [sum] [sum_squares] std::vector storage; // [uint8_t values (dim)] [min] [delta] [sum] [sum_squares] @@ -537,9 +520,7 @@ BuildSQ8_FP32_L2_TestVectorsFromValues(const std::vector &query_values, return v; } -// Draws one base (query, storage) float pair from N(0, spread) and returns it alongside the -// same pair with `offset` added to every element -- a real shift of the identical values, not -// two independently-seeded draws over different ranges. +// One base (query, storage) pair, plus the same pair with `offset` added to every element. std::pair BuildSQ8_FP32_L2_ShiftedPair(size_t dim, float spread, float offset, int seed) { std::vector query_values(dim), storage_values(dim); @@ -566,9 +547,8 @@ TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_TranslationInvariance) { const float large_offset = 100000.0f; // offset/spread == 100000, far past the ~4000 threshold // where the old identity started to catastrophically // cancel in FP32. - // Dims deliberately mix exact multiples of the SIMD chunk widths with awkward remainders - // (129, 145, 513, 527), so the large-offset cancellation case also runs through the - // residual/masked-lane tails, not just the aligned main loops. + // Mixes exact multiples of the SIMD chunk widths with awkward remainders, so the + // large-offset case also runs through the residual/masked-lane tails. for (const size_t dim : {128UL, 129UL, 145UL, 512UL, 513UL, 527UL, 768UL}) { auto [no_offset, with_offset] = BuildSQ8_FP32_L2_ShiftedPair(dim, spread, large_offset, 1234); @@ -578,11 +558,9 @@ TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_TranslationInvariance) { const double expected_with_offset = SQ8_FP32_L2Sqr_DoubleReference( with_offset.storage.data(), with_offset.query.data(), dim); - // Sanity check on the reference itself: the two builds only differ by a shared additive - // constant on both vectors, so their true L2^2 must match closely. A small tolerance - // (rather than exact equality) accounts for the shifted float32 inputs themselves losing - // a bit of precision at the larger magnitude (float32 spacing near 1e5 is coarser than - // near 0), which perturbs which byte a value quantizes to right at a code boundary. + // The reference itself must be translation-invariant. Not exact equality: float32 + // spacing near 1e5 is coarser than near 0, which can shift a value to an adjacent + // quantization code. ASSERT_NEAR(expected_no_offset, expected_with_offset, 1e-2 * std::max(expected_no_offset, 1.0)) << "dim " << dim << ": independent reference is not translation-invariant"; @@ -590,11 +568,8 @@ TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_TranslationInvariance) { ExpectSQ8_FP32_L2SqrNear(no_offset.storage.data(), no_offset.query.data(), dim, expected_no_offset, 1e-3, "dim " + std::to_string(dim) + ", C=0"); - // The direct-diff kernel computes `diff = delta*q + (min - y)`. `min` and `y` both carry - // the large shared offset and are within a factor of 2 of each other, so by Sterbenz's - // lemma `min - y` is computed *exactly* in FP32 (no rounding at all) -- it's the small - // `delta*q` correction, not the offset, that's added afterward. So there's no large-scale - // cancellation left to tolerate here; the tolerance stays as tight as the no-offset case. + // `min - y` is exact in FP32 by Sterbenz (both large, within a factor of 2), so no + // large-scale cancellation remains and the tolerance stays as tight as the C=0 case. ExpectSQ8_FP32_L2SqrNear( with_offset.storage.data(), with_offset.query.data(), dim, expected_with_offset, 1e-3, "dim " + std::to_string(dim) + ", C=" + std::to_string(large_offset)); @@ -623,17 +598,11 @@ TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_TicketReproShape) { } TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_ScalarAssociativityRegression) { - // Pins a scalar-only regression found in review: `diff = (min_val + delta*q) - y` and - // `diff = delta*q + (min_val - y)` are the same expression algebraically, but not in FP32 -- - // the first rounds the dequantized value to the large offset's precision *before* subtracting - // y, silently discarding the residual the direct-diff kernel exists to preserve; the second - // computes the (Sterbenz-exact) min-y cancellation first and adds the small correction after. - // - // Concrete numbers where this bites: min=100000, delta=8/255, q=1 dequantizes to - // 100000.031372549..., which FP32 rounds to 100000.03125 (float32 ULP near 1e5 is 1/64). - // With y=100000.03125 chosen to land on exactly that rounded value, the buggy order gives - // diff=0 (100000.03125 - 100000.03125), while the correct order gives the true residual - // delta*1 + (100000 - 100000.03125) ~= 0.0001225, matching the double reference. + // Pins a scalar-only ordering bug: `(min_val + delta*q) - y` rounds the dequantized value at + // the large offset's precision before subtracting, discarding the residual, while + // `delta*q + (min_val - y)` keeps it. With min=100000, delta=8/255, q=1 the dequantized value + // is 100000.031372549..., which FP32 rounds to exactly y=100000.03125, so the buggy order + // yields diff=0 instead of ~0.0001225. const size_t dim = 8; // storage: min=100000 (index 0), max=100008 (index 7) => delta = 8/255, so index 1 // quantizes to q=1 (dequantizes to 100000 + 8/255 = 100000.031372549...). @@ -655,9 +624,7 @@ TEST_F(SpacesTest, SQ8_FP32_L2Sqr_MOD17526_ScalarAssociativityRegression) { const double expected = SQ8_FP32_L2Sqr_DoubleReference(storage.data(), query.data(), dim); ASSERT_GT(expected, 0.0) << "test construction should produce a nonzero true distance"; - // Absolute tolerance, not relative: the true value is tiny (~1.5e-8), and the bug this - // pins is "kernel returns exactly 0 instead of a small nonzero value", not a rounding-scale - // discrepancy -- a relative check against a near-zero expected value isn't meaningful here. + // Absolute, not relative: the true value is ~1.5e-8 and the bug returns exactly 0. EXPECT_NEAR(SQ8_FP32_L2Sqr(storage.data(), query.data(), dim), expected, 1e-9) << "scalar kernel: expected=" << expected; } From b57df308ef37ca7a56e9b6979cd31c9f80c852fb Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 2 Sep 2026 12:57:25 +0300 Subject: [PATCH 09/10] style(sq8): give the x86 SQ8-FP32 L2 kernels the same for-loop shape The three x86 kernels drove their main loop off a pEnd1 sentinel while their NEON and SVE counterparts count chunks, so the family this PR touches read three different ways. They now all count chunks. Scoped to SQ8-FP32 L2: the sentinel is still the shape in ~32 other x86 kernels, and churning those is unrelated to this fix. Costs nothing: the residual blocks consume dimension % 32, leaving exactly 32 * floor(dimension / 32), and `dimension / 32` divides by a literal, so it is a shift rather than the runtime udiv the SVE version had to avoid. Accumulator counts deliberately unchanged -- AVX2 and AVX2+FMA already use four, and AVX512 stays at two because that was measured faster than four (95.0ns vs 95.9ns at dim 1536). 595/595 SQ8 tests, and 595/595 under ASAN, on x86. Benchmarks land within run-to-run noise, as expected for an unchanged instruction sequence. --- src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h | 4 ++-- src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h | 4 ++-- src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h index 285bb2ab7..2776b9a52 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h @@ -55,7 +55,6 @@ template // 0..31 float SQ8_FP32_L2SqrSIMD16_AVX2_FMA(const void *pVect1v, const void *pVect2v, size_t dimension) { const uint8_t *pVect1 = static_cast(pVect1v); // SQ8 storage const float *pVect2 = static_cast(pVect2v); // FP32 query - const uint8_t *pEnd1 = pVect1 + dimension; // Get quantization parameters from stored vector (after quantized data) const uint8_t *pVect1Base = static_cast(pVect1v); @@ -112,7 +111,8 @@ float SQ8_FP32_L2SqrSIMD16_AVX2_FMA(const void *pVect1v, const void *pVect2v, si // We dealt with the residual part. We are left with some multiple of 32 elements. // In each iteration we calculate 32 elements = 4 chunks of 8. The loop may run zero times // (dim can be as small as 8). - while (pVect1 < pEnd1) { + const size_t num_of_chunks = dimension / 32; + for (size_t i = 0; i < num_of_chunks; i++) { L2StepSQ8_FP32_FMA(pVect1, pVect2, sum0, min_val_vec, delta_vec); L2StepSQ8_FP32_FMA(pVect1, pVect2, sum1, min_val_vec, delta_vec); L2StepSQ8_FP32_FMA(pVect1, pVect2, sum2, min_val_vec, delta_vec); diff --git a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h index f3f1a31f8..7a8152f27 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h @@ -51,7 +51,6 @@ template // 0..31 float SQ8_FP32_L2SqrSIMD16_AVX2(const void *pVect1v, const void *pVect2v, size_t dimension) { const uint8_t *pVect1 = static_cast(pVect1v); // SQ8 storage const float *pVect2 = static_cast(pVect2v); // FP32 query - const uint8_t *pEnd1 = pVect1 + dimension; // Get quantization parameters from stored vector (after quantized data) const uint8_t *pVect1Base = static_cast(pVect1v); @@ -108,7 +107,8 @@ float SQ8_FP32_L2SqrSIMD16_AVX2(const void *pVect1v, const void *pVect2v, size_t // We dealt with the residual part. We are left with some multiple of 32 elements. // In each iteration we calculate 32 elements = 4 chunks of 8. The loop may run zero times // (dim can be as small as 8). - while (pVect1 < pEnd1) { + const size_t num_of_chunks = dimension / 32; + for (size_t i = 0; i < num_of_chunks; i++) { L2StepSQ8_FP32_AVX2(pVect1, pVect2, sum0, min_val_vec, delta_vec); L2StepSQ8_FP32_AVX2(pVect1, pVect2, sum1, min_val_vec, delta_vec); L2StepSQ8_FP32_AVX2(pVect1, pVect2, sum2, min_val_vec, delta_vec); diff --git a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h index 0c321a228..da5b99df0 100644 --- a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h @@ -56,7 +56,6 @@ float SQ8_FP32_L2SqrSIMD16_AVX512F_BW_VL_VNNI(const void *pVec1v, const void *pV size_t dimension) { const uint8_t *pVec1 = static_cast(pVec1v); // SQ8 storage const float *pVec2 = static_cast(pVec2v); // FP32 query - const uint8_t *pEnd1 = pVec1 + dimension; // Get quantization parameters from stored vector (after quantized data) const uint8_t *pVec1Base = static_cast(pVec1v); @@ -105,7 +104,8 @@ float SQ8_FP32_L2SqrSIMD16_AVX512F_BW_VL_VNNI(const void *pVec1v, const void *pV // We dealt with the residual part. We are left with some multiple of 32 elements. // In each iteration we calculate 32 elements = 2 chunks of 16. The loop may run zero times // (dim can be as small as 8). - while (pVec1 < pEnd1) { + const size_t num_of_chunks = dimension / 32; + for (size_t i = 0; i < num_of_chunks; i++) { L2StepSQ8_FP32_AVX512(pVec1, pVec2, sum0, min_val_vec, delta_vec); L2StepSQ8_FP32_AVX512(pVec1, pVec2, sum1, min_val_vec, delta_vec); } From 63a23092799f21befcb5f95fe9293bebdcb26f42 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 2 Sep 2026 13:09:39 +0300 Subject: [PATCH 10/10] Revert "style(sq8): give the x86 SQ8-FP32 L2 kernels the same for-loop shape" This reverts commit b57df308ef37ca7a56e9b6979cd31c9f80c852fb. --- src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h | 4 ++-- src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h | 4 ++-- src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h index 2776b9a52..285bb2ab7 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h @@ -55,6 +55,7 @@ template // 0..31 float SQ8_FP32_L2SqrSIMD16_AVX2_FMA(const void *pVect1v, const void *pVect2v, size_t dimension) { const uint8_t *pVect1 = static_cast(pVect1v); // SQ8 storage const float *pVect2 = static_cast(pVect2v); // FP32 query + const uint8_t *pEnd1 = pVect1 + dimension; // Get quantization parameters from stored vector (after quantized data) const uint8_t *pVect1Base = static_cast(pVect1v); @@ -111,8 +112,7 @@ float SQ8_FP32_L2SqrSIMD16_AVX2_FMA(const void *pVect1v, const void *pVect2v, si // We dealt with the residual part. We are left with some multiple of 32 elements. // In each iteration we calculate 32 elements = 4 chunks of 8. The loop may run zero times // (dim can be as small as 8). - const size_t num_of_chunks = dimension / 32; - for (size_t i = 0; i < num_of_chunks; i++) { + while (pVect1 < pEnd1) { L2StepSQ8_FP32_FMA(pVect1, pVect2, sum0, min_val_vec, delta_vec); L2StepSQ8_FP32_FMA(pVect1, pVect2, sum1, min_val_vec, delta_vec); L2StepSQ8_FP32_FMA(pVect1, pVect2, sum2, min_val_vec, delta_vec); diff --git a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h index 7a8152f27..f3f1a31f8 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h @@ -51,6 +51,7 @@ template // 0..31 float SQ8_FP32_L2SqrSIMD16_AVX2(const void *pVect1v, const void *pVect2v, size_t dimension) { const uint8_t *pVect1 = static_cast(pVect1v); // SQ8 storage const float *pVect2 = static_cast(pVect2v); // FP32 query + const uint8_t *pEnd1 = pVect1 + dimension; // Get quantization parameters from stored vector (after quantized data) const uint8_t *pVect1Base = static_cast(pVect1v); @@ -107,8 +108,7 @@ float SQ8_FP32_L2SqrSIMD16_AVX2(const void *pVect1v, const void *pVect2v, size_t // We dealt with the residual part. We are left with some multiple of 32 elements. // In each iteration we calculate 32 elements = 4 chunks of 8. The loop may run zero times // (dim can be as small as 8). - const size_t num_of_chunks = dimension / 32; - for (size_t i = 0; i < num_of_chunks; i++) { + while (pVect1 < pEnd1) { L2StepSQ8_FP32_AVX2(pVect1, pVect2, sum0, min_val_vec, delta_vec); L2StepSQ8_FP32_AVX2(pVect1, pVect2, sum1, min_val_vec, delta_vec); L2StepSQ8_FP32_AVX2(pVect1, pVect2, sum2, min_val_vec, delta_vec); diff --git a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h index da5b99df0..0c321a228 100644 --- a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h @@ -56,6 +56,7 @@ float SQ8_FP32_L2SqrSIMD16_AVX512F_BW_VL_VNNI(const void *pVec1v, const void *pV size_t dimension) { const uint8_t *pVec1 = static_cast(pVec1v); // SQ8 storage const float *pVec2 = static_cast(pVec2v); // FP32 query + const uint8_t *pEnd1 = pVec1 + dimension; // Get quantization parameters from stored vector (after quantized data) const uint8_t *pVec1Base = static_cast(pVec1v); @@ -104,8 +105,7 @@ float SQ8_FP32_L2SqrSIMD16_AVX512F_BW_VL_VNNI(const void *pVec1v, const void *pV // We dealt with the residual part. We are left with some multiple of 32 elements. // In each iteration we calculate 32 elements = 2 chunks of 16. The loop may run zero times // (dim can be as small as 8). - const size_t num_of_chunks = dimension / 32; - for (size_t i = 0; i < num_of_chunks; i++) { + while (pVec1 < pEnd1) { L2StepSQ8_FP32_AVX512(pVec1, pVec2, sum0, min_val_vec, delta_vec); L2StepSQ8_FP32_AVX512(pVec1, pVec2, sum1, min_val_vec, delta_vec); }