From 6a4803a495c844cb445cf8741f706698fdd9eeba Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Mon, 10 Aug 2026 16:57:41 +0300 Subject: [PATCH 01/17] Make SQ8 metadata exact and fix the symmetric L2 and IP formulations WIP: source complete, test fixtures not yet updated, not yet compiled. Stored SQ8 metadata now describes the reconstruction rather than the input, and the sums are integers rather than FP32. * sq8.h: SUM and SUM_SQUARES become Q_SUM and Q_SUM_SQUARES, uint32 sums over the quantized bytes. Same slot count and blob size. Adds unaligned readers and four exact-derivation helpers (reconstructed_sum, reconstructed_sum_squares, reconstructed_ip, reconstructed_l2_sqr) so the layout has one interpretation. MAX_EXACT_DIM records the dimension past which sum(a^2) leaves uint32. * preprocessors.h: accumulates the quantized bytes in integers and packs mixed FP32/uint32 metadata. * The symmetric IP kernels previously combined sums taken over the original input with a dot product of the reconstructions, mixing two different vectors into one formula. All five now derive every term from the integer sums. * The symmetric L2 kernels used ||x||^2 + ||y||^2 - 2*IP in FP32, which cancels away the answer when the vectors share a large offset: stored [100000, 100008] against [100000, 100000] returned 0 where the truth is 64. They now use the expanded difference, which keeps the offset in a (min1 - min2) term formed once and exactly, and combine in double. Measured relative error drops to ~2e-8, including at dim 16384 where the FP32 combination was 16% off. * The asymmetric L2 kernels derive ||x||^2 exactly from the integer sums, which removes the negative distances (stored [0, 0.25, 1] against [0, 0.2501, 1] returned -0.000490427). Their remaining cancellation needs the difference taken inside the loop and is the next change in this stack. * UINT8_InnerProductImp now returns the exact integer dot on all four ISAs. NEON and SVE accumulated exactly in integer lanes and then discarded it by returning float, exact only to dimension 258; AVX512 reduced 16 int32 lanes into an int, which overflows from dimension 33,027. Callers that wanted a float now convert explicitly, and one that computed 1 - imp() as int is now float arithmetic. This is the widening asked for separately, needed here because the L2 expansion depends on the dot being exact. Verified so far: clang-format clean; -Wall -Werror -fsyntax-only clean on both scalar TUs and on 17 of 19 x86 SQ8 kernels. The other two are not standalone compilable in an ad-hoc TU, on main as much as here, so they are covered only by a real build. Nothing has been executed yet. Co-Authored-By: Claude Opus 5 --- src/VecSim/spaces/IP/IP.cpp | 28 ++--- .../spaces/IP/IP_AVX512F_BW_VL_VNNI_SQ8_SQ8.h | 27 ++--- .../spaces/IP/IP_AVX512F_BW_VL_VNNI_UINT8.h | 15 ++- .../spaces/IP/IP_NEON_DOTPROD_SQ8_SQ8.h | 27 ++--- src/VecSim/spaces/IP/IP_NEON_DOTPROD_UINT8.h | 12 +-- src/VecSim/spaces/IP/IP_NEON_SQ8_SQ8.h | 28 ++--- src/VecSim/spaces/IP/IP_NEON_UINT8.h | 13 ++- src/VecSim/spaces/IP/IP_SVE_SQ8_SQ8.h | 27 ++--- src/VecSim/spaces/IP/IP_SVE_UINT8.h | 16 +-- src/VecSim/spaces/L2/L2.cpp | 49 +++++---- src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP16.h | 5 +- src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h | 6 +- src/VecSim/spaces/L2/L2_AVX2_SQ8_FP16.h | 5 +- src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h | 6 +- .../L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h | 6 +- .../spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_SQ8.h | 25 ++--- src/VecSim/spaces/L2/L2_AVX512F_SQ8_FP16.h | 5 +- .../spaces/L2/L2_NEON_DOTPROD_SQ8_SQ8.h | 25 +++-- src/VecSim/spaces/L2/L2_NEON_SQ8_FP16.h | 10 +- src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h | 6 +- src/VecSim/spaces/L2/L2_NEON_SQ8_SQ8.h | 24 ++--- src/VecSim/spaces/L2/L2_SSE4_SQ8_FP16.h | 5 +- src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h | 6 +- src/VecSim/spaces/L2/L2_SVE2_SQ8_FP16.h | 5 +- src/VecSim/spaces/L2/L2_SVE_SQ8_FP16.h | 5 +- src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h | 6 +- src/VecSim/spaces/L2/L2_SVE_SQ8_SQ8.h | 26 ++--- src/VecSim/spaces/computer/preprocessors.h | 82 ++++++++------ src/VecSim/types/sq8.h | 102 +++++++++++++++++- 29 files changed, 338 insertions(+), 264 deletions(-) diff --git a/src/VecSim/spaces/IP/IP.cpp b/src/VecSim/spaces/IP/IP.cpp index 2140c2345..6894f7c60 100644 --- a/src/VecSim/spaces/IP/IP.cpp +++ b/src/VecSim/spaces/IP/IP.cpp @@ -146,28 +146,18 @@ float SQ8_SQ8_InnerProduct_Impl(const void *pVect1v, const void *pVect2v, size_t const auto *pVect1 = static_cast(pVect1v); const auto *pVect2 = static_cast(pVect2v); - // Compute inner product of quantized values: Σ(q1[i]*q2[i]) - float product = 0; + // Inner product of the quantized bytes: Σ(a[i]*b[i]). Exact in integer arithmetic. + uint64_t q_dot = 0; for (size_t i = 0; i < dimension; i++) { - product += pVect1[i] * pVect2[i]; + q_dot += static_cast(pVect1[i]) * static_cast(pVect2[i]); } - // Metadata follows byte payloads and is not necessarily float-aligned. - const auto *params1 = pVect1 + dimension; - const float min_val1 = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); - const float delta1 = load_unaligned(params1 + sq8::DELTA * sizeof(float)); - const float sum1 = load_unaligned(params1 + sq8::SUM * sizeof(float)); - - // Get quantization parameters from pVect2 - const auto *params2 = pVect2 + dimension; - const float min_val2 = load_unaligned(params2 + sq8::MIN_VAL * sizeof(float)); - const float delta2 = load_unaligned(params2 + sq8::DELTA * sizeof(float)); - const float sum2 = load_unaligned(params2 + sq8::SUM * sizeof(float)); - - // Apply the algebraic formula using precomputed sums: - // IP = min1*sum2 + min2*sum1 + delta1*delta2*Σ(q1[i]*q2[i]) - dim*min1*min2 - return min_val1 * sum2 + min_val2 * sum1 - static_cast(dimension) * min_val1 * min_val2 + - delta1 * delta2 * product; + // Both operands describe their reconstructions, so every term is derived from the stored + // integer sums rather than from sums over the original inputs. See sq8.h. + return static_cast(sq8::reconstructed_ip( + sq8::min_val_of(pVect1, dimension), sq8::delta_of(pVect1, dimension), + sq8::q_sum_of(pVect1, dimension), sq8::min_val_of(pVect2, dimension), + sq8::delta_of(pVect2, dimension), sq8::q_sum_of(pVect2, dimension), q_dot, dimension)); } // SQ8-to-SQ8: Both vectors are uint8 quantized with precomputed sum diff --git a/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_SQ8_SQ8.h b/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_SQ8_SQ8.h index ae6f96ea2..58480dcf5 100644 --- a/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_SQ8_SQ8.h +++ b/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_SQ8_SQ8.h @@ -38,29 +38,16 @@ using sq8 = vecsim_types::sq8; // Uses UINT8_InnerProductImp for efficient dot product computation with VNNI template // 0..63 float SQ8_SQ8_InnerProductImp(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Compute raw dot product using efficient UINT8 AVX512 VNNI implementation - // UINT8_InnerProductImp uses _mm512_dpwssd_epi32 for native integer dot product - int dot_product = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); + // Exact integer dot product of the quantized bytes. + const uint64_t q_dot = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); - // Get dequantization parameters and precomputed values from the end of vectors - // Layout: [data (dim)] [min (float)] [delta (float)] [sum (float)] + // Every term describes the reconstructions, derived from the stored integer sums. See sq8.h. const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - - const auto *params1 = pVec1 + dimension; - const float min1 = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); - const float delta1 = load_unaligned(params1 + sq8::DELTA * sizeof(float)); - const float sum1 = load_unaligned(params1 + sq8::SUM * sizeof(float)); - - const auto *params2 = pVec2 + dimension; - const float min2 = load_unaligned(params2 + sq8::MIN_VAL * sizeof(float)); - const float delta2 = load_unaligned(params2 + sq8::DELTA * sizeof(float)); - const float sum2 = load_unaligned(params2 + sq8::SUM * sizeof(float)); - - // Apply the algebraic formula using precomputed sums: - // IP = min1*sum2 + min2*sum1 + δ1*δ2 * Σ(q1[i]*q2[i]) - dim*min1*min2 - return min1 * sum2 + min2 * sum1 + delta1 * delta2 * static_cast(dot_product) - - static_cast(dimension) * min1 * min2; + return static_cast(sq8::reconstructed_ip( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::min_val_of(pVec2, dimension), + sq8::delta_of(pVec2, dimension), sq8::q_sum_of(pVec2, dimension), q_dot, dimension)); } // SQ8-to-SQ8 Inner Product distance function diff --git a/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_UINT8.h b/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_UINT8.h index bd43bc901..3d95e53c9 100644 --- a/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_UINT8.h +++ b/src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_UINT8.h @@ -31,8 +31,8 @@ static inline void InnerProductStep(uint8_t *&pVect1, uint8_t *&pVect2, __m512i } template // 0..63 -static inline int UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, - size_t dimension) { +static inline uint64_t UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, + size_t dimension) { uint8_t *pVect1 = (uint8_t *)pVect1v; uint8_t *pVect2 = (uint8_t *)pVect2v; @@ -87,19 +87,24 @@ static inline int UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v } while (pVect1 < pEnd1); } - return _mm512_reduce_add_epi32(sum); + // Reduce in 64-bit. Each int32 lane is in range, but their total reaches 255*255*dim, which + // passes INT_MAX from dimension 33,027, so a 32-bit horizontal sum would wrap. Callers rely + // on this being the exact integer dot product. + const __m512i lo = _mm512_cvtepu32_epi64(_mm512_extracti64x4_epi64(sum, 0)); + const __m512i hi = _mm512_cvtepu32_epi64(_mm512_extracti64x4_epi64(sum, 1)); + return static_cast(_mm512_reduce_add_epi64(_mm512_add_epi64(lo, hi))); } template // 0..63 float UINT8_InnerProductSIMD64_AVX512F_BW_VL_VNNI(const void *pVect1v, const void *pVect2v, size_t dimension) { - return 1 - UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + return 1.0f - static_cast(UINT8_InnerProductImp(pVect1v, pVect2v, dimension)); } template // 0..63 float UINT8_CosineSIMD64_AVX512F_BW_VL_VNNI(const void *pVect1v, const void *pVect2v, size_t dimension) { - float ip = UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + float ip = static_cast(UINT8_InnerProductImp(pVect1v, pVect2v, dimension)); const float norm_v1 = load_unaligned(static_cast(pVect1v) + dimension); const float norm_v2 = load_unaligned(static_cast(pVect2v) + dimension); return 1.0f - ip / (norm_v1 * norm_v2); diff --git a/src/VecSim/spaces/IP/IP_NEON_DOTPROD_SQ8_SQ8.h b/src/VecSim/spaces/IP/IP_NEON_DOTPROD_SQ8_SQ8.h index d7f5b444e..b16b3fc99 100644 --- a/src/VecSim/spaces/IP/IP_NEON_DOTPROD_SQ8_SQ8.h +++ b/src/VecSim/spaces/IP/IP_NEON_DOTPROD_SQ8_SQ8.h @@ -39,29 +39,16 @@ using sq8 = vecsim_types::sq8; template // 0..63 float SQ8_SQ8_InnerProductSIMD64_NEON_DOTPROD_IMP(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Compute raw dot product using efficient UINT8 DOTPROD implementation - // UINT8_InnerProductImp uses vdotq_u32 for native uint8 dot product - float dot_product = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); + // Exact integer dot product of the quantized bytes. + const uint64_t q_dot = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); - // Get dequantization parameters and precomputed values from the end of vectors - // Layout: [data (dim)] [min (float)] [delta (float)] [sum (float)] + // Every term describes the reconstructions, derived from the stored integer sums. See sq8.h. const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - - const auto *params1 = pVec1 + dimension; - const float min1 = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); - const float delta1 = load_unaligned(params1 + sq8::DELTA * sizeof(float)); - const float sum1 = load_unaligned(params1 + sq8::SUM * sizeof(float)); - - const auto *params2 = pVec2 + dimension; - const float min2 = load_unaligned(params2 + sq8::MIN_VAL * sizeof(float)); - const float delta2 = load_unaligned(params2 + sq8::DELTA * sizeof(float)); - const float sum2 = load_unaligned(params2 + sq8::SUM * sizeof(float)); - - // Apply algebraic formula using precomputed sums: - // IP = min1*sum2 + min2*sum1 + δ1*δ2 * Σ(q1*q2) - dim*min1*min2 - return min1 * sum2 + min2 * sum1 + delta1 * delta2 * dot_product - - static_cast(dimension) * min1 * min2; + return static_cast(sq8::reconstructed_ip( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::min_val_of(pVec2, dimension), + sq8::delta_of(pVec2, dimension), sq8::q_sum_of(pVec2, dimension), q_dot, dimension)); } // SQ8-to-SQ8 Inner Product distance function diff --git a/src/VecSim/spaces/IP/IP_NEON_DOTPROD_UINT8.h b/src/VecSim/spaces/IP/IP_NEON_DOTPROD_UINT8.h index 3abbd9bba..71620f816 100644 --- a/src/VecSim/spaces/IP/IP_NEON_DOTPROD_UINT8.h +++ b/src/VecSim/spaces/IP/IP_NEON_DOTPROD_UINT8.h @@ -27,7 +27,7 @@ InnerProductStep(uint8_t *&pVect1, uint8_t *&pVect2, uint32x4_t &sum) { } template // 0..63 -float UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dimension) { +uint64_t UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dimension) { uint8_t *pVect1 = (uint8_t *)pVect1v; uint8_t *pVect2 = (uint8_t *)pVect2v; @@ -97,20 +97,20 @@ float UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dim uint32x4_t total_sum = vaddq_u32(sum0, sum1); - int32_t result = vaddvq_u32(total_sum); - - return static_cast(result); + // Widening horizontal sum: the lane total reaches 255*255*dim, so reducing into 32 bits + // would wrap from dimension 33,027. Callers rely on this being the exact integer dot. + return vaddlvq_u32(total_sum); } template // 0..63 float UINT8_InnerProductSIMD16_NEON_DOTPROD(const void *pVect1v, const void *pVect2v, size_t dimension) { - return 1.0f - UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + return 1.0f - static_cast(UINT8_InnerProductImp(pVect1v, pVect2v, dimension)); } template // 0..63 float UINT8_CosineSIMD_NEON_DOTPROD(const void *pVect1v, const void *pVect2v, size_t dimension) { - float ip = UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + float ip = static_cast(UINT8_InnerProductImp(pVect1v, pVect2v, dimension)); const float norm_v1 = load_unaligned(static_cast(pVect1v) + dimension); const float norm_v2 = load_unaligned(static_cast(pVect2v) + dimension); return 1.0f - ip / (norm_v1 * norm_v2); diff --git a/src/VecSim/spaces/IP/IP_NEON_SQ8_SQ8.h b/src/VecSim/spaces/IP/IP_NEON_SQ8_SQ8.h index 3e931ee0e..cbdac623b 100644 --- a/src/VecSim/spaces/IP/IP_NEON_SQ8_SQ8.h +++ b/src/VecSim/spaces/IP/IP_NEON_SQ8_SQ8.h @@ -39,30 +39,16 @@ using sq8 = vecsim_types::sq8; template // 0..63 float SQ8_SQ8_InnerProductSIMD64_NEON_IMP(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Compute raw dot product using efficient UINT8 implementation - // UINT8_InnerProductImp processes 16 elements at a time using native uint8 instructions - float dot_product = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); + // Exact integer dot product of the quantized bytes. + const uint64_t q_dot = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); - // Get dequantization parameters and precomputed values from the end of pVec1 - // Layout: [data (dim)] [min (float)] [delta (float)] [sum (float)] + // Every term describes the reconstructions, derived from the stored integer sums. See sq8.h. const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - - const auto *params1 = pVec1 + dimension; - const float min1 = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); - const float delta1 = load_unaligned(params1 + sq8::DELTA * sizeof(float)); - const float sum1 = load_unaligned(params1 + sq8::SUM * sizeof(float)); - - // Get dequantization parameters and precomputed values from the end of pVec2 - const auto *params2 = pVec2 + dimension; - const float min2 = load_unaligned(params2 + sq8::MIN_VAL * sizeof(float)); - const float delta2 = load_unaligned(params2 + sq8::DELTA * sizeof(float)); - const float sum2 = load_unaligned(params2 + sq8::SUM * sizeof(float)); - - // Apply algebraic formula using precomputed sums: - // IP = min1*sum2 + min2*sum1 + δ1*δ2 * Σ(q1*q2) - dim*min1*min2 - return min1 * sum2 + min2 * sum1 + delta1 * delta2 * dot_product - - static_cast(dimension) * min1 * min2; + return static_cast(sq8::reconstructed_ip( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::min_val_of(pVec2, dimension), + sq8::delta_of(pVec2, dimension), sq8::q_sum_of(pVec2, dimension), q_dot, dimension)); } // SQ8-to-SQ8 Inner Product distance function diff --git a/src/VecSim/spaces/IP/IP_NEON_UINT8.h b/src/VecSim/spaces/IP/IP_NEON_UINT8.h index 2d2b3f555..cb302a04e 100644 --- a/src/VecSim/spaces/IP/IP_NEON_UINT8.h +++ b/src/VecSim/spaces/IP/IP_NEON_UINT8.h @@ -35,7 +35,7 @@ InnerProductStep(uint8_t *&pVect1, uint8_t *&pVect2, uint32x4_t &sum) { } template // 0..63 -float UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dimension) { +uint64_t UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dimension) { uint8_t *pVect1 = (uint8_t *)pVect1v; uint8_t *pVect2 = (uint8_t *)pVect2v; @@ -105,20 +105,19 @@ float UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dim uint32x4_t total_sum = vaddq_u32(sum0, sum1); - // Horizontal sum of the 4 elements in the combined sum register - int32_t result = vaddvq_u32(total_sum); - - return static_cast(result); + // Widening horizontal sum: the lane total reaches 255*255*dim, so reducing into 32 bits + // would wrap from dimension 33,027. Callers rely on this being the exact integer dot. + return vaddlvq_u32(total_sum); } template // 0..15 float UINT8_InnerProductSIMD16_NEON(const void *pVect1v, const void *pVect2v, size_t dimension) { - return 1.0f - UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + return 1.0f - static_cast(UINT8_InnerProductImp(pVect1v, pVect2v, dimension)); } template // 0..63 float UINT8_CosineSIMD_NEON(const void *pVect1v, const void *pVect2v, size_t dimension) { - float ip = UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + float ip = static_cast(UINT8_InnerProductImp(pVect1v, pVect2v, dimension)); const float norm_v1 = load_unaligned(static_cast(pVect1v) + dimension); const float norm_v2 = load_unaligned(static_cast(pVect2v) + dimension); return 1.0f - ip / (norm_v1 * norm_v2); diff --git a/src/VecSim/spaces/IP/IP_SVE_SQ8_SQ8.h b/src/VecSim/spaces/IP/IP_SVE_SQ8_SQ8.h index 93ddb76cb..a61dc1c08 100644 --- a/src/VecSim/spaces/IP/IP_SVE_SQ8_SQ8.h +++ b/src/VecSim/spaces/IP/IP_SVE_SQ8_SQ8.h @@ -38,30 +38,17 @@ using sq8 = vecsim_types::sq8; // Uses UINT8_InnerProductImp for efficient dot product computation with SVE template float SQ8_SQ8_InnerProductSIMD_SVE_IMP(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Compute raw dot product using efficient UINT8 SVE implementation - // UINT8_InnerProductImp uses svdot_u32 for native uint8 dot product - float dot_product = + // Exact integer dot product of the quantized bytes. + const uint64_t q_dot = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); - // Get dequantization parameters and precomputed values from the end of vectors - // Layout: [data (dim)] [min (float)] [delta (float)] [sum (float)] + // Every term describes the reconstructions, derived from the stored integer sums. See sq8.h. const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - - const auto *params1 = pVec1 + dimension; - const float min1 = load_unaligned(params1 + sq8::MIN_VAL * sizeof(float)); - const float delta1 = load_unaligned(params1 + sq8::DELTA * sizeof(float)); - const float sum1 = load_unaligned(params1 + sq8::SUM * sizeof(float)); - - const auto *params2 = pVec2 + dimension; - const float min2 = load_unaligned(params2 + sq8::MIN_VAL * sizeof(float)); - const float delta2 = load_unaligned(params2 + sq8::DELTA * sizeof(float)); - const float sum2 = load_unaligned(params2 + sq8::SUM * sizeof(float)); - - // Apply algebraic formula with float conversion only at the end: - // IP = min1*sum2 + min2*sum1 + δ1*δ2 * Σ(q1*q2) - dim*min1*min2 - return min1 * sum2 + min2 * sum1 + delta1 * delta2 * dot_product - - static_cast(dimension) * min1 * min2; + return static_cast(sq8::reconstructed_ip( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::min_val_of(pVec2, dimension), + sq8::delta_of(pVec2, dimension), sq8::q_sum_of(pVec2, dimension), q_dot, dimension)); } // SQ8-to-SQ8 Inner Product distance function diff --git a/src/VecSim/spaces/IP/IP_SVE_UINT8.h b/src/VecSim/spaces/IP/IP_SVE_UINT8.h index f6c6af3b8..9c5c0cb73 100644 --- a/src/VecSim/spaces/IP/IP_SVE_UINT8.h +++ b/src/VecSim/spaces/IP/IP_SVE_UINT8.h @@ -24,7 +24,7 @@ inline void InnerProductStep(const uint8_t *&pVect1, const uint8_t *&pVect2, siz } template -float UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dimension) { +uint64_t UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dimension) { const uint8_t *pVect1 = reinterpret_cast(pVect1v); const uint8_t *pVect2 = reinterpret_cast(pVect2v); @@ -83,20 +83,22 @@ float UINT8_InnerProductImp(const void *pVect1v, const void *pVect2v, size_t dim sum2 = svadd_u32_x(svptrue_b32(), sum2, sum3); // Perform vector addition in parallel and Horizontal sum - int32_t sum_all = svaddv_u32(svptrue_b32(), svadd_u32_x(svptrue_b32(), sum0, sum2)); - - return sum_all; + // svaddv_u32 already reduces into a 64-bit scalar; the previous int32_t truncated it, which + // wrapped from dimension 33,027 since the total reaches 255*255*dim. Callers rely on this + // being the exact integer dot product. + return svaddv_u32(svptrue_b32(), svadd_u32_x(svptrue_b32(), sum0, sum2)); } template float UINT8_InnerProductSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t dimension) { - return 1.0f - - UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + return 1.0f - static_cast(UINT8_InnerProductImp( + pVect1v, pVect2v, dimension)); } template float UINT8_CosineSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t dimension) { - float ip = UINT8_InnerProductImp(pVect1v, pVect2v, dimension); + float ip = static_cast( + UINT8_InnerProductImp(pVect1v, pVect2v, dimension)); const float norm_v1 = load_unaligned(static_cast(pVect1v) + dimension); const float norm_v2 = load_unaligned(static_cast(pVect2v) + dimension); return 1.0f - ip / (norm_v1 * norm_v2); diff --git a/src/VecSim/spaces/L2/L2.cpp b/src/VecSim/spaces/L2/L2.cpp index 015f2200d..3d3f682d5 100644 --- a/src/VecSim/spaces/L2/L2.cpp +++ b/src/VecSim/spaces/L2/L2.cpp @@ -31,17 +31,21 @@ 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. + // ||x||² describes the reconstruction, derived exactly from the stored integer sums. const auto *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares( + sq8::min_val_of(pVect1, dimension), sq8::delta_of(pVect1, dimension), + sq8::q_sum_of(pVect1, dimension), sq8::q_sum_squares_of(pVect1, dimension), dimension); // 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; + // L2² = ||x||² + ||y||² - 2*IP(x, y). Note ip is accumulated in FP32 by the shared kernel, so + // this remains vulnerable to cancellation when ||x||² greatly exceeds the true distance; the + // asymmetric path needs the difference taken inside the loop to fix that. + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } /* @@ -62,15 +66,18 @@ float SQ8_FP16_L2Sqr(const void *pVect1v, const void *pVect2v, size_t dimension) // byte offsets that are not guaranteed 4-byte aligned for odd `dimension`, so use // load_unaligned to avoid alignment UB. const auto *pVect1 = static_cast(pVect1v); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares( + sq8::min_val_of(pVect1, dimension), sq8::delta_of(pVect1, dimension), + sq8::q_sum_of(pVect1, dimension), sq8::q_sum_squares_of(pVect1, dimension), dimension); const auto *pVect2 = static_cast(pVect2v); const auto *query_meta_bytes = reinterpret_cast(pVect2 + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + // L2² = ||x||² + ||y||² - 2*IP(x, y). Still cancellation-prone for the same reason as the FP32 + // variant: ip is accumulated in FP32 by the shared kernel. + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } float FP32_L2Sqr(const void *pVect1v, const void *pVect2v, size_t dimension) { @@ -186,16 +193,18 @@ float SQ8_SQ8_L2Sqr(const void *pVect1v, const void *pVect2v, size_t dimension) const auto *pVect1 = static_cast(pVect1v); const auto *pVect2 = static_cast(pVect2v); - // Get precomputed sum of squares from both vectors - // Layout: [uint8_t values (dim)] [min_val] [delta] [sum] [sum_of_squares] - const float sum_sq_1 = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); - const float sum_sq_2 = - load_unaligned(pVect2 + dimension + sq8::SUM_SQUARES * sizeof(float)); - - // Use the common inner product implementation - const float ip = SQ8_SQ8_InnerProduct_Impl(pVect1v, pVect2v, dimension); + // Integer dot product of the quantized bytes, exact. + uint64_t q_dot = 0; + for (size_t i = 0; i < dimension; i++) { + q_dot += static_cast(pVect1[i]) * static_cast(pVect2[i]); + } - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return sum_sq_1 + sum_sq_2 - 2.0f * ip; + // Expanded difference rather than ||x||² + ||y||² - 2*IP: see sq8::reconstructed_l2_sqr for + // why the two are not interchangeable in floating point. + return static_cast(sq8::reconstructed_l2_sqr( + sq8::min_val_of(pVect1, dimension), sq8::delta_of(pVect1, dimension), + sq8::q_sum_of(pVect1, dimension), sq8::q_sum_squares_of(pVect1, dimension), + sq8::min_val_of(pVect2, dimension), sq8::delta_of(pVect2, dimension), + sq8::q_sum_of(pVect2, dimension), sq8::q_sum_squares_of(pVect2, dimension), q_dot, + dimension)); } diff --git a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP16.h b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP16.h index c855b62ca..b03d8d8d4 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP16.h +++ b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP16.h @@ -21,12 +21,13 @@ float SQ8_FP16_L2SqrSIMD16_AVX2_FMA(const void *pVect1v, const void *pVect2v, si const uint8_t *pVect1 = static_cast(pVect1v); const uint8_t *params_bytes = pVect1 + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const float16 *pVect2 = static_cast(pVect2v); const auto *query_meta_bytes = reinterpret_cast(pVect2 + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } 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..d98b2c7a3 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h @@ -35,12 +35,12 @@ float SQ8_FP32_L2SqrSIMD16_AVX2_FMA(const void *pVect1v, const void *pVect2v, si // 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)); + const double x_sum_sq = sq8::reconstructed_sum_squares_of(pVect1, dimension); // 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]; // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP16.h b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP16.h index 7c2cbfcd8..7d5d35b8d 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP16.h +++ b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP16.h @@ -21,12 +21,13 @@ float SQ8_FP16_L2SqrSIMD16_AVX2(const void *pVect1v, const void *pVect2v, size_t const uint8_t *pVect1 = static_cast(pVect1v); const uint8_t *params_bytes = pVect1 + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const float16 *pVect2 = static_cast(pVect2v); const auto *query_meta_bytes = reinterpret_cast(pVect2 + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h index d08474f71..491ed2e98 100644 --- a/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h @@ -35,12 +35,12 @@ float SQ8_FP32_L2SqrSIMD16_AVX2(const void *pVect1v, const void *pVect2v, size_t // 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)); + const double x_sum_sq = sq8::reconstructed_sum_squares_of(pVect1, dimension); // 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]; // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } 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..61cd7e45b 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 @@ -35,12 +35,12 @@ float SQ8_FP32_L2SqrSIMD16_AVX512F_BW_VL_VNNI(const void *pVect1v, const void *p // 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)); + const double x_sum_sq = sq8::reconstructed_sum_squares_of(pVect1, dimension); // 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]; // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_SQ8.h b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_SQ8.h index 6a862b2ed..7a94d59d5 100644 --- a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_SQ8.h +++ b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_SQ8.h @@ -8,6 +8,7 @@ */ #pragma once #include "VecSim/spaces/space_includes.h" +#include "VecSim/types/sq8.h" #include "VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_SQ8_SQ8.h" /** @@ -26,19 +27,19 @@ template // 0..63 float SQ8_SQ8_L2SqrSIMD64_AVX512F_BW_VL_VNNI(const void *pVec1v, const void *pVec2v, size_t dimension) { - - // Use the common inner product implementation (returns raw IP, not distance) - const float ip = SQ8_SQ8_InnerProductImp(pVec1v, pVec2v, dimension); + // Exact integer dot product of the quantized bytes. Taken from the shared UINT8 helper rather + // than from the SQ8 inner-product wrapper, whose FP32 return would reintroduce the rounding + // this formulation exists to avoid. + const uint64_t q_dot = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - // Get precomputed sum of squares from both vectors - // Layout: [uint8_t values (dim)] [min_val] [delta] [sum] [sum_of_squares] - const float sum_sq_1 = - load_unaligned(pVec1 + dimension + sq8::SUM_SQUARES * sizeof(float)); - const float sum_sq_2 = - load_unaligned(pVec2 + dimension + sq8::SUM_SQUARES * sizeof(float)); - - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return sum_sq_1 + sum_sq_2 - 2.0f * ip; + // Expanded difference rather than ||x||^2 + ||y||^2 - 2*IP: algebraically identical, but the + // expansion keeps any common offset in (min1 - min2) instead of letting it cancel the answer. + return static_cast(sq8::reconstructed_l2_sqr( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::q_sum_squares_of(pVec1, dimension), + sq8::min_val_of(pVec2, dimension), sq8::delta_of(pVec2, dimension), + sq8::q_sum_of(pVec2, dimension), sq8::q_sum_squares_of(pVec2, dimension), q_dot, + dimension)); } diff --git a/src/VecSim/spaces/L2/L2_AVX512F_SQ8_FP16.h b/src/VecSim/spaces/L2/L2_AVX512F_SQ8_FP16.h index 9d7b1569f..0f4ed27f9 100644 --- a/src/VecSim/spaces/L2/L2_AVX512F_SQ8_FP16.h +++ b/src/VecSim/spaces/L2/L2_AVX512F_SQ8_FP16.h @@ -21,12 +21,13 @@ float SQ8_FP16_L2SqrSIMD16_AVX512F(const void *pVect1v, const void *pVect2v, siz const uint8_t *pVect1 = static_cast(pVect1v); const uint8_t *params_bytes = pVect1 + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const float16 *pVect2 = static_cast(pVect2v); const auto *query_meta_bytes = reinterpret_cast(pVect2 + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_NEON_DOTPROD_SQ8_SQ8.h b/src/VecSim/spaces/L2/L2_NEON_DOTPROD_SQ8_SQ8.h index 6c55b3f03..b9d8a19dd 100644 --- a/src/VecSim/spaces/L2/L2_NEON_DOTPROD_SQ8_SQ8.h +++ b/src/VecSim/spaces/L2/L2_NEON_DOTPROD_SQ8_SQ8.h @@ -28,20 +28,19 @@ using sq8 = vecsim_types::sq8; // L2 squared distance using the common inner product implementation template // 0..63 float SQ8_SQ8_L2SqrSIMD64_NEON_DOTPROD(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Use the common inner product implementation (returns raw IP, not distance) - const float ip = - SQ8_SQ8_InnerProductSIMD64_NEON_DOTPROD_IMP(pVec1v, pVec2v, dimension); + // Exact integer dot product of the quantized bytes. Taken from the shared UINT8 helper rather + // than from the SQ8 inner-product wrapper, whose FP32 return would reintroduce the rounding + // this formulation exists to avoid. + const uint64_t q_dot = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - - // Get precomputed sum of squares from both vectors - // Layout: [uint8_t values (dim)] [min_val] [delta] [sum] [sum_of_squares] - const float sum_sq_1 = - load_unaligned(pVec1 + dimension + sq8::SUM_SQUARES * sizeof(float)); - const float sum_sq_2 = - load_unaligned(pVec2 + dimension + sq8::SUM_SQUARES * sizeof(float)); - - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return sum_sq_1 + sum_sq_2 - 2.0f * ip; + // Expanded difference rather than ||x||^2 + ||y||^2 - 2*IP: algebraically identical, but the + // expansion keeps any common offset in (min1 - min2) instead of letting it cancel the answer. + return static_cast(sq8::reconstructed_l2_sqr( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::q_sum_squares_of(pVec1, dimension), + sq8::min_val_of(pVec2, dimension), sq8::delta_of(pVec2, dimension), + sq8::q_sum_of(pVec2, dimension), sq8::q_sum_squares_of(pVec2, dimension), q_dot, + dimension)); } diff --git a/src/VecSim/spaces/L2/L2_NEON_SQ8_FP16.h b/src/VecSim/spaces/L2/L2_NEON_SQ8_FP16.h index 2964c1cee..a5922361c 100644 --- a/src/VecSim/spaces/L2/L2_NEON_SQ8_FP16.h +++ b/src/VecSim/spaces/L2/L2_NEON_SQ8_FP16.h @@ -24,14 +24,15 @@ float SQ8_FP16_L2SqrSIMD16_NEON_HP(const void *pVect1v, const void *pVect2v, siz const float ip = SQ8_FP16_InnerProductSIMD16_NEON_HP_IMP(pVect1v, pVect2v, dimension); const uint8_t *params_bytes = static_cast(pVect1v) + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const uint8_t *query_meta_bytes = reinterpret_cast(static_cast(pVect2v) + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } // FMLAL (FEAT_FHM) variant — same identity, FMLAL widening-FMA IP core. @@ -41,12 +42,13 @@ float SQ8_FP16_L2SqrSIMD16_NEON_FHM(const void *pVect1v, const void *pVect2v, si SQ8_FP16_InnerProductSIMD16_NEON_HP_IMP(pVect1v, pVect2v, dimension); const uint8_t *params_bytes = static_cast(pVect1v) + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const uint8_t *query_meta_bytes = reinterpret_cast(static_cast(pVect2v) + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h index f6f7a6bc0..273bf4d74 100644 --- a/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h @@ -36,12 +36,12 @@ float SQ8_FP32_L2SqrSIMD16_NEON(const void *pVect1v, const void *pVect2v, size_t // 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)); + const double x_sum_sq = sq8::reconstructed_sum_squares_of(pVect1, dimension); // 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]; // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_NEON_SQ8_SQ8.h b/src/VecSim/spaces/L2/L2_NEON_SQ8_SQ8.h index 2750997f7..bcb651db3 100644 --- a/src/VecSim/spaces/L2/L2_NEON_SQ8_SQ8.h +++ b/src/VecSim/spaces/L2/L2_NEON_SQ8_SQ8.h @@ -28,19 +28,19 @@ using sq8 = vecsim_types::sq8; // L2 squared distance using the common inner product implementation template // 0..63 float SQ8_SQ8_L2SqrSIMD64_NEON(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Use the common inner product implementation (returns raw IP, not distance) - const float ip = SQ8_SQ8_InnerProductSIMD64_NEON_IMP(pVec1v, pVec2v, dimension); + // Exact integer dot product of the quantized bytes. Taken from the shared UINT8 helper rather + // than from the SQ8 inner-product wrapper, whose FP32 return would reintroduce the rounding + // this formulation exists to avoid. + const uint64_t q_dot = UINT8_InnerProductImp(pVec1v, pVec2v, dimension); const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - - // Get precomputed sum of squares from both vectors - // Layout: [uint8_t values (dim)] [min_val] [delta] [sum] [sum_of_squares] - const float sum_sq_1 = - load_unaligned(pVec1 + dimension + sq8::SUM_SQUARES * sizeof(float)); - const float sum_sq_2 = - load_unaligned(pVec2 + dimension + sq8::SUM_SQUARES * sizeof(float)); - - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return sum_sq_1 + sum_sq_2 - 2.0f * ip; + // Expanded difference rather than ||x||^2 + ||y||^2 - 2*IP: algebraically identical, but the + // expansion keeps any common offset in (min1 - min2) instead of letting it cancel the answer. + return static_cast(sq8::reconstructed_l2_sqr( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::q_sum_squares_of(pVec1, dimension), + sq8::min_val_of(pVec2, dimension), sq8::delta_of(pVec2, dimension), + sq8::q_sum_of(pVec2, dimension), sq8::q_sum_squares_of(pVec2, dimension), q_dot, + dimension)); } diff --git a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP16.h b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP16.h index d0a0fea06..3d80a39ef 100644 --- a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP16.h +++ b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP16.h @@ -20,12 +20,13 @@ float SQ8_FP16_L2SqrSIMD16_SSE4(const void *pVect1v, const void *pVect2v, size_t const uint8_t *pVect1 = static_cast(pVect1v); const uint8_t *params_bytes = pVect1 + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const float16 *pVect2 = static_cast(pVect2v); const auto *query_meta_bytes = reinterpret_cast(pVect2 + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h index 3a3a4d12d..943ecf2b0 100644 --- a/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h @@ -35,12 +35,12 @@ float SQ8_FP32_L2SqrSIMD16_SSE4(const void *pVect1v, const void *pVect2v, size_t // 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)); + const double x_sum_sq = sq8::reconstructed_sum_squares_of(pVect1, dimension); // 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]; // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_SVE2_SQ8_FP16.h b/src/VecSim/spaces/L2/L2_SVE2_SQ8_FP16.h index d9451fe2a..3c0cdf343 100644 --- a/src/VecSim/spaces/L2/L2_SVE2_SQ8_FP16.h +++ b/src/VecSim/spaces/L2/L2_SVE2_SQ8_FP16.h @@ -22,11 +22,12 @@ float SQ8_FP16_L2SqrSIMD_SVE2(const void *pVect1v, const void *pVect2v, size_t d pVect1v, pVect2v, dimension); const uint8_t *params_bytes = static_cast(pVect1v) + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const uint8_t *query_meta_bytes = reinterpret_cast(static_cast(pVect2v) + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP16.h b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP16.h index f70ef493d..4848b5ceb 100644 --- a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP16.h +++ b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP16.h @@ -22,11 +22,12 @@ float SQ8_FP16_L2SqrSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t di pVect1v, pVect2v, dimension); const uint8_t *params_bytes = static_cast(pVect1v) + dimension; - const float x_sum_sq = load_unaligned(params_bytes + sq8::SUM_SQUARES * sizeof(float)); + const double x_sum_sq = sq8::reconstructed_sum_squares_at(params_bytes, dimension); const uint8_t *query_meta_bytes = reinterpret_cast(static_cast(pVect2v) + dimension); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h index 3f95c4ae4..cd1c979dc 100644 --- a/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h +++ b/src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h @@ -37,12 +37,12 @@ float SQ8_FP32_L2SqrSIMD_SVE(const void *pVect1v, const void *pVect2v, size_t di // 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)); + const double x_sum_sq = sq8::reconstructed_sum_squares_of(pVect1, dimension); // 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]; // L2² = ||x||² + ||y||² - 2*IP(x, y) - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } diff --git a/src/VecSim/spaces/L2/L2_SVE_SQ8_SQ8.h b/src/VecSim/spaces/L2/L2_SVE_SQ8_SQ8.h index 5fe194d80..a8a906b03 100644 --- a/src/VecSim/spaces/L2/L2_SVE_SQ8_SQ8.h +++ b/src/VecSim/spaces/L2/L2_SVE_SQ8_SQ8.h @@ -28,20 +28,20 @@ using sq8 = vecsim_types::sq8; // L2 squared distance using the common inner product implementation template float SQ8_SQ8_L2SqrSIMD_SVE(const void *pVec1v, const void *pVec2v, size_t dimension) { - // Use the common inner product implementation (returns raw IP, not distance) - const float ip = SQ8_SQ8_InnerProductSIMD_SVE_IMP( - pVec1v, pVec2v, dimension); + // Exact integer dot product of the quantized bytes. Taken from the shared UINT8 helper rather + // than from the SQ8 inner-product wrapper, whose FP32 return would reintroduce the rounding + // this formulation exists to avoid. + const uint64_t q_dot = + UINT8_InnerProductImp(pVec1v, pVec2v, dimension); const uint8_t *pVec1 = static_cast(pVec1v); const uint8_t *pVec2 = static_cast(pVec2v); - - // Get precomputed sum of squares from both vectors - // Layout: [uint8_t values (dim)] [min_val] [delta] [sum] [sum_of_squares] - const float sum_sq_1 = - load_unaligned(pVec1 + dimension + sq8::SUM_SQUARES * sizeof(float)); - const float sum_sq_2 = - load_unaligned(pVec2 + dimension + sq8::SUM_SQUARES * sizeof(float)); - - // L2² = ||x||² + ||y||² - 2*IP(x, y) - return sum_sq_1 + sum_sq_2 - 2.0f * ip; + // Expanded difference rather than ||x||^2 + ||y||^2 - 2*IP: algebraically identical, but the + // expansion keeps any common offset in (min1 - min2) instead of letting it cancel the answer. + return static_cast(sq8::reconstructed_l2_sqr( + sq8::min_val_of(pVec1, dimension), sq8::delta_of(pVec1, dimension), + sq8::q_sum_of(pVec1, dimension), sq8::q_sum_squares_of(pVec1, dimension), + sq8::min_val_of(pVec2, dimension), sq8::delta_of(pVec2, dimension), + sq8::q_sum_of(pVec2, dimension), sq8::q_sum_squares_of(pVec2, dimension), q_dot, + dimension)); } diff --git a/src/VecSim/spaces/computer/preprocessors.h b/src/VecSim/spaces/computer/preprocessors.h index 195be1418..acdc4bfe0 100644 --- a/src/VecSim/spaces/computer/preprocessors.h +++ b/src/VecSim/spaces/computer/preprocessors.h @@ -274,13 +274,15 @@ class QuantPreprocessor : public PreprocessorInterface { const MetadataType delta = (diff == 0.0f) ? MetadataType{1} : diff / MetadataType{255}; const MetadataType inv_delta = MetadataType{1} / delta; - // Compute sum (and sum of squares for L2) while quantizing. - // Accumulators are FP32 to preserve metadata precision for FP16 inputs. - // 4 independent accumulators (sum) - float s0{}, s1{}, s2{}, s3{}; + // Sum the *quantized* bytes, not the input floats. A stored blob only ever describes its + // reconstruction min + delta * a[i], so metadata derived from it has to describe that same + // vector; summing the input instead mixes two different vectors into one distance formula. + // Integer accumulation is also exact, which the L2 identity depends on (see sq8.h). + // 4 independent accumulators each, mirroring the unrolled loop below. + uint32_t s0{}, s1{}, s2{}, s3{}; - // 4 independent accumulators (sum of squares), only used for L2 - float q0{}, q1{}, q2{}, q3{}; + // Sum of squares, only used for L2. + uint32_t q0{}, q1{}, q2{}, q3{}; size_t i = 0; // round dim down to the nearest multiple of 4 @@ -296,50 +298,66 @@ class QuantPreprocessor : public PreprocessorInterface { // We know (input - min) => 0 // If min == max, all values are the same and should be quantized to 0. // reconstruction will yield the same original value for all vectors. - quantized[i] = static_cast(std::round((x0 - min_val) * inv_delta)); - quantized[i + 1] = static_cast(std::round((x1 - min_val) * inv_delta)); - quantized[i + 2] = static_cast(std::round((x2 - min_val) * inv_delta)); - quantized[i + 3] = static_cast(std::round((x3 - min_val) * inv_delta)); + const uint32_t a0 = static_cast(std::round((x0 - min_val) * inv_delta)); + const uint32_t a1 = static_cast(std::round((x1 - min_val) * inv_delta)); + const uint32_t a2 = static_cast(std::round((x2 - min_val) * inv_delta)); + const uint32_t a3 = static_cast(std::round((x3 - min_val) * inv_delta)); + quantized[i] = static_cast(a0); + quantized[i + 1] = static_cast(a1); + quantized[i + 2] = static_cast(a2); + quantized[i + 3] = static_cast(a3); // Accumulate sum for all metrics - s0 += x0; - s1 += x1; - s2 += x2; - s3 += x3; + s0 += a0; + s1 += a1; + s2 += a2; + s3 += a3; // Accumulate sum of squares only for L2 metric if constexpr (Metric == VecSimMetric_L2) { - q0 += x0 * x0; - q1 += x1 * x1; - q2 += x2 * x2; - q3 += x3 * x3; + q0 += a0 * a0; + q1 += a1 * a1; + q2 += a2 * a2; + q3 += a3 * a3; } } // Tail: 0..3 remaining elements (still the same pass, just finishing work). - // Sum/sum_squares become metadata, so they are MetadataType. - MetadataType sum = (s0 + s1) + (s2 + s3); - MetadataType sum_squares = (q0 + q1) + (q2 + q3); + uint32_t q_sum = (s0 + s1) + (s2 + s3); + uint32_t q_sum_squares = (q0 + q1) + (q2 + q3); for (; i < this->dim; ++i) { const float x = transformed_value(input, i); - quantized[i] = static_cast(std::round((x - min_val) * inv_delta)); - sum += x; + const uint32_t a = static_cast(std::round((x - min_val) * inv_delta)); + quantized[i] = static_cast(a); + q_sum += a; if constexpr (Metric == VecSimMetric_L2) { - sum_squares += x * x; + q_sum_squares += a * a; } } - // Metadata uses MetadataType. Use memcpy because the metadata offset - // (dim * sizeof(uint8_t)) is not guaranteed to be sizeof(MetadataType)-aligned. - void *meta_dst = quantized + this->dim; - MetadataType buf[5] = {min_val, delta, sum}; - size_t n = 3; + // sum(a[i]^2) is exact in uint32 only up to sq8::MAX_EXACT_DIM; past it this wraps and the + // L2 distance silently degrades. Callers that build quantized indexes must fence dim. + assert(this->dim <= sq8::MAX_EXACT_DIM && "SQ8 dimension exceeds exact metadata range"); + + // Metadata is a packed run of 4-byte fields: min and delta as FP32, the quantized sums as + // uint32, and x_mean_ip (FP32) last for IP with norm. memcpy throughout because the offset + // (dim * sizeof(uint8_t)) is not guaranteed to be 4-byte aligned. + unsigned char buf[5 * sizeof(uint32_t)]; + size_t n = 0; + const auto put = [&buf, &n](const auto &value) { + static_assert(sizeof(value) == sizeof(uint32_t)); + memcpy(buf + n, &value, sizeof(value)); + n += sizeof(value); + }; + put(min_val); + put(delta); + put(q_sum); if constexpr (Metric == VecSimMetric_L2) - buf[n++] = sum_squares; + put(q_sum_squares); if constexpr (WithNorm && Metric == VecSimMetric_IP) - buf[n++] = x_mean_ip; - memcpy(meta_dst, buf, n * sizeof(MetadataType)); + put(x_mean_ip); + memcpy(quantized + this->dim, buf, n); } // Computes and writes query metadata (FP32) in a single pass over the query values. diff --git a/src/VecSim/types/sq8.h b/src/VecSim/types/sq8.h index c1e9c40b8..1721fde26 100644 --- a/src/VecSim/types/sq8.h +++ b/src/VecSim/types/sq8.h @@ -11,6 +11,7 @@ #include #include #include "VecSim/vec_sim_common.h" +#include "VecSim/utils/alignment.h" namespace vecsim_types { @@ -18,14 +19,29 @@ namespace vecsim_types { struct sq8 { using value_type = uint8_t; - // Metadata layout indices (stored after quantized values) + // Metadata layout indices (stored after quantized values). Each slot is 4 bytes. + // + // MIN_VAL and DELTA are FP32. Q_SUM and Q_SUM_SQUARES are uint32 sums over the *quantized* + // bytes, not over the input floats. That distinction is the whole point: a stored blob only + // ever describes its reconstruction x_r[i] = min + delta * a[i], so every quantity derived + // from it has to describe x_r as well. Summing the original input instead mixes two different + // vectors into one formula and produces distances that can even be negative. + // + // Integer sums are also exact, which matters more than it looks: the L2 identity + // ||x||^2 + ||y||^2 - 2*IP subtracts large nearly-equal quantities, so any rounding in the + // inputs lands directly on the answer. Keeping the sums exact and combining them in double + // (see the helpers below) is what makes that subtraction trustworthy. enum MetadataIndex : size_t { MIN_VAL = 0, DELTA = 1, - SUM = 2, - SUM_SQUARES = 3 // Only for L2 + Q_SUM = 2, // uint32: sum(a[i]) + Q_SUM_SQUARES = 3 // uint32: sum(a[i]^2). Only for L2 }; + // sum(a[i]^2) <= 255^2 * dim, which is exactly representable in uint32 up to this dimension. + // Past it the accumulation wraps, so callers that build quantized indexes must fence dim. + static constexpr size_t MAX_EXACT_DIM = 33026; + enum QueryMetadataIndex : size_t { SUM_QUERY = 0, SUM_SQUARES_QUERY = 1 // Only for L2 @@ -57,6 +73,86 @@ struct sq8 { static constexpr size_t query_mean_ip_index() { return query_metadata_count() - 1; } + + // Readers for the packed metadata run that follows the dim quantized bytes. The run is not + // guaranteed to be 4-byte aligned, so every field goes through an unaligned load. Kernels use + // these rather than open-coding the offsets, so a layout change has one place to happen. + // Taking the metadata pointer, for kernels that have already formed it. + static float min_val_at(const value_type *meta) { + return load_unaligned(meta + MIN_VAL * sizeof(float)); + } + static float delta_at(const value_type *meta) { + return load_unaligned(meta + DELTA * sizeof(float)); + } + static uint32_t q_sum_at(const value_type *meta) { + return load_unaligned(meta + Q_SUM * sizeof(uint32_t)); + } + static uint32_t q_sum_squares_at(const value_type *meta) { + return load_unaligned(meta + Q_SUM_SQUARES * sizeof(uint32_t)); + } + + // Taking the blob pointer, for kernels that have not. + static float min_val_of(const value_type *blob, size_t dim) { return min_val_at(blob + dim); } + static float delta_of(const value_type *blob, size_t dim) { return delta_at(blob + dim); } + static uint32_t q_sum_of(const value_type *blob, size_t dim) { return q_sum_at(blob + dim); } + static uint32_t q_sum_squares_of(const value_type *blob, size_t dim) { + return q_sum_squares_at(blob + dim); + } + + // Single source of truth for turning stored metadata back into properties of the + // reconstruction. Every kernel derives its terms here so the storage layout has exactly one + // interpretation. All of these take double and return double on purpose: the integer sums are + // exact, and the point is to keep them exact through the combination. + + // sum(x_r[i]) where x_r[i] = min + delta * a[i]. + static double reconstructed_sum(double min_val, double delta, uint32_t q_sum, size_t dim) { + return static_cast(dim) * min_val + delta * static_cast(q_sum); + } + + // sum(x_r[i]^2), expanded so that no term is ever formed from a rounded input. + static double reconstructed_sum_squares(double min_val, double delta, uint32_t q_sum, + uint32_t q_sum_squares, size_t dim) { + return static_cast(dim) * min_val * min_val + + 2.0 * min_val * delta * static_cast(q_sum) + + delta * delta * static_cast(q_sum_squares); + } + + // ||x_r||^2 for one stored blob. Two entry points because some kernels have the blob pointer + // and some have already formed the metadata pointer. + static double reconstructed_sum_squares_at(const value_type *meta, size_t dim) { + return reconstructed_sum_squares(min_val_at(meta), delta_at(meta), q_sum_at(meta), + q_sum_squares_at(meta), dim); + } + static double reconstructed_sum_squares_of(const value_type *blob, size_t dim) { + return reconstructed_sum_squares_at(blob + dim, dim); + } + + // IP(x_r, y_r) for two quantized blobs, from their exact integer sums and integer dot product. + static double reconstructed_ip(double min1, double delta1, uint32_t q_sum1, double min2, + double delta2, uint32_t q_sum2, uint64_t q_dot, size_t dim) { + return static_cast(dim) * min1 * min2 + + min1 * delta2 * static_cast(q_sum2) + + min2 * delta1 * static_cast(q_sum1) + + delta1 * delta2 * static_cast(q_dot); + } + + // ||x_r - y_r||^2 for two quantized blobs. Expanding the difference before summing keeps any + // common offset in the (min1 - min2) term, where it is formed once and exactly, instead of + // letting it inflate the summed terms and cancel away the answer. That is why this is not + // written as reconstructed_sum_squares + reconstructed_sum_squares - 2 * reconstructed_ip: + // algebraically identical, numerically not. + static double reconstructed_l2_sqr(double min1, double delta1, uint32_t q_sum1, + uint32_t q_sum_squares1, double min2, double delta2, + uint32_t q_sum2, uint32_t q_sum_squares2, uint64_t q_dot, + size_t dim) { + const double c = min1 - min2; + return static_cast(dim) * c * c + + delta1 * delta1 * static_cast(q_sum_squares1) + + delta2 * delta2 * static_cast(q_sum_squares2) + + 2.0 * c * delta1 * static_cast(q_sum1) - + 2.0 * c * delta2 * static_cast(q_sum2) - + 2.0 * delta1 * delta2 * static_cast(q_dot); + } }; } // namespace vecsim_types From 2b84cc2fba5e4a7226d029b2bd75c2c0bb5ad1c4 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Mon, 10 Aug 2026 17:04:34 +0300 Subject: [PATCH 02/17] Update the SQ8 test fixtures for the exact integer metadata The fixtures hand-build SQ8 blobs, so they encoded the old semantics. * tests_utils.h: the blob builder sums the quantized bytes as integers, and the FP16 reference L2 derives ||x||^2 through sq8::reconstructed_sum_squares_of and combines in double, matching the kernel it is a reference for. * test_spaces.cpp: the unaligned-metadata case stores uint32 sums. Its values are unchanged, since min 0 and delta 1 make the reconstruction equal to the bytes, so the expected distance does not move. It still asserts the metadata addresses are not float-aligned, which the uint32 loads must also tolerate. * test_components.cpp: adds a uint32 metadata reader beside the float one, and compares the FP16 path against the FP32 baseline on the integer sums. * The degenerate min == max case now expects both sums to be 0, because every value collapses to byte 0, and additionally asserts reconstructed_sum still recovers 3.5 * dim. That expectation moved because the quantity changed meaning, not because the number was adjusted to fit. Co-Authored-By: Claude Opus 5 --- .sites.txt | 25 +++++++++++++++++++++++ tests/unit/test_components.cpp | 37 +++++++++++++++++++++------------- tests/unit/test_spaces.cpp | 14 ++++++++----- tests/utils/tests_utils.h | 35 ++++++++++++++++---------------- 4 files changed, 75 insertions(+), 36 deletions(-) create mode 100644 .sites.txt diff --git a/.sites.txt b/.sites.txt new file mode 100644 index 000000000..9bf88f993 --- /dev/null +++ b/.sites.txt @@ -0,0 +1,25 @@ +src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_SQ8_SQ8.h +src/VecSim/spaces/IP/IP_NEON_DOTPROD_SQ8_SQ8.h +src/VecSim/spaces/IP/IP_NEON_SQ8_SQ8.h +src/VecSim/spaces/IP/IP_SVE_SQ8_SQ8.h +src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP16.h +src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h +src/VecSim/spaces/L2/L2_AVX2_SQ8_FP16.h +src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h +src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h +src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_SQ8.h +src/VecSim/spaces/L2/L2_AVX512F_SQ8_FP16.h +src/VecSim/spaces/L2/L2.cpp +src/VecSim/spaces/L2/L2_NEON_DOTPROD_SQ8_SQ8.h +src/VecSim/spaces/L2/L2_NEON_SQ8_FP16.h +src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h +src/VecSim/spaces/L2/L2_NEON_SQ8_SQ8.h +src/VecSim/spaces/L2/L2_SSE4_SQ8_FP16.h +src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h +src/VecSim/spaces/L2/L2_SVE2_SQ8_FP16.h +src/VecSim/spaces/L2/L2_SVE_SQ8_FP16.h +src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h +src/VecSim/spaces/L2/L2_SVE_SQ8_SQ8.h +tests/unit/test_components.cpp +tests/unit/test_spaces.cpp +tests/utils/tests_utils.h diff --git a/tests/unit/test_components.cpp b/tests/unit/test_components.cpp index efb39615d..5a5ddbabf 100644 --- a/tests/unit/test_components.cpp +++ b/tests/unit/test_components.cpp @@ -1233,12 +1233,12 @@ TEST(PreprocessorsTest, QuantizationTestAllEntriesEqual) { ASSERT_FLOAT_EQ(min_val, 3.5f); ASSERT_FLOAT_EQ(delta, 1.0f); - // Verify sum and sum_squares for L2 metric - float expected_sum = 3.5f * dim; - float expected_sum_squares = 3.5f * 3.5f * dim; - ASSERT_FLOAT_EQ(load_unaligned(metadata + sq8::SUM * sizeof(float)), expected_sum); - ASSERT_FLOAT_EQ(load_unaligned(metadata + sq8::SUM_SQUARES * sizeof(float)), - expected_sum_squares); + // The sums are over the quantized bytes, and every value collapsed to 0 here, so both are 0. + // The vector itself is still recovered exactly: reconstructed_sum gives dim * min = 3.5 * dim, + // which is what any consumer of the metadata actually asks for. + ASSERT_EQ(load_unaligned(metadata + sq8::Q_SUM * sizeof(uint32_t)), 0u); + ASSERT_EQ(load_unaligned(metadata + sq8::Q_SUM_SQUARES * sizeof(uint32_t)), 0u); + ASSERT_DOUBLE_EQ(sq8::reconstructed_sum(min_val, delta, 0u, dim), 3.5 * dim); // Reconstruct and verify: min + quantized * delta = 3.5 + 0 * 1 = 3.5 for (size_t i = 0; i < dim; ++i) { @@ -1444,6 +1444,13 @@ class QuantPreprocessorFP16MetricTest : public testing::TestWithParam(base) + byte_offset, sizeof(v)); + return v; + } + template void runQuantizationTest() { const size_t expected_storage_size = getExpectedStorageSize(); @@ -1459,9 +1466,10 @@ class QuantPreprocessorFP16MetricTest : public testing::TestWithParam(allocator, dim); @@ -1496,12 +1504,13 @@ class QuantPreprocessorFP16MetricTest : public testing::TestWithParam(i + 2)); @@ -518,7 +522,7 @@ TEST_F(SpacesTest, SQ8_FP16_l2sqr_odd_dim_unaligned_metadata_test) { store_float(query_meta + sq8::SUM_QUERY * sizeof(float), query_sum); store_float(query_meta + sq8::SUM_SQUARES_QUERY * sizeof(float), query_sum_squares); - const auto *storage_sum_squares_addr = storage_meta + sq8::SUM_SQUARES * sizeof(float); + const auto *storage_sum_squares_addr = storage_meta + sq8::Q_SUM_SQUARES * sizeof(uint32_t); const auto *query_sum_squares_addr = query_meta + sq8::SUM_SQUARES_QUERY * sizeof(float); ASSERT_NE(reinterpret_cast(storage_sum_squares_addr) % alignof(float), 0u); ASSERT_NE(reinterpret_cast(query_sum_squares_addr) % alignof(float), 0u); diff --git a/tests/utils/tests_utils.h b/tests/utils/tests_utils.h index 67d93177c..d7fa5ff62 100644 --- a/tests/utils/tests_utils.h +++ b/tests/utils/tests_utils.h @@ -181,31 +181,32 @@ static void quantize_float_vec_to_sq8_with_metadata(const float *v, size_t dim, max_val = std::max(max_val, v[i]); } - float sum = 0.0f; - float square_sum = 0.0f; - for (size_t i = 0; i < dim; i++) { - sum += v[i]; - square_sum += v[i] * v[i]; - } - // Calculate delta float delta = (max_val - min_val) / 255.0f; if (delta == 0) delta = 1.0f; // Avoid division by zero - // Quantize each value + // Quantize each value, summing the quantized bytes as we go. The sums describe the + // reconstruction rather than the input, and they are integers, matching the layout the kernels + // read (see sq8.h). + uint32_t q_sum = 0; + uint32_t q_sum_squares = 0; for (size_t i = 0; i < dim; i++) { float normalized = (v[i] - min_val) / delta; normalized = std::max(0.0f, std::min(255.0f, normalized)); - qv[i] = static_cast(std::round(normalized)); + const uint32_t q = static_cast(std::round(normalized)); + qv[i] = static_cast(q); + q_sum += q; + q_sum_squares += q * q; } - // Store parameters: [min, delta, sum, square_sum] + // Store parameters: [min, delta, q_sum, q_sum_squares] auto *params = qv + dim; std::memcpy(params + sq8::MIN_VAL * sizeof(float), &min_val, sizeof(float)); std::memcpy(params + sq8::DELTA * sizeof(float), &delta, sizeof(float)); - std::memcpy(params + sq8::SUM * sizeof(float), &sum, sizeof(float)); - std::memcpy(params + sq8::SUM_SQUARES * sizeof(float), &square_sum, sizeof(float)); + std::memcpy(params + sq8::Q_SUM * sizeof(uint32_t), &q_sum, sizeof(q_sum)); + std::memcpy(params + sq8::Q_SUM_SQUARES * sizeof(uint32_t), &q_sum_squares, + sizeof(q_sum_squares)); } // Preprocess fp32 query for SQ8 IP/Cosine/L2 space. @@ -315,10 +316,9 @@ static float SQ8_FP16_NotOptimized_L2Sqr(const void *pVect1v, const void *pVect2 // Storage and query metadata sit at byte offsets that are not guaranteed 4-byte aligned // for odd `dimension`; use load_unaligned to avoid alignment UB. - const float min_val = load_unaligned(pVect1 + dimension + sq8::MIN_VAL * sizeof(float)); - const float delta = load_unaligned(pVect1 + dimension + sq8::DELTA * sizeof(float)); - const float x_sum_sq = - load_unaligned(pVect1 + dimension + sq8::SUM_SQUARES * sizeof(float)); + const float min_val = sq8::min_val_of(pVect1, dimension); + const float delta = sq8::delta_of(pVect1, dimension); + const double x_sum_sq = sq8::reconstructed_sum_squares_of(pVect1, dimension); const auto *query_meta_bytes = pVect2 + dimension * sizeof(vecsim_types::float16); const float y_sum_sq = load_unaligned(query_meta_bytes + sq8::SUM_SQUARES_QUERY * sizeof(float)); @@ -330,7 +330,8 @@ static float SQ8_FP16_NotOptimized_L2Sqr(const void *pVect1v, const void *pVect2 const float dequantized = pVect1[i] * delta + min_val; ip += dequantized * vecsim_types::FP16_to_FP32(vecsim_types::float16{raw}); } - return x_sum_sq + y_sum_sq - 2.0f * ip; + return static_cast(x_sum_sq + static_cast(y_sum_sq) - + 2.0 * static_cast(ip)); } // Preprocess FP16 query for SQ8 IP/Cosine/L2 space. From 50a395e20e31be850ce6a85074b8219fe741dc72 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Mon, 10 Aug 2026 17:06:22 +0300 Subject: [PATCH 03/17] Give the query-metadata assertions their own baseline The FP16 preprocessor test compared query metadata against the storage sums. That only worked because both were taken over the input values; the storage sums now describe the reconstruction, so the query side needs its own expectation computed from the widened query values. Queries are never quantized, so their metadata stays FP32. Co-Authored-By: Claude Opus 5 --- tests/unit/test_components.cpp | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/tests/unit/test_components.cpp b/tests/unit/test_components.cpp index 5a5ddbabf..482623bba 100644 --- a/tests/unit/test_components.cpp +++ b/tests/unit/test_components.cpp @@ -1471,6 +1471,17 @@ class QuantPreprocessorFP16MetricTest : public testing::TestWithParam(allocator, dim); @@ -1520,11 +1531,11 @@ class QuantPreprocessorFP16MetricTest : public testing::TestWithParamfree_allocation(storage_blob); @@ -1542,11 +1553,11 @@ class QuantPreprocessorFP16MetricTest : public testing::TestWithParam(static_cast(blob), original_blob, dim)); ASSERT_FLOAT_EQ(load_meta(blob, query_meta_offset + sq8::SUM_QUERY * sizeof(float)), - baseline_sum); + baseline_y_sum); if constexpr (Metric == VecSimMetric_L2) { ASSERT_FLOAT_EQ( load_meta(blob, query_meta_offset + sq8::SUM_SQUARES_QUERY * sizeof(float)), - baseline_sum_sq); + baseline_y_sum_sq); } allocator->free_allocation(blob); } From e405dd76c8e87b2931aa08e792f5e6b48b46e6f3 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Mon, 10 Aug 2026 17:10:06 +0300 Subject: [PATCH 04/17] Sum the quantized bytes in the reference quantizer too ComputeSQ8Quantization builds the expected blob that the preprocessor tests byte-compare against, and it still wrote FP32 sums over the input. It did not reference the layout enum, so renaming the slots did not catch it; the byte comparison did. Co-Authored-By: Claude Opus 5 --- tests/unit/unit_test_utils.h | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/tests/unit/unit_test_utils.h b/tests/unit/unit_test_utils.h index 0039dea4b..fa8312595 100644 --- a/tests/unit/unit_test_utils.h +++ b/tests/unit/unit_test_utils.h @@ -250,20 +250,24 @@ inline void ComputeSQ8Quantization(const float *original_blob, size_t dim, uint8 float diff = max_val - min_val; float delta = (diff == 0.0f) ? 1.0f : diff / 255.0f; - // Calculate quantized values, sum and sum_squares - float sum = 0.0f; - float sum_squares = 0.0f; + // Calculate quantized values, and sum them as integers: the stored sums describe the + // reconstruction rather than the input, and are exact (see sq8.h). + uint32_t q_sum = 0; + uint32_t q_sum_squares = 0; for (size_t i = 0; i < dim; i++) { float normalized = (original_blob[i] - min_val) / delta; - output[i] = static_cast(std::round(normalized)); - sum += original_blob[i]; - sum_squares += original_blob[i] * original_blob[i]; + const uint32_t q = static_cast(std::round(normalized)); + output[i] = static_cast(q); + q_sum += q; + q_sum_squares += q * q; } - // Store metadata: min_val, delta, sum, sum_squares. Use memcpy because the metadata region - // (output + dim) is not guaranteed to be 4-byte aligned for arbitrary dim values. - const float metadata[4] = {min_val, delta, sum, sum_squares}; - std::memcpy(output + dim, metadata, sizeof(metadata)); + // Store metadata: min_val and delta as FP32, then the two integer sums. Use memcpy because the + // metadata region (output + dim) is not guaranteed to be 4-byte aligned for arbitrary dim. + const float head[2] = {min_val, delta}; + const uint32_t sums[2] = {q_sum, q_sum_squares}; + std::memcpy(output + dim, head, sizeof(head)); + std::memcpy(output + dim + sizeof(head), sums, sizeof(sums)); } // TODO: Move all test_utils to this namespace From 5ad4f03b0ac93a458a53beba7711f6321047cbd5 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Mon, 10 Aug 2026 17:27:38 +0300 Subject: [PATCH 05/17] Fix the NaN-to-uint8 UB, and pin all four defects with regressions Folds in the quantization UB, since it lives in the same function and is the same question of whether the preprocessor's output can be trusted. Finite FP32 input whose range is not representable, [-FLT_MAX, +FLT_MAX], made max - min overflow to inf, then delta inf, inv_delta 0, and inf * 0 = NaN, whose conversion to an integer is undefined behaviour. The range and the per-element scaling are now computed in double, which cannot overflow for any pair of finite floats, so the class disappears instead of being special-cased. The result is also clamped to [0, 255] before conversion, so a value the arithmetic should never produce would be a wrong byte rather than UB. Note the reformulation x * inv_delta - min * inv_delta was rejected: it avoids the overflow but reintroduces cancellation for ordinary vectors with a large offset, quantizing them coarsely. Regressions, all written against an FP64 reference over the reconstructed vectors rather than against numbers recorded from a run: * SQ8_SQ8_L2_survives_large_common_offset: stored [100000, 100008] against [100000, 100000], which returned 0 where the truth is 64. * SQ8_FP32_L2_is_non_negative_for_off_grid_query: stored [0, 0.25, 1] against [0, 0.2501, 1], which returned -0.000490427 and so matched a radius-0 query. * SQ8_SQ8_L2_matches_fp64_reference_at_high_dim: dim 16384 near-identical vectors, the regime where combining even the exact sums in FP32 was 16% off. * QuantizationHandlesNonRepresentableRange: the UB case above, under UBSan in CI. Also brought the QuantPreprocessor class documentation back in line with the code: it still described sums over the original values and presented the L2 identity as the symmetric formulation. Co-Authored-By: Claude Opus 5 --- src/VecSim/spaces/computer/preprocessors.h | 74 +++++++++++------ tests/unit/test_components.cpp | 41 ++++++++++ tests/unit/test_spaces.cpp | 94 ++++++++++++++++++++++ 3 files changed, 186 insertions(+), 23 deletions(-) diff --git a/src/VecSim/spaces/computer/preprocessors.h b/src/VecSim/spaces/computer/preprocessors.h index acdc4bfe0..a670d52f2 100644 --- a/src/VecSim/spaces/computer/preprocessors.h +++ b/src/VecSim/spaces/computer/preprocessors.h @@ -151,10 +151,15 @@ class CosinePreprocessor : public PreprocessorInterface { * and then scaling the values to fit in the range of [0, 255]. * * Storage layout: - * | quantized_values[dim] | min_val | delta | x_sum | (x_sum_squares for L2 only) | + * | quantized_values[dim] | min_val | delta | q_sum | (q_sum_squares for L2 only) | * where: - * x_sum = Σx_i: sum of the original values, - * x_sum_squares = Σx_i²: sum of squares of the original values. + * q_sum = Σa_i: sum of the quantized bytes (uint32), + * q_sum_squares = Σa_i²: sum of their squares (uint32). + * + * The sums are over the quantized bytes, not the original values, because a stored blob only ever + * describes its reconstruction min + delta * a_i. Anything derived from the blob must describe that + * same vector; sums over the input belong to a different one. They are integers so that they stay + * exact, which the L2 formulation depends on. See sq8.h for the derivations that consume them. * * Storage metadata is always FP32 (independent of DataType) to match the asymmetric distance * kernels. The quantized blob size is: @@ -195,11 +200,14 @@ class CosinePreprocessor : public PreprocessorInterface { * * For L2: * ||x - y||² = Σx_i² - 2*Σ(x_i * y_i) + Σy_i² - * = x_sum_squares - 2 * IP(x, y) + y_sum_squares * where: - * - x_sum_squares = Σx_i² is precomputed and stored in the storage blob + * - Σx_i² describes the reconstruction and is derived from q_sum and q_sum_squares by + * sq8::reconstructed_sum_squares, in double * - IP(x, y) is computed using the formula above - * - y_sum_squares = Σy_i² is precomputed and stored in the query blob + * - y_sum_squares = Σy_i² is precomputed and stored in the query blob, in FP32, because the + * query is never quantized + * Note this identity is only safe while Σx_i² is close to the answer. It is not when the two + * vectors share a large offset, which is why the symmetric case below does not use it. * For normalized L2, x and y in this formula are the centered values x' and y'; their distance * is identical to the distance between the original vectors. * @@ -213,16 +221,20 @@ class CosinePreprocessor : public PreprocessorInterface { * = dim * min_x * min_y * + min_x * (sum_y - dim * min_y) + min_y * (sum_x - dim * min_x) * + delta_x * delta_y * Σ(qx_i * qy_i) - * = min_x * sum_y + min_y * sum_x - dim * min_x * min_y - * + delta_x * delta_y * Σ(qx_i * qy_i) - * where: - * - sum_x, sum_y are precomputed sums of the values represented by each blob - * - Σqx_i = (sum_x - dim * min_x) / delta_x (sum of quantized values, derived from stored sum) - * - Σqy_i = (sum_y - dim * min_y) / delta_y + * where Σqx_i and Σqy_i are read directly from each blob's q_sum, and Σ(qx_i * qy_i) is the + * integer dot product of the quantized bytes. Every term is exact until the final combination, + * which happens in double (sq8::reconstructed_ip). * * For L2: - * ||x - y||² = sum_sq_x + sum_sq_y - 2 * IP(x, y) - * where sum_sq_x, sum_sq_y are precomputed sums of the squared represented values. + * ||x - y||² = Σ((min_x - min_y) + delta_x * qx_i - delta_y * qy_i)² + * = dim * c² + delta_x² * Σqx_i² + delta_y² * Σqy_i² + * + 2c * delta_x * Σqx_i - 2c * delta_y * Σqy_i + * - 2 * delta_x * delta_y * Σ(qx_i * qy_i), with c = min_x - min_y + * Deliberately not sum_sq_x + sum_sq_y - 2 * IP(x, y), which is algebraically identical but + * subtracts two large nearly-equal numbers: for vectors sharing a large offset the answer is + * smaller than one ULP of either term and is lost entirely. Expanding the difference first keeps + * the offset in c, formed once and exactly, so the summed terms scale with the spread instead. + * See sq8::reconstructed_l2_sqr. */ // Input types accepted by QuantPreprocessor. Opt-in via std::same_as so unrelated types // (e.g. integers, double, bfloat16) are rejected at the template head with a named constraint. @@ -269,10 +281,26 @@ class QuantPreprocessor : public PreprocessorInterface { float x_mean_ip = 0.0f; // only used for WithNorm auto [min_val, max_val] = find_min_max(input, x_mean_ip); - // Calculate scaling factor (typed as MetadataType because they end up as metadata). - const MetadataType diff = (max_val - min_val); - const MetadataType delta = (diff == 0.0f) ? MetadataType{1} : diff / MetadataType{255}; - const MetadataType inv_delta = MetadataType{1} / delta; + // The range is computed in double, and so is the per-element scaling below, because the + // FP32 arithmetic could leave the representable range for input that is entirely valid. + // [-FLT_MAX, +FLT_MAX] made max - min overflow to inf, then delta inf, inv_delta 0, and + // finally inf * 0 = NaN, whose conversion to an integer is undefined behaviour. Doubles + // cannot overflow for any pair of finite floats, so the whole class disappears rather than + // being special-cased. delta and inv_delta stay MetadataType, since delta is stored. + const double diff = static_cast(max_val) - static_cast(min_val); + const MetadataType delta = + (diff == 0.0) ? MetadataType{1} : static_cast(diff / 255.0); + const double inv_delta = 1.0 / static_cast(delta); + const double min_val_d = static_cast(min_val); + + // Scale one input value to its quantized byte. Kept in double for the reason above, and + // clamped so that a value outside [0, 255] can never reach the conversion: the arithmetic + // should not produce one, and if it ever does the result must be a wrong byte rather than + // undefined behaviour. + const auto to_byte = [inv_delta, min_val_d](float x) { + const double scaled = (static_cast(x) - min_val_d) * inv_delta; + return static_cast(std::round(std::clamp(scaled, 0.0, 255.0))); + }; // Sum the *quantized* bytes, not the input floats. A stored blob only ever describes its // reconstruction min + delta * a[i], so metadata derived from it has to describe that same @@ -298,10 +326,10 @@ class QuantPreprocessor : public PreprocessorInterface { // We know (input - min) => 0 // If min == max, all values are the same and should be quantized to 0. // reconstruction will yield the same original value for all vectors. - const uint32_t a0 = static_cast(std::round((x0 - min_val) * inv_delta)); - const uint32_t a1 = static_cast(std::round((x1 - min_val) * inv_delta)); - const uint32_t a2 = static_cast(std::round((x2 - min_val) * inv_delta)); - const uint32_t a3 = static_cast(std::round((x3 - min_val) * inv_delta)); + const uint32_t a0 = to_byte(x0); + const uint32_t a1 = to_byte(x1); + const uint32_t a2 = to_byte(x2); + const uint32_t a3 = to_byte(x3); quantized[i] = static_cast(a0); quantized[i + 1] = static_cast(a1); quantized[i + 2] = static_cast(a2); @@ -328,7 +356,7 @@ class QuantPreprocessor : public PreprocessorInterface { for (; i < this->dim; ++i) { const float x = transformed_value(input, i); - const uint32_t a = static_cast(std::round((x - min_val) * inv_delta)); + const uint32_t a = to_byte(x); quantized[i] = static_cast(a); q_sum += a; if constexpr (Metric == VecSimMetric_L2) { diff --git a/tests/unit/test_components.cpp b/tests/unit/test_components.cpp index 482623bba..764bc69cd 100644 --- a/tests/unit/test_components.cpp +++ b/tests/unit/test_components.cpp @@ -8,6 +8,8 @@ */ #include +#include +#include #include "gtest/gtest.h" #include "VecSim/vec_sim.h" #include "VecSim/spaces/computer/preprocessor_container.h" @@ -1200,6 +1202,45 @@ TEST(PreprocessorsTest, QuantizationAsymmetricAlignment) { } } +// Finite input whose range is not representable in FP32. max - min overflowed to inf, which made +// delta inf and inv_delta 0, and inf * 0 is NaN; converting that NaN to an integer is undefined +// behaviour, so this reproduced UB on AddVector for input the API accepts. Runs under UBSan in CI. +TEST(PreprocessorsTest, QuantizationHandlesNonRepresentableRange) { + std::shared_ptr allocator = VecSimAllocator::newVecsimAllocator(); + constexpr size_t dim = 4; + constexpr unsigned char alignment = 0; + constexpr float huge = std::numeric_limits::max(); + float original_blob[dim] = {-huge, 0.0f, huge, 1.0f}; + + auto quant_preprocessor = + new (allocator) QuantPreprocessor(allocator, dim); + + void *storage_blob = nullptr; + void *query_blob = nullptr; + size_t storage_blob_size = dim * sizeof(float); + size_t query_blob_size = dim * sizeof(float); + quant_preprocessor->preprocess(original_blob, storage_blob, query_blob, storage_blob_size, + query_blob_size, alignment, alignment); + ASSERT_NE(storage_blob, nullptr); + + // The extremes must land on the ends of the range and everything must stay in [0, 255]. The + // middle values collapse to the same byte, which is expected: one byte cannot resolve a range + // this wide, and losing resolution is the correct outcome where UB was not. + const auto *quantized = static_cast(storage_blob); + EXPECT_EQ(quantized[0], 0); + EXPECT_EQ(quantized[2], 255); + + // Metadata must be finite, since the kernels derive every term from it. + const auto *metadata = quantized + dim; + EXPECT_TRUE(std::isfinite(load_unaligned(metadata + sq8::MIN_VAL * sizeof(float)))); + EXPECT_TRUE(std::isfinite(load_unaligned(metadata + sq8::DELTA * sizeof(float)))); + + allocator->free_allocation(storage_blob); + if (query_blob != storage_blob) { + allocator->free_allocation(query_blob); + } +} + // Test edge case where all entries are equal TEST(PreprocessorsTest, QuantizationTestAllEntriesEqual) { std::shared_ptr allocator = VecSimAllocator::newVecsimAllocator(); diff --git a/tests/unit/test_spaces.cpp b/tests/unit/test_spaces.cpp index c599afc64..c7fe454fb 100644 --- a/tests/unit/test_spaces.cpp +++ b/tests/unit/test_spaces.cpp @@ -371,6 +371,100 @@ TEST_F(SpacesTest, SQ8_FP32_l2sqr_no_optimization_func_test) { ASSERT_NEAR(dist, baseline, 0.01) << "SQ8_FP32_L2Sqr failed to match expected distance"; } +// The three cases below are the defects this formulation exists to fix. They are written against +// an FP64 reference over the reconstructed vectors, which is the quantity the kernels are trying to +// compute, rather than against numbers recorded from a previous run. +namespace { +// ||x_r - y_r||^2 in FP64, where each operand is read back from its own quantized blob. +double ReferenceL2SqrFromBlobs(const uint8_t *a, const uint8_t *b, size_t dim) { + const double min_a = sq8::min_val_of(a, dim), delta_a = sq8::delta_of(a, dim); + const double min_b = sq8::min_val_of(b, dim), delta_b = sq8::delta_of(b, dim); + double acc = 0.0; + for (size_t i = 0; i < dim; i++) { + const double d = (min_a + delta_a * a[i]) - (min_b + delta_b * b[i]); + acc += d * d; + } + return acc; +} +} // namespace + +// A large common offset used to cancel the answer away entirely. Quantization is exact for this +// input (q = [0, 255]), so a wrong result is the summation form's fault, not the 8 bits'. +TEST_F(SpacesTest, SQ8_SQ8_L2_survives_large_common_offset) { + constexpr size_t dim = 2; + const std::vector x = {100000.0f, 100008.0f}; + const std::vector y = {100000.0f, 100000.0f}; + + const size_t blob_size = + dim * sizeof(uint8_t) + sq8::storage_metadata_count() * sizeof(float); + std::vector qx(blob_size), qy(blob_size); + test_utils::quantize_float_vec_to_sq8_with_metadata(x.data(), dim, qx.data()); + test_utils::quantize_float_vec_to_sq8_with_metadata(y.data(), dim, qy.data()); + + const double expected = ReferenceL2SqrFromBlobs(qx.data(), qy.data(), dim); + ASSERT_NEAR(expected, 64.0, 1e-6) << "the reference itself should see 64"; + + // Previously returned 0: ||x||^2 and 2*IP are both about 4e10, one FP32 step there is ~4768, + // and the answer is 64. + EXPECT_NEAR(SQ8_SQ8_L2Sqr(qx.data(), qy.data(), dim), expected, 1e-3); +} + +// A query that does not land on a quantization grid point. The stored sums used to be taken over +// the original input while the inner product described the reconstruction, mixing two vectors into +// one formula, which drove the result negative and made a radius-0 range query return a match. +TEST_F(SpacesTest, SQ8_FP32_L2_is_non_negative_for_off_grid_query) { + constexpr size_t dim = 3; + const std::vector stored = {0.0f, 0.25f, 1.0f}; + + const size_t blob_size = + dim * sizeof(uint8_t) + sq8::storage_metadata_count() * sizeof(float); + std::vector qs(blob_size); + test_utils::quantize_float_vec_to_sq8_with_metadata(stored.data(), dim, qs.data()); + + std::vector query(dim + sq8::query_metadata_count()); + query[0] = 0.0f; + query[1] = 0.2501f; + query[2] = 1.0f; + test_utils::preprocess_sq8_fp32_query(query.data(), dim); + + // Reference: reconstruction of the stored blob against the exact query values. + const double min_s = sq8::min_val_of(qs.data(), dim), delta_s = sq8::delta_of(qs.data(), dim); + double expected = 0.0; + for (size_t i = 0; i < dim; i++) { + const double d = (min_s + delta_s * qs[i]) - static_cast(query[i]); + expected += d * d; + } + + const float dist = SQ8_FP32_L2Sqr(qs.data(), query.data(), dim); + EXPECT_GE(dist, 0.0f) << "a squared distance must never be negative"; + EXPECT_NEAR(dist, expected, 1e-6); +} + +// High dimension, near-identical vectors: the regime where combining the exact integer sums in FP32 +// still loses the answer, because the three large sums cancel. +TEST_F(SpacesTest, SQ8_SQ8_L2_matches_fp64_reference_at_high_dim) { + constexpr size_t dim = 16384; + std::vector x(dim), y(dim); + std::mt19937 gen(4242); + std::uniform_real_distribution dist_gen(0.0f, 1.0f); + for (size_t i = 0; i < dim; i++) { + x[i] = dist_gen(gen); + y[i] = x[i] + (i % 7 == 0 ? 0.0005f : 0.0f); + } + + const size_t blob_size = + dim * sizeof(uint8_t) + sq8::storage_metadata_count() * sizeof(float); + std::vector qx(blob_size), qy(blob_size); + test_utils::quantize_float_vec_to_sq8_with_metadata(x.data(), dim, qx.data()); + test_utils::quantize_float_vec_to_sq8_with_metadata(y.data(), dim, qy.data()); + + const double expected = ReferenceL2SqrFromBlobs(qx.data(), qy.data(), dim); + ASSERT_GT(expected, 0.0); + const double got = SQ8_SQ8_L2Sqr(qx.data(), qy.data(), dim); + EXPECT_LT(std::abs(got - expected) / expected, 1e-5) + << "got " << got << " expected " << expected; +} + TEST_F(SpacesTest, SQ8_FP32_odd_dim_unaligned_metadata_test) { for (const size_t dim : {1UL, 5UL, 7UL, 15UL}) { const size_t query_size = dim + sq8::query_metadata_count(); From a82713b01f08bea821fce3198ef6b61f77ed6a5f Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Mon, 10 Aug 2026 17:30:45 +0300 Subject: [PATCH 06/17] Loosen a sanity bound to the quantization residue The large-offset regression asserted that the FP64 reference equals 64 within 1e-6 before comparing the kernel against it. The reference is 64.0000076: reconstructing through FP32 min and delta does not land exactly on 100008, so that residue is quantization error and belongs in the bound. The assertion that matters, kernel against reference, passed unchanged. Co-Authored-By: Claude Opus 5 --- tests/unit/test_spaces.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/tests/unit/test_spaces.cpp b/tests/unit/test_spaces.cpp index c7fe454fb..0c7ef4d70 100644 --- a/tests/unit/test_spaces.cpp +++ b/tests/unit/test_spaces.cpp @@ -401,8 +401,12 @@ TEST_F(SpacesTest, SQ8_SQ8_L2_survives_large_common_offset) { test_utils::quantize_float_vec_to_sq8_with_metadata(x.data(), dim, qx.data()); test_utils::quantize_float_vec_to_sq8_with_metadata(y.data(), dim, qy.data()); + // 64 is the distance between the inputs; the reference sees 64.0000076 because reconstructing + // through FP32 min and delta does not land exactly on 100008. That residue is quantization, + // which is expected and tiny. What matters is that the kernel agrees with the reference rather + // than returning 0. const double expected = ReferenceL2SqrFromBlobs(qx.data(), qy.data(), dim); - ASSERT_NEAR(expected, 64.0, 1e-6) << "the reference itself should see 64"; + ASSERT_NEAR(expected, 64.0, 1e-4) << "the reference itself should see 64"; // Previously returned 0: ||x||^2 and 2*IP are both about 4e10, one FP32 step there is ~4768, // and the answer is 64. From b711c0c9578920bf9e71dbd823364f8b0c59620c Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 12 Aug 2026 14:40:46 +0300 Subject: [PATCH 07/17] fix(uint8): widen and unsign the integer accumulators The uint8 kernels accumulate products of bytes, so the total reaches 255*255*dim = 65025*dim. Three paths could not hold that: * IP.cpp / L2.cpp: ret_t for a 1-byte element type was int, so the scalar UINT8_InnerProduct and UINT8_L2Sqr hit signed-overflow UB from dimension 33,026, while the comment claimed support to 2^16. ret_t is now 64-bit for every element type. Keeping it signed means the "1 - ip" in the wrappers stays signed arithmetic and cannot underflow, and the int8 paths are unaffected in behaviour. * L2_AVX512F_BW_VL_VNNI_UINT8: the horizontal reduce was read back as a signed int, so the distance went negative from dimension 33,026. It is now read as uint32_t, which costs no instructions. * L2_NEON_UINT8: same defect, int32_t receiving vaddvq_u32. Now uint32_t. L2_NEON_DOTPROD_UINT8 and L2_SVE_UINT8 were already unsigned and are unchanged. This is the same defect class as MOD-17527, which this branch already fixes on the inner product side; these are the L2 and scalar paths it did not cover. Also fixes the uint8 spaces benchmark fixture, which paired new[] with delete and stored the trailing norms through unaligned float casts, so that measurements taken from it are trustworthy. --- src/VecSim/spaces/IP/IP.cpp | 11 ++++++----- src/VecSim/spaces/L2/L2.cpp | 2 +- src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_UINT8.h | 4 +++- src/VecSim/spaces/L2/L2_NEON_UINT8.h | 4 +++- tests/benchmark/spaces_benchmarks/bm_spaces_uint8.cpp | 11 +++++++---- 5 files changed, 20 insertions(+), 12 deletions(-) diff --git a/src/VecSim/spaces/IP/IP.cpp b/src/VecSim/spaces/IP/IP.cpp index 6894f7c60..b78ee4fe2 100644 --- a/src/VecSim/spaces/IP/IP.cpp +++ b/src/VecSim/spaces/IP/IP.cpp @@ -228,11 +228,12 @@ float FP16_InnerProduct(const void *pVect1, const void *pVect2, size_t dimension } // Return type for the inner product functions. -// The type should be able to hold `dimension * MAX_VAL(int_elem_t) * MAX_VAL(int_elem_t)`. -// To support dimension up to 2^16, we need the difference between the type and int_elem_t to be at -// least 2 bytes. We assert that in the implementation. +// The type must hold `dimension * MAX_VAL(int_elem_t) * MAX_VAL(int_elem_t)`. For uint8 that is +// 65025 * dimension, which overflows a 32-bit int from dimension 33,026, so this is 64-bit for +// every element type. Keeping it signed also means the `1 - ip` in the wrappers below stays +// signed arithmetic and cannot underflow. We assert the width in the implementation. template -using ret_t = std::conditional_t; +using ret_t = long long; template static inline ret_t @@ -263,7 +264,7 @@ float INT8_Cosine(const void *pVect1v, const void *pVect2v, size_t dimension) { float UINT8_InnerProduct(const void *pVect1v, const void *pVect2v, size_t dimension) { const auto *pVect1 = static_cast(pVect1v); const auto *pVect2 = static_cast(pVect2v); - return 1 - INTEGER_InnerProductImp(pVect1, pVect2, dimension); + return 1.0f - static_cast(INTEGER_InnerProductImp(pVect1, pVect2, dimension)); } float UINT8_Cosine(const void *pVect1v, const void *pVect2v, size_t dimension) { diff --git a/src/VecSim/spaces/L2/L2.cpp b/src/VecSim/spaces/L2/L2.cpp index 3d3f682d5..37f371174 100644 --- a/src/VecSim/spaces/L2/L2.cpp +++ b/src/VecSim/spaces/L2/L2.cpp @@ -144,7 +144,7 @@ float FP16_L2Sqr(const void *pVect1, const void *pVect2, size_t dimension) { // To support dimension up to 2^16, we need the difference between the type and int_elem_t to be at // least 2 bytes. We assert that in the implementation. template -using ret_t = std::conditional_t; +using ret_t = long long; // Difference type for the L2 functions. // The type should be able to hold `MIN_VAL(int_elem_t)-MAX_VAL(int_elem_t)`, and should be signed diff --git a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_UINT8.h b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_UINT8.h index 350b759ea..f906e7dc4 100644 --- a/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_UINT8.h +++ b/src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_UINT8.h @@ -92,5 +92,7 @@ float UINT8_L2SqrSIMD64_AVX512F_BW_VL_VNNI(const void *pVect1v, const void *pVec } while (pVect1 < pEnd1); } - return _mm512_reduce_add_epi32(sum); + // The lanes hold sums of squared byte differences, so the horizontal total is unsigned and + // reaches 255*255*dim. Reading it as a signed int wrapped it negative from dimension 33,026. + return static_cast(static_cast(_mm512_reduce_add_epi32(sum))); } diff --git a/src/VecSim/spaces/L2/L2_NEON_UINT8.h b/src/VecSim/spaces/L2/L2_NEON_UINT8.h index aa3769867..bd280a750 100644 --- a/src/VecSim/spaces/L2/L2_NEON_UINT8.h +++ b/src/VecSim/spaces/L2/L2_NEON_UINT8.h @@ -126,7 +126,9 @@ float UINT8_L2SqrSIMD16_NEON(const void *pVect1v, const void *pVect2v, size_t di total_sum = vaddq_u32(total_sum, sum3); // Horizontal sum of the 4 elements in the combined sum register - int32_t result = vaddvq_u32(total_sum); + // Unsigned: the total is a sum of squared byte differences, reaching 255*255*dim. As a signed + // int32 this wrapped negative from dimension 33,026. + uint32_t result = vaddvq_u32(total_sum); // Return the L2 squared distance as a float return static_cast(result); diff --git a/tests/benchmark/spaces_benchmarks/bm_spaces_uint8.cpp b/tests/benchmark/spaces_benchmarks/bm_spaces_uint8.cpp index 602fff719..33f819936 100644 --- a/tests/benchmark/spaces_benchmarks/bm_spaces_uint8.cpp +++ b/tests/benchmark/spaces_benchmarks/bm_spaces_uint8.cpp @@ -31,12 +31,15 @@ class BM_VecSimSpaces_Integers_UINT8 : public benchmark::Fixture { test_utils::populate_uint8_vec(v2, dim, 1234); // Store the norm in the extra space for cosine calculations - *(float *)(v1 + dim) = test_utils::integral_compute_norm(v1, dim); - *(float *)(v2 + dim) = test_utils::integral_compute_norm(v2, dim); + // memcpy because v1 + dim is not guaranteed to be 4-byte aligned for arbitrary dim. + const float norm1 = test_utils::integral_compute_norm(v1, dim); + const float norm2 = test_utils::integral_compute_norm(v2, dim); + memcpy(v1 + dim, &norm1, sizeof(norm1)); + memcpy(v2 + dim, &norm2, sizeof(norm2)); } void TearDown(const ::benchmark::State &state) { - delete v1; - delete v2; + delete[] v1; + delete[] v2; } }; From 43250bf6f7ec37c13ea258f36fb66aafa4cfec27 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 12 Aug 2026 14:41:20 +0300 Subject: [PATCH 08/17] refactor(sq8): drop MAX_EXACT_DIM The constant was not load-bearing: it appeared only in its own definition and in one debug assert, so release builds behaved identically with or without it, and it told the reader nothing that was true. Its comment claimed 33,026 was the uint32 exactness limit for sum(a[i]^2), but 65025 * dim stays exact in uint32 through dimension 66,051; 33,026 is where a *signed* 32-bit accumulator would have wrapped, which is a bound this code no longer has anywhere. For the record, the real per-field ceilings are: q_sum_squares (L2 only) exact through dim 66,051, q_sum through 16,843,009, and the AVX512 per-lane accumulation through roughly 528,000. The first is the binding one, and it is far above any dimension the library is used at. --- src/VecSim/spaces/computer/preprocessors.h | 3 --- src/VecSim/types/sq8.h | 3 --- 2 files changed, 6 deletions(-) diff --git a/src/VecSim/spaces/computer/preprocessors.h b/src/VecSim/spaces/computer/preprocessors.h index a670d52f2..be77a0345 100644 --- a/src/VecSim/spaces/computer/preprocessors.h +++ b/src/VecSim/spaces/computer/preprocessors.h @@ -364,9 +364,6 @@ class QuantPreprocessor : public PreprocessorInterface { } } - // sum(a[i]^2) is exact in uint32 only up to sq8::MAX_EXACT_DIM; past it this wraps and the - // L2 distance silently degrades. Callers that build quantized indexes must fence dim. - assert(this->dim <= sq8::MAX_EXACT_DIM && "SQ8 dimension exceeds exact metadata range"); // Metadata is a packed run of 4-byte fields: min and delta as FP32, the quantized sums as // uint32, and x_mean_ip (FP32) last for IP with norm. memcpy throughout because the offset diff --git a/src/VecSim/types/sq8.h b/src/VecSim/types/sq8.h index 1721fde26..c546e694d 100644 --- a/src/VecSim/types/sq8.h +++ b/src/VecSim/types/sq8.h @@ -38,9 +38,6 @@ struct sq8 { Q_SUM_SQUARES = 3 // uint32: sum(a[i]^2). Only for L2 }; - // sum(a[i]^2) <= 255^2 * dim, which is exactly representable in uint32 up to this dimension. - // Past it the accumulation wraps, so callers that build quantized indexes must fence dim. - static constexpr size_t MAX_EXACT_DIM = 33026; enum QueryMetadataIndex : size_t { SUM_QUERY = 0, From c090e1514bde53c36ef5efe658203606bece5ff7 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 12 Aug 2026 14:41:20 +0300 Subject: [PATCH 09/17] perf(sq8): avoid a libm round call per element when quantizing std::round compiles to an out-of-line call to round() here, because this translation unit is built at the x86-64 baseline and so has no roundsd available; nm confirms the symbol and the object contains no rounding instruction. That call ran once per element. scaled is already clamped to [0, 255] and therefore non-negative, so adding 0.5 and truncating is std::round, and compiles to a single cvttsd2si. Measured on Ice Lake-SP, this makes the quantization path about 26% faster than before this branch, which also absorbs the ~18% this branch had added to it. The two forms differ only where x + 0.5 itself double-rounds, at scaled == 0.49999999999999994, where this yields byte 1 and std::round yields 0. No divergence at any of the 255 .5 boundaries, across 2M uniform random values, or across 3M realistic (x - min) * inv_delta draws. Note that lrint and nearbyint are NOT substitutes: they round ties to even, which would move 178.5 to 178 and change bytes in common cases. --- src/VecSim/spaces/computer/preprocessors.h | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/VecSim/spaces/computer/preprocessors.h b/src/VecSim/spaces/computer/preprocessors.h index be77a0345..654192c53 100644 --- a/src/VecSim/spaces/computer/preprocessors.h +++ b/src/VecSim/spaces/computer/preprocessors.h @@ -299,7 +299,11 @@ class QuantPreprocessor : public PreprocessorInterface { // undefined behaviour. const auto to_byte = [inv_delta, min_val_d](float x) { const double scaled = (static_cast(x) - min_val_d) * inv_delta; - return static_cast(std::round(std::clamp(scaled, 0.0, 255.0))); + // Clamped first, so the value is non-negative and adding 0.5 before truncating is + // std::round without the library call. This translation unit is compiled at the + // x86-64 baseline, where std::round is an out-of-line call to round() rather than a + // roundsd instruction, and this runs once per element. + return static_cast(std::clamp(scaled, 0.0, 255.0) + 0.5); }; // Sum the *quantized* bytes, not the input floats. A stored blob only ever describes its From f6780aa37f7fe82bbd481c345209e219b8ee4da1 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 12 Aug 2026 14:42:00 +0300 Subject: [PATCH 10/17] fix(sq8): make the symmetric L2 cancellation exact reconstructed_l2_sqr returned a slightly negative squared distance for a blob against itself: about -6e-14 at dimension 64 and -1.7e-12 at 512 on the AVX512 kernel. The six terms cancel exactly in real arithmetic, but the compiler contracts some of the products into fused multiply-adds and not others, so they no longer round identically and no longer cancel. Harmless to ranking, but a caller taking a square root of it gets NaN, and a function named L2Sqr should not return a negative number. Regrouping the quadratic part fixes it structurally: d1^2*S1 + d2^2*S2 - 2*d1*d2*Q == d1*d2*(S1 + S2 - 2Q) + (d1 - d2)*(d1*S1 - d2*S2) S1 + S2 - 2Q is sum((a[i] - b[i])^2), an exact non-negative integer, so for equal deltas the answer is a non-negative integer scaled by delta^2 and cannot be negative, and for identical blobs both factors are exactly zero. Multiplying by an exact zero is exact, so contraction cannot break it. Measured against a long double reference with -mfma: self-distance exactly zero at dimensions 3 through 4096, mean relative error on random pairs 5.6e-17 versus 2.9e-16 before, and in the regime this branch exists to fix, near-duplicates sharing a large offset, exactly zero error versus 1e-9 before. Deliberately not fixed by clamping at zero. A clamp would have turned the -201.5 the old metadata produced into a clean 0 and hidden the very defect this branch fixes. It also shortens the dependency chain, since the integer combination is independent of the delta terms. --- src/VecSim/types/sq8.h | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/src/VecSim/types/sq8.h b/src/VecSim/types/sq8.h index c546e694d..41a088a3b 100644 --- a/src/VecSim/types/sq8.h +++ b/src/VecSim/types/sq8.h @@ -142,13 +142,26 @@ struct sq8 { uint32_t q_sum_squares1, double min2, double delta2, uint32_t q_sum2, uint32_t q_sum_squares2, uint64_t q_dot, size_t dim) { + // The quadratic part is regrouped so that the integer combination is formed first: + // d1^2*S1 + d2^2*S2 - 2*d1*d2*Q == d1*d2*(S1 + S2 - 2Q) + (d1 - d2)*(d1*S1 - d2*S2) + // S1 + S2 - 2Q is sum((a[i] - b[i])^2), an exact non-negative integer. Two consequences: + // for two blobs sharing a delta the answer is delta^2 times a non-negative integer and so + // cannot come out negative, and for a blob against itself both factors are exactly zero, so + // the distance is exactly zero. Written as six independent floating point terms it was not: + // the compiler contracts some of them into fused multiply-adds and not others, so the + // products stop rounding identically and stop cancelling, which left about -1.7e-12 at + // dimension 512. Note this is deliberately not fixed by clamping at zero, which would also + // have hidden the much larger negative values the old metadata produced. + const int64_t sq_diff_sum = static_cast(q_sum_squares1) + + static_cast(q_sum_squares2) - + 2 * static_cast(q_dot); + const double quadratic = delta1 * delta2 * static_cast(sq_diff_sum) + + (delta1 - delta2) * (delta1 * static_cast(q_sum_squares1) - + delta2 * static_cast(q_sum_squares2)); const double c = min1 - min2; return static_cast(dim) * c * c + - delta1 * delta1 * static_cast(q_sum_squares1) + - delta2 * delta2 * static_cast(q_sum_squares2) + 2.0 * c * delta1 * static_cast(q_sum1) - - 2.0 * c * delta2 * static_cast(q_sum2) - - 2.0 * delta1 * delta2 * static_cast(q_dot); + 2.0 * c * delta2 * static_cast(q_sum2) + quadratic; } }; From 62929137e0cf2e61392c5a863ce4116915100df8 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 12 Aug 2026 14:43:06 +0300 Subject: [PATCH 11/17] perf(sq8): force-inline the reconstruction helpers GCC outlines reconstructed_l2_sqr, plausibly because the 64 residual instantiations of each SQ8 kernel all call it. The cost is not just the call: in the AVX512 kernel it also sets up a realigned 64-byte stack frame and a vzeroupper inside what is otherwise a leaf SIMD function. SQ8_SQ8_L2SqrSIMD64_AVX512F_BW_VL_VNNI<0> went from 15 instructions to 54 for this reason. Marking the six helpers always_inline removes all 96 out-of-line calls in that translation unit and 33 of the 48 stack realignments, and recovers 24% to 30% of the per-call regression this branch introduced (median +2.15 ns to +1.44 ns on spaces_sq8_sq8, replicated across two runs per side on Ice Lake-SP). No case regressed against the unpatched branch. The cost is about 11 extra instructions per instantiation, since a 24-instruction body is now duplicated rather than shared. --- src/VecSim/types/sq8.h | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/VecSim/types/sq8.h b/src/VecSim/types/sq8.h index 41a088a3b..70b710640 100644 --- a/src/VecSim/types/sq8.h +++ b/src/VecSim/types/sq8.h @@ -38,7 +38,6 @@ struct sq8 { Q_SUM_SQUARES = 3 // uint32: sum(a[i]^2). Only for L2 }; - enum QueryMetadataIndex : size_t { SUM_QUERY = 0, SUM_SQUARES_QUERY = 1 // Only for L2 @@ -102,11 +101,13 @@ struct sq8 { // exact, and the point is to keep them exact through the combination. // sum(x_r[i]) where x_r[i] = min + delta * a[i]. + [[gnu::always_inline]] static double reconstructed_sum(double min_val, double delta, uint32_t q_sum, size_t dim) { return static_cast(dim) * min_val + delta * static_cast(q_sum); } // sum(x_r[i]^2), expanded so that no term is ever formed from a rounded input. + [[gnu::always_inline]] static double reconstructed_sum_squares(double min_val, double delta, uint32_t q_sum, uint32_t q_sum_squares, size_t dim) { return static_cast(dim) * min_val * min_val + @@ -116,15 +117,18 @@ struct sq8 { // ||x_r||^2 for one stored blob. Two entry points because some kernels have the blob pointer // and some have already formed the metadata pointer. + [[gnu::always_inline]] static double reconstructed_sum_squares_at(const value_type *meta, size_t dim) { return reconstructed_sum_squares(min_val_at(meta), delta_at(meta), q_sum_at(meta), q_sum_squares_at(meta), dim); } + [[gnu::always_inline]] static double reconstructed_sum_squares_of(const value_type *blob, size_t dim) { return reconstructed_sum_squares_at(blob + dim, dim); } // IP(x_r, y_r) for two quantized blobs, from their exact integer sums and integer dot product. + [[gnu::always_inline]] static double reconstructed_ip(double min1, double delta1, uint32_t q_sum1, double min2, double delta2, uint32_t q_sum2, uint64_t q_dot, size_t dim) { return static_cast(dim) * min1 * min2 + @@ -138,6 +142,7 @@ struct sq8 { // letting it inflate the summed terms and cancel away the answer. That is why this is not // written as reconstructed_sum_squares + reconstructed_sum_squares - 2 * reconstructed_ip: // algebraically identical, numerically not. + [[gnu::always_inline]] static double reconstructed_l2_sqr(double min1, double delta1, uint32_t q_sum1, uint32_t q_sum_squares1, double min2, double delta2, uint32_t q_sum2, uint32_t q_sum_squares2, uint64_t q_dot, @@ -159,8 +164,7 @@ struct sq8 { (delta1 - delta2) * (delta1 * static_cast(q_sum_squares1) - delta2 * static_cast(q_sum_squares2)); const double c = min1 - min2; - return static_cast(dim) * c * c + - 2.0 * c * delta1 * static_cast(q_sum1) - + return static_cast(dim) * c * c + 2.0 * c * delta1 * static_cast(q_sum1) - 2.0 * c * delta2 * static_cast(q_sum2) + quadratic; } }; From 983c742db4e88475aff66a010b350842664cc79a Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 12 Aug 2026 14:45:11 +0300 Subject: [PATCH 12/17] test(sq8): add regression tests for the exact-metadata contract Six tests, each of which fails on the commit before this series and passes after it. They were run on both revisions to confirm that. They use only the quantizer test helper, sq8::MIN_VAL, sq8::DELTA and the public kernels, so the same source compiles on either side. That is deliberate: a test written against Q_SUM or sq8::reconstructed_* would not compile before the change and so could not demonstrate anything. * self-distance is exactly zero, scalar and dispatched, five dimensions. Was nonzero at every dimension, sign varying with the rounding of individual elements. * the same property with hand-checkable numbers, where the old result is exactly -201.5 for a vector against itself. * two inputs that quantize to identical bytes produce identical distances. The distance must be a function of the stored blob and nothing else; it differed by 2.1e-3 for L2 and 1.5e-3 for IP. * the inner product equals the dot product of the reconstructions, which was off by 2.3e-3 relative. * plain uint8 L2 and inner product stay exact past INT_MAX, at dimensions 33,026 and 40,000 with worst-case all-255 bytes. The scalar path was UB there and the AVX512 and NEON L2 reduces returned negative distances. The dimensions are chosen so that both the scalar kernels and the SIMD kernels are exercised: 3, 5 and 8 dispatch to the naive implementation and 64 and above to AVX512F_BW_VL_VNNI on an Ice Lake host. The third test keeps min nonzero on purpose. With min == 0 the old inner product formula collapses to delta_x * delta_y * q_dot, the input-derived sums drop out, and the inner product half of the test cannot detect anything. --- tests/unit/test_spaces.cpp | 166 +++++++++++++++++++++++++++++++++++++ 1 file changed, 166 insertions(+) diff --git a/tests/unit/test_spaces.cpp b/tests/unit/test_spaces.cpp index 0c7ef4d70..d6d39cdf4 100644 --- a/tests/unit/test_spaces.cpp +++ b/tests/unit/test_spaces.cpp @@ -4869,3 +4869,169 @@ TEST_F(SpacesTest, SQ8_SQ8_DispatcherAlignmentHints) { check("Cosine", &spaces::Cosine_SQ8_SQ8_GetDistFunc); } #endif // CPU_FEATURES_ARCH_X86_64 + +// =================================================================================== +// Regression tests for the exact-metadata contract. Every one of these fails on the +// commit before this series and passes after it. They deliberately use only the +// quantizer helper, sq8::MIN_VAL, sq8::DELTA and the public kernels, so the same +// source compiles on both revisions and the before/after comparison is meaningful. +// =================================================================================== + +#include + +namespace { + +// min_val and delta occupy metadata slots 0 and 1 as FP32 in every revision. +float sq8_meta_float(const uint8_t *blob, size_t dim, size_t slot) { + float v; + std::memcpy(&v, blob + dim + slot * sizeof(float), sizeof(v)); + return v; +} + +// FP64 dot product of the two reconstructions, derived only from the bytes, min and delta. +// This is the quantity every SQ8-to-SQ8 kernel is trying to compute. +double sq8_reference_dot(const uint8_t *a, const uint8_t *b, size_t dim) { + const double mn_a = sq8_meta_float(a, dim, sq8::MIN_VAL); + const double d_a = sq8_meta_float(a, dim, sq8::DELTA); + const double mn_b = sq8_meta_float(b, dim, sq8::MIN_VAL); + const double d_b = sq8_meta_float(b, dim, sq8::DELTA); + double acc = 0.0; + for (size_t i = 0; i < dim; i++) { + acc += (mn_a + d_a * a[i]) * (mn_b + d_b * b[i]); + } + return acc; +} + +std::vector sq8_quantize_l2(const std::vector &v) { + const size_t dim = v.size(); + std::vector blob(dim * sizeof(uint8_t) + + sq8::storage_metadata_count() * sizeof(float)); + test_utils::quantize_float_vec_to_sq8_with_metadata(v.data(), dim, blob.data()); + return blob; +} + +} // namespace + +// A blob against itself must be exactly zero, on every kernel. Both factors of the regrouped +// quadratic term are exactly zero for identical blobs, and multiplying by an exact zero is exact, +// so this holds regardless of how the compiler contracts the surrounding arithmetic. Before the +// metadata change the residue was 4 * sum((x_r[i] - min) * rounding_error[i]), nonzero and of +// arbitrary sign. +TEST_F(SpacesTest, SQ8_SQ8_L2_self_distance_is_exactly_zero) { + std::mt19937 gen(20260812); + std::uniform_real_distribution value_gen(-3.0f, 5.0f); + + for (const size_t dim : {3UL, 5UL, 8UL, 64UL, 512UL}) { + std::vector x(dim); + for (size_t i = 0; i < dim; i++) { + x[i] = value_gen(gen); + } + const std::vector q = sq8_quantize_l2(x); + + EXPECT_EQ(SQ8_SQ8_L2Sqr(q.data(), q.data(), dim), 0.0f) << "scalar kernel, dim " << dim; + + // Also exercise whatever this CPU dispatches to. + unsigned char alignment = 0; + auto dispatched = L2_SQ8_SQ8_GetDistFunc(dim, &alignment, nullptr); + EXPECT_EQ(dispatched(q.data(), q.data(), dim), 0.0f) << "dispatched kernel, dim " << dim; + } +} + +// The same property with numbers small enough to check by hand. x = [0, 100.5, 255] gives min 0 and +// delta exactly 1, so the bytes are [0, 101, 255] and the stored vector is [0, 101, 255]. Summing +// the input gave x_sum_squares 75125.25 while the dot product of the bytes gives 75226, and the old +// formula used the sum of squares twice, so it returned 75125.25 + 75125.25 - 2 * 75226 = -201.5. +TEST_F(SpacesTest, SQ8_SQ8_L2_self_distance_whiteboard_case) { + const std::vector x = {0.0f, 100.5f, 255.0f}; + const std::vector q = sq8_quantize_l2(x); + + ASSERT_EQ(q[0], 0); + ASSERT_EQ(q[1], 101) << "expected 100.5 to round up to byte 101"; + ASSERT_EQ(q[2], 255); + + const float dist = SQ8_SQ8_L2Sqr(q.data(), q.data(), 3); + EXPECT_GE(dist, 0.0f) << "a squared distance must never be negative"; + EXPECT_EQ(dist, 0.0f); +} + +// Two inputs that quantize to the same bytes are the same stored vector, so they must answer a +// query identically. The distance has to be a function of the stored blob and nothing else. +TEST_F(SpacesTest, SQ8_SQ8_distance_ignores_sub_quantum_input_changes) { + constexpr size_t dim = 4; + // Same min and max, and both interior values perturbed well under half a delta and away from a + // rounding boundary, so both inputs land on the identical byte string. min is nonzero: with + // min == 0 the old inner product formula collapses to delta_x * delta_y * q_dot and the + // input-derived sums drop out, which would make the IP half of this test inert. + const std::vector x = {2.0f, 2.4f, 2.72f, 3.0f}; + const std::vector x_perturbed = {2.0f, 2.4005f, 2.7205f, 3.0f}; + const std::vector y = {1.5f, 1.6f, 2.4f, 2.5f}; + + const std::vector qx = sq8_quantize_l2(x); + const std::vector qxp = sq8_quantize_l2(x_perturbed); + const std::vector qy = sq8_quantize_l2(y); + + ASSERT_TRUE(std::equal(qx.begin(), qx.begin() + dim, qxp.begin())) + << "inputs no longer quantize identically; pick new values"; + ASSERT_EQ(sq8_meta_float(qx.data(), dim, sq8::MIN_VAL), + sq8_meta_float(qxp.data(), dim, sq8::MIN_VAL)); + ASSERT_EQ(sq8_meta_float(qx.data(), dim, sq8::DELTA), + sq8_meta_float(qxp.data(), dim, sq8::DELTA)); + + EXPECT_EQ(SQ8_SQ8_L2Sqr(qx.data(), qy.data(), dim), SQ8_SQ8_L2Sqr(qxp.data(), qy.data(), dim)); + EXPECT_EQ(SQ8_SQ8_InnerProduct(qx.data(), qy.data(), dim), + SQ8_SQ8_InnerProduct(qxp.data(), qy.data(), dim)); +} + +// The inner product must equal the dot product of the two reconstructions. +TEST_F(SpacesTest, SQ8_SQ8_IP_matches_reconstruction_dot) { + const std::vector> inputs = { + {-1.0f, 0.3f, 2.0f}, + {-2.5f, -1.1f, 0.4f, 0.9f, 3.3f, 7.7f, 8.0f, 9.9f}, + }; + for (const auto &x : inputs) { + const size_t dim = x.size(); + const std::vector q = sq8_quantize_l2(x); + const double expected = sq8_reference_dot(q.data(), q.data(), dim); + + // SQ8_SQ8_InnerProduct returns 1 - IP. + const double got = 1.0 - static_cast(SQ8_SQ8_InnerProduct(q.data(), q.data(), dim)); + EXPECT_LT(std::abs(got - expected) / expected, 1e-5) + << "dim " << dim << ": got " << got << " expected " << expected; + } +} + +// Plain uint8, not SQ8. The accumulated total is 255 * 255 * dim, which passes INT_MAX from +// dimension 33,026: the scalar path was signed-overflow UB there and the AVX512 and NEON L2 +// reduces read their unsigned total back as a signed int and went negative. All-255 bytes are the +// worst case and make the expected value an exact integer. +TEST_F(SpacesTest, UINT8_L2Sqr_and_InnerProduct_are_exact_past_int32) { + for (const size_t dim : {33026UL, 40000UL}) { + std::vector v1(dim + sizeof(float), 255); + std::vector v2(dim + sizeof(float), 0); + + // L2 between all-255 and all-0 is 255^2 * dim. + const double expected_l2 = 255.0 * 255.0 * static_cast(dim); + const float l2 = UINT8_L2Sqr(v1.data(), v2.data(), dim); + EXPECT_GT(l2, 0.0f) << "dim " << dim << ": squared distance went negative"; + EXPECT_LT(std::abs(static_cast(l2) - expected_l2) / expected_l2, 1e-6) + << "scalar L2, dim " << dim; + + unsigned char alignment = 0; + auto dispatched_l2 = L2_UINT8_GetDistFunc(dim, &alignment, nullptr); + const float l2_simd = dispatched_l2(v1.data(), v2.data(), dim); + EXPECT_GT(l2_simd, 0.0f) << "dim " << dim << ": SIMD squared distance went negative"; + EXPECT_LT(std::abs(static_cast(l2_simd) - expected_l2) / expected_l2, 1e-6) + << "dispatched L2, dim " << dim; + + // IP between two all-255 vectors is 255^2 * dim, and the kernel returns 1 - IP. + const double expected_ip = 1.0 - 255.0 * 255.0 * static_cast(dim); + const double ip = static_cast(UINT8_InnerProduct(v1.data(), v1.data(), dim)); + EXPECT_LT(std::abs(ip - expected_ip) / std::abs(expected_ip), 1e-6) + << "scalar IP, dim " << dim; + + auto dispatched_ip = IP_UINT8_GetDistFunc(dim, &alignment, nullptr); + const double ip_simd = static_cast(dispatched_ip(v1.data(), v1.data(), dim)); + EXPECT_LT(std::abs(ip_simd - expected_ip) / std::abs(expected_ip), 1e-6) + << "dispatched IP, dim " << dim; + } +} From bfbcd90666db9b1925112da249def6c9b14556f5 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 12 Aug 2026 15:01:38 +0300 Subject: [PATCH 13/17] style: drop the blank line left by removing the MAX_EXACT_DIM assert make check-format was failing on it. --- src/VecSim/spaces/computer/preprocessors.h | 1 - 1 file changed, 1 deletion(-) diff --git a/src/VecSim/spaces/computer/preprocessors.h b/src/VecSim/spaces/computer/preprocessors.h index 654192c53..5c13681dc 100644 --- a/src/VecSim/spaces/computer/preprocessors.h +++ b/src/VecSim/spaces/computer/preprocessors.h @@ -368,7 +368,6 @@ class QuantPreprocessor : public PreprocessorInterface { } } - // Metadata is a packed run of 4-byte fields: min and delta as FP32, the quantized sums as // uint32, and x_mean_ip (FP32) last for IP with norm. memcpy throughout because the offset // (dim * sizeof(uint8_t)) is not guaranteed to be 4-byte aligned. From 2b66b5f50f8ebd60e22b401469e911ae8b8a2886 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 12 Aug 2026 15:04:37 +0300 Subject: [PATCH 14/17] chore: remove .sites.txt A list of the source paths touched by this change set, committed by accident in 2b84cc2f. Nothing in the repo, build or tests references it. Reported by Cursor Bugbot. --- .sites.txt | 25 ------------------------- 1 file changed, 25 deletions(-) delete mode 100644 .sites.txt diff --git a/.sites.txt b/.sites.txt deleted file mode 100644 index 9bf88f993..000000000 --- a/.sites.txt +++ /dev/null @@ -1,25 +0,0 @@ -src/VecSim/spaces/IP/IP_AVX512F_BW_VL_VNNI_SQ8_SQ8.h -src/VecSim/spaces/IP/IP_NEON_DOTPROD_SQ8_SQ8.h -src/VecSim/spaces/IP/IP_NEON_SQ8_SQ8.h -src/VecSim/spaces/IP/IP_SVE_SQ8_SQ8.h -src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP16.h -src/VecSim/spaces/L2/L2_AVX2_FMA_SQ8_FP32.h -src/VecSim/spaces/L2/L2_AVX2_SQ8_FP16.h -src/VecSim/spaces/L2/L2_AVX2_SQ8_FP32.h -src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_FP32.h -src/VecSim/spaces/L2/L2_AVX512F_BW_VL_VNNI_SQ8_SQ8.h -src/VecSim/spaces/L2/L2_AVX512F_SQ8_FP16.h -src/VecSim/spaces/L2/L2.cpp -src/VecSim/spaces/L2/L2_NEON_DOTPROD_SQ8_SQ8.h -src/VecSim/spaces/L2/L2_NEON_SQ8_FP16.h -src/VecSim/spaces/L2/L2_NEON_SQ8_FP32.h -src/VecSim/spaces/L2/L2_NEON_SQ8_SQ8.h -src/VecSim/spaces/L2/L2_SSE4_SQ8_FP16.h -src/VecSim/spaces/L2/L2_SSE4_SQ8_FP32.h -src/VecSim/spaces/L2/L2_SVE2_SQ8_FP16.h -src/VecSim/spaces/L2/L2_SVE_SQ8_FP16.h -src/VecSim/spaces/L2/L2_SVE_SQ8_FP32.h -src/VecSim/spaces/L2/L2_SVE_SQ8_SQ8.h -tests/unit/test_components.cpp -tests/unit/test_spaces.cpp -tests/utils/tests_utils.h From 9594c8a76947519609524e12090cf736bf209140 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Wed, 12 Aug 2026 15:49:59 +0300 Subject: [PATCH 15/17] test: free the preprocessor in QuantizationHandlesNonRepresentableRange The test freed both blobs but never deleted the QuantPreprocessor it allocated, leaking 88 bytes in 3 allocations. LeakSanitizer failed it in the asan and coverage jobs; a plain debug ctest run does not catch this, which is why it passed locally. Every sibling test in the file already deletes its preprocessor. --- tests/unit/test_components.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/unit/test_components.cpp b/tests/unit/test_components.cpp index 764bc69cd..4ca297fe9 100644 --- a/tests/unit/test_components.cpp +++ b/tests/unit/test_components.cpp @@ -1239,6 +1239,7 @@ TEST(PreprocessorsTest, QuantizationHandlesNonRepresentableRange) { if (query_blob != storage_blob) { allocator->free_allocation(query_blob); } + delete quant_preprocessor; } // Test edge case where all entries are equal From 6c8d6e0dfe48656f7ac95ceb94a06fa9cf7b11b0 Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Thu, 13 Aug 2026 16:48:21 +0300 Subject: [PATCH 16/17] fix(sq8): close the two remaining paths to NaN in quantize Widening the range to double closed the path MOD-17528 described and left two others open, both reported in review and both reaching the same conversion of NaN to an integer. * find_min_max centers in FP32, so input FLT_MAX against mean -FLT_MAX centers to 6.8e38, i.e. inf, before the double range is ever computed. min_val is stored as FP32, so a non-finite endpoint could not be stored under any arithmetic; the endpoints are now clamped to the float range, which is both the storage limit and the guard. Only the endpoints are clamped, so the per-element loop stays FP32 and pays nothing. * delta is stored as FP32, and (float)(diff / 255.0) underflows to zero for any diff below about 1.8e-43. diff itself is not zero there, so testing diff did not catch it. 1/delta was then inf and the minimum element, whose numerator is exactly zero, scaled to 0 * inf = NaN. The guard now tests delta, which is the value that gets stored and inverted, and subsumes the old min == max check. With min_val finite and inv_delta finite and positive, NaN is unreachable in to_byte. The bound there is still moved from std::clamp to std::fmin/std::fmax, which return the non-NaN operand where clamp propagates NaN, so the conversion is defined on its own terms rather than on a proof spanning two functions. This function has now produced three separate conversion-UB defects, which is what makes the backstop worth one instruction. Both cases are pinned by UBSan regressions next to the existing one. Also drops the claim that widening the subtraction removed the whole class, which is what the review was disputing, and names the two places that finish the job instead. Co-Authored-By: Claude Opus 5 (1M context) --- src/VecSim/spaces/computer/preprocessors.h | 58 ++++++++++++-- tests/unit/test_components.cpp | 93 ++++++++++++++++++++++ 2 files changed, 143 insertions(+), 8 deletions(-) diff --git a/src/VecSim/spaces/computer/preprocessors.h b/src/VecSim/spaces/computer/preprocessors.h index 5c13681dc..a2972f7ed 100644 --- a/src/VecSim/spaces/computer/preprocessors.h +++ b/src/VecSim/spaces/computer/preprocessors.h @@ -15,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -285,25 +286,55 @@ class QuantPreprocessor : public PreprocessorInterface { // FP32 arithmetic could leave the representable range for input that is entirely valid. // [-FLT_MAX, +FLT_MAX] made max - min overflow to inf, then delta inf, inv_delta 0, and // finally inf * 0 = NaN, whose conversion to an integer is undefined behaviour. Doubles - // cannot overflow for any pair of finite floats, so the whole class disappears rather than - // being special-cased. delta and inv_delta stay MetadataType, since delta is stored. + // cannot overflow for any pair of finite floats, so the subtraction itself is now safe. + // delta and inv_delta stay MetadataType, since delta is stored. + // + // Widening this subtraction is not on its own sufficient, and the two places that finish + // the job are worth naming here because they are easy to undo: find_min_max clamps its + // endpoints, since centering can produce an inf before this line ever runs, and the delta + // guard below catches the FP32 underflow that happens after it. What makes NaN unreachable + // is those three together, not this line. const double diff = static_cast(max_val) - static_cast(min_val); - const MetadataType delta = - (diff == 0.0) ? MetadataType{1} : static_cast(diff / 255.0); + + // The guard is on delta itself, not on the range it came from, because delta is what gets + // stored and inverted and it has a narrower type than the range. Two ways it goes wrong: + // * (float)(diff / 255.0) underflows to zero for any diff below about 1.8e-43, since the + // smallest positive float is 1.4e-45. A diff of 1e-44 is not zero, so testing diff + // would not catch it. Then 1/delta is inf and the minimum element, whose numerator is + // exactly zero, scales to 0 * inf = NaN. + // * diff is inf or NaN. find_min_max clamps the endpoints so this cannot happen today, + // but delta should not depend on that proof holding elsewhere. + // delta = 1 is the handling min == max already gets: every value collapses to byte 0 and + // the reconstruction is min everywhere, which is the best 8 bits can do for a range this + // narrow, since delta has to be a representable positive float. + MetadataType delta = static_cast(diff / 255.0); + if (!std::isfinite(delta) || delta <= MetadataType{0}) + delta = MetadataType{1}; + const double inv_delta = 1.0 / static_cast(delta); const double min_val_d = static_cast(min_val); // Scale one input value to its quantized byte. Kept in double for the reason above, and - // clamped so that a value outside [0, 255] can never reach the conversion: the arithmetic + // bounded so that a value outside [0, 255] can never reach the conversion: the arithmetic // should not produce one, and if it ever does the result must be a wrong byte rather than // undefined behaviour. + // + // With min_val finite and inv_delta finite and positive, both guaranteed above, NaN is + // unreachable here: x is finite or +/-inf but never NaN, so the numerator is finite or + // +/-inf, and multiplying that by a finite non-zero can give inf but not NaN. An inf + // element (only possible from FP32 centering) therefore lands on 255, which is where the + // largest value belongs. The bound is still written NaN-tolerantly so that the conversion + // is defined on its own terms rather than on that argument holding. const auto to_byte = [inv_delta, min_val_d](float x) { const double scaled = (static_cast(x) - min_val_d) * inv_delta; - // Clamped first, so the value is non-negative and adding 0.5 before truncating is + // fmin/fmax are IEEE minNum/maxNum: they return the non-NaN operand, so a NaN would + // land on 0 rather than reaching the conversion. std::clamp is written as comparisons, + // every comparison against NaN is false, and it therefore propagates NaN. + // Bounded first, so the value is non-negative and adding 0.5 before truncating is // std::round without the library call. This translation unit is compiled at the // x86-64 baseline, where std::round is an out-of-line call to round() rather than a // roundsd instruction, and this runs once per element. - return static_cast(std::clamp(scaled, 0.0, 255.0) + 0.5); + return static_cast(std::fmin(std::fmax(scaled, 0.0), 255.0) + 0.5); }; // Sum the *quantized* bytes, not the input floats. A stored blob only ever describes its @@ -631,7 +662,18 @@ class QuantPreprocessor : public PreprocessorInterface { if constexpr (Metric == VecSimMetric_IP) x_mean_ip += input_value * mean[i]; } - return {min_val, max_val}; + // Centering is an FP32 subtraction of two finite floats, which can still leave the + // FP32 range: input FLT_MAX against mean -FLT_MAX centers to 6.8e38, i.e. +inf. The + // range is what delta is derived from, and min_val is stored as FP32, so a non-finite + // endpoint could not be stored under any arithmetic. Clamping it here is therefore + // both the storage limit and the guard: it keeps inf out of delta, where 1/inf = 0 + // would turn the per-element scaling into inf * 0 = NaN. Only the endpoints are + // clamped, so the per-element loop stays FP32 and pays nothing. min_val and max_val + // are finite or +/-inf but never NaN here, since finite minus finite cannot be NaN, + // so std::clamp (which would propagate a NaN) is safe on them. + constexpr float representable = std::numeric_limits::max(); + return {std::clamp(min_val, -representable, representable), + std::clamp(max_val, -representable, representable)}; } } diff --git a/tests/unit/test_components.cpp b/tests/unit/test_components.cpp index 4ca297fe9..983e826b3 100644 --- a/tests/unit/test_components.cpp +++ b/tests/unit/test_components.cpp @@ -1242,6 +1242,99 @@ TEST(PreprocessorsTest, QuantizationHandlesNonRepresentableRange) { delete quant_preprocessor; } +// The other direction of the same defect: a range too narrow to derive a usable delta from, rather +// than too wide. diff is not zero, so the min == max guard does not fire, but (float)(diff / 255) +// underflows to zero because the smallest positive float is 1.4e-45. 1 / 0 is inf, and the minimum +// element, whose numerator is exactly zero, then scales to 0 * inf = NaN, whose conversion to an +// integer is undefined behaviour. Runs under UBSan in CI. +TEST(PreprocessorsTest, QuantizationHandlesSubnormalRange) { + std::shared_ptr allocator = VecSimAllocator::newVecsimAllocator(); + constexpr size_t dim = 4; + constexpr unsigned char alignment = 0; + // 1e-44 is subnormal, so diff / 255 has no representable non-zero FP32 value. + float original_blob[dim] = {0.0f, 1e-44f, 0.0f, 1e-44f}; + + auto quant_preprocessor = + new (allocator) QuantPreprocessor(allocator, dim); + + void *storage_blob = nullptr; + void *query_blob = nullptr; + size_t storage_blob_size = dim * sizeof(float); + size_t query_blob_size = dim * sizeof(float); + quant_preprocessor->preprocess(original_blob, storage_blob, query_blob, storage_blob_size, + query_blob_size, alignment, alignment); + ASSERT_NE(storage_blob, nullptr); + + // delta must be a usable positive float, not zero, because every kernel divides by it or + // multiplies the sums by it. + const auto *quantized = static_cast(storage_blob); + const auto *metadata = quantized + dim; + const float delta = load_unaligned(metadata + sq8::DELTA * sizeof(float)); + EXPECT_TRUE(std::isfinite(delta)); + EXPECT_GT(delta, 0.0f); + EXPECT_TRUE(std::isfinite(load_unaligned(metadata + sq8::MIN_VAL * sizeof(float)))); + + // The range is far below one quantization step, so collapsing to a single byte is the correct + // outcome. Which byte is not the point; not being undefined behaviour is. + for (size_t i = 0; i < dim; ++i) { + EXPECT_EQ(quantized[i], quantized[0]) << "index " << i; + } + + allocator->free_allocation(storage_blob); + if (query_blob != storage_blob) { + allocator->free_allocation(query_blob); + } + delete quant_preprocessor; +} + +// The same defect reached through centering instead of through the raw range. find_min_max and +// transformed_value compute input[i] - mean[i] in FP32, so finite input against a finite mean can +// center to +/-inf: FLT_MAX against -FLT_MAX is 6.8e38. That inf used to reach delta, where +// 1 / inf = 0 turned the per-element scaling into inf * 0 = NaN. Widening the range to double does +// not help, because the value is already lost upstream in FP32. Runs under UBSan in CI. +TEST(PreprocessorsTest, QuantizationHandlesNonRepresentableCenteredRange) { + std::shared_ptr allocator = VecSimAllocator::newVecsimAllocator(); + constexpr size_t dim = 4; + constexpr unsigned char alignment = 0; + constexpr float huge = std::numeric_limits::max(); + + float original_blob[dim] = {huge, 0.0f, -huge, 1.0f}; + vecsim_stl::vector mean_vec(allocator); + // Centers element 0 to +inf and element 2 to -inf, from entirely finite operands. + mean_vec.push_back(-huge); + mean_vec.push_back(0.0f); + mean_vec.push_back(huge); + mean_vec.push_back(0.0f); + + auto quant_preprocessor = + new (allocator) QuantPreprocessor(allocator, dim, mean_vec); + + void *storage_blob = nullptr; + size_t storage_blob_size = dim * sizeof(float); + quant_preprocessor->preprocessForStorage(original_blob, storage_blob, storage_blob_size, + alignment); + ASSERT_NE(storage_blob, nullptr); + + // Metadata must be finite, since the kernels derive every term from it, and delta must be + // positive so the reconstruction is not degenerate. + const auto *quantized = static_cast(storage_blob); + const auto *metadata = quantized + dim; + const float min_val = load_unaligned(metadata + sq8::MIN_VAL * sizeof(float)); + const float delta = load_unaligned(metadata + sq8::DELTA * sizeof(float)); + EXPECT_TRUE(std::isfinite(min_val)); + EXPECT_TRUE(std::isfinite(delta)); + EXPECT_GT(delta, 0.0f); + + // The centered extremes land on the ends of the byte range, which is where they belong. The + // two middle values collapse together, as in the non-centered case: one byte cannot resolve a + // range this wide, and losing resolution is the correct outcome where UB was not. + EXPECT_EQ(quantized[0], 255); + EXPECT_EQ(quantized[2], 0); + + allocator->free_allocation(storage_blob); + delete quant_preprocessor; +} + // Test edge case where all entries are equal TEST(PreprocessorsTest, QuantizationTestAllEntriesEqual) { std::shared_ptr allocator = VecSimAllocator::newVecsimAllocator(); From 729f6214043dd527a335bb57d24ab5ff64ab39fd Mon Sep 17 00:00:00 2001 From: Dor Forer Date: Thu, 13 Aug 2026 17:23:21 +0300 Subject: [PATCH 17/17] docs(sq8): state the quantize numerics once instead of per line The reasoning had accumulated as a block above nearly every statement, which buried the code it was explaining. It is now one description above quantize(), covering why the scaling is in double and which three guards keep NaN out of the byte conversion. Two inline comments are kept: the +0.5 trick, which reads like an oversight and invites being changed back to std::round, and one naming what transformed_value does. Also drops the pre-existing "We know (input - min) => 0", which stopped being true once the endpoints were clamped and would have someone read the fmax bound in to_byte as dead code. Comments only. With comments and blank lines stripped, the file is identical to the previous commit. Co-Authored-By: Claude Opus 5 (1M context) --- src/VecSim/spaces/computer/preprocessors.h | 93 ++++++---------------- 1 file changed, 25 insertions(+), 68 deletions(-) diff --git a/src/VecSim/spaces/computer/preprocessors.h b/src/VecSim/spaces/computer/preprocessors.h index a2972f7ed..532eff487 100644 --- a/src/VecSim/spaces/computer/preprocessors.h +++ b/src/VecSim/spaces/computer/preprocessors.h @@ -274,39 +274,29 @@ class QuantPreprocessor : public PreprocessorInterface { static_assert(!WithNorm || Metric != VecSimMetric_Cosine, "WithNorm does not support Cosine metric."); - // Helper function to perform quantization. This function is used by the storage preprocessing - // methods. + // Quantizes one input vector into dim bytes followed by its metadata run. Used by the storage + // preprocessing methods. + // + // Scaling is in double: [-FLT_MAX, +FLT_MAX] is valid input and max - min overflows FP32. + // Three things together keep NaN out of the conversion to a byte, and none is sufficient + // alone. find_min_max clamps its endpoints, because FP32 centering can already hold an inf and + // because min is stored as FP32 and could not carry one anyway. delta is guarded rather than + // diff, since (float)(diff / 255) underflows to zero below diff ~1.8e-43 while diff itself is + // still nonzero, leaving 1/delta = inf and scaling the minimum element, whose numerator is + // exactly zero, to 0 * inf. And to_byte bounds with fmin/fmax, which return the non-NaN + // operand where std::clamp propagates it. The first two make NaN unreachable on their own, so + // the third is a backstop, kept because this function has produced three separate + // conversion-UB defects. delta = 1 is the degenerate handling min == max already got. + // + // The sums are exact uint32 over the quantized bytes rather than the input floats; see the + // class documentation above for why, and sq8.h for the derivations that consume them. void quantize(const DataType *input, OUTPUT_TYPE *quantized) const { assert(input && quantized); float x_mean_ip = 0.0f; // only used for WithNorm auto [min_val, max_val] = find_min_max(input, x_mean_ip); - // The range is computed in double, and so is the per-element scaling below, because the - // FP32 arithmetic could leave the representable range for input that is entirely valid. - // [-FLT_MAX, +FLT_MAX] made max - min overflow to inf, then delta inf, inv_delta 0, and - // finally inf * 0 = NaN, whose conversion to an integer is undefined behaviour. Doubles - // cannot overflow for any pair of finite floats, so the subtraction itself is now safe. - // delta and inv_delta stay MetadataType, since delta is stored. - // - // Widening this subtraction is not on its own sufficient, and the two places that finish - // the job are worth naming here because they are easy to undo: find_min_max clamps its - // endpoints, since centering can produce an inf before this line ever runs, and the delta - // guard below catches the FP32 underflow that happens after it. What makes NaN unreachable - // is those three together, not this line. const double diff = static_cast(max_val) - static_cast(min_val); - - // The guard is on delta itself, not on the range it came from, because delta is what gets - // stored and inverted and it has a narrower type than the range. Two ways it goes wrong: - // * (float)(diff / 255.0) underflows to zero for any diff below about 1.8e-43, since the - // smallest positive float is 1.4e-45. A diff of 1e-44 is not zero, so testing diff - // would not catch it. Then 1/delta is inf and the minimum element, whose numerator is - // exactly zero, scales to 0 * inf = NaN. - // * diff is inf or NaN. find_min_max clamps the endpoints so this cannot happen today, - // but delta should not depend on that proof holding elsewhere. - // delta = 1 is the handling min == max already gets: every value collapses to byte 0 and - // the reconstruction is min everywhere, which is the best 8 bits can do for a range this - // narrow, since delta has to be a representable positive float. MetadataType delta = static_cast(diff / 255.0); if (!std::isfinite(delta) || delta <= MetadataType{0}) delta = MetadataType{1}; @@ -314,53 +304,26 @@ class QuantPreprocessor : public PreprocessorInterface { const double inv_delta = 1.0 / static_cast(delta); const double min_val_d = static_cast(min_val); - // Scale one input value to its quantized byte. Kept in double for the reason above, and - // bounded so that a value outside [0, 255] can never reach the conversion: the arithmetic - // should not produce one, and if it ever does the result must be a wrong byte rather than - // undefined behaviour. - // - // With min_val finite and inv_delta finite and positive, both guaranteed above, NaN is - // unreachable here: x is finite or +/-inf but never NaN, so the numerator is finite or - // +/-inf, and multiplying that by a finite non-zero can give inf but not NaN. An inf - // element (only possible from FP32 centering) therefore lands on 255, which is where the - // largest value belongs. The bound is still written NaN-tolerantly so that the conversion - // is defined on its own terms rather than on that argument holding. const auto to_byte = [inv_delta, min_val_d](float x) { const double scaled = (static_cast(x) - min_val_d) * inv_delta; - // fmin/fmax are IEEE minNum/maxNum: they return the non-NaN operand, so a NaN would - // land on 0 rather than reaching the conversion. std::clamp is written as comparisons, - // every comparison against NaN is false, and it therefore propagates NaN. - // Bounded first, so the value is non-negative and adding 0.5 before truncating is - // std::round without the library call. This translation unit is compiled at the - // x86-64 baseline, where std::round is an out-of-line call to round() rather than a - // roundsd instruction, and this runs once per element. + // Bounded first, so +0.5 before truncating is std::round without the library call: + // this TU is compiled at the x86-64 baseline, where round() is an out-of-line call. return static_cast(std::fmin(std::fmax(scaled, 0.0), 255.0) + 0.5); }; - // Sum the *quantized* bytes, not the input floats. A stored blob only ever describes its - // reconstruction min + delta * a[i], so metadata derived from it has to describe that same - // vector; summing the input instead mixes two different vectors into one distance formula. - // Integer accumulation is also exact, which the L2 identity depends on (see sq8.h). - // 4 independent accumulators each, mirroring the unrolled loop below. + // 4 independent accumulators each, mirroring the unrolled loop below. q* is L2 only. uint32_t s0{}, s1{}, s2{}, s3{}; - - // Sum of squares, only used for L2. uint32_t q0{}, q1{}, q2{}, q3{}; size_t i = 0; - // round dim down to the nearest multiple of 4 size_t dim_round_down = this->dim & ~size_t(3); - // Quantize the values for (; i < dim_round_down; i += 4) { - // Load once (widened to FP32 if DataType is FP16). + // transformed_value widens to FP32 if DataType is FP16, and centers if WithNorm. const float x0 = transformed_value(input, i); const float x1 = transformed_value(input, i + 1); const float x2 = transformed_value(input, i + 2); const float x3 = transformed_value(input, i + 3); - // We know (input - min) => 0 - // If min == max, all values are the same and should be quantized to 0. - // reconstruction will yield the same original value for all vectors. const uint32_t a0 = to_byte(x0); const uint32_t a1 = to_byte(x1); const uint32_t a2 = to_byte(x2); @@ -370,13 +333,11 @@ class QuantPreprocessor : public PreprocessorInterface { quantized[i + 2] = static_cast(a2); quantized[i + 3] = static_cast(a3); - // Accumulate sum for all metrics s0 += a0; s1 += a1; s2 += a2; s3 += a3; - // Accumulate sum of squares only for L2 metric if constexpr (Metric == VecSimMetric_L2) { q0 += a0 * a0; q1 += a1 * a1; @@ -662,15 +623,11 @@ class QuantPreprocessor : public PreprocessorInterface { if constexpr (Metric == VecSimMetric_IP) x_mean_ip += input_value * mean[i]; } - // Centering is an FP32 subtraction of two finite floats, which can still leave the - // FP32 range: input FLT_MAX against mean -FLT_MAX centers to 6.8e38, i.e. +inf. The - // range is what delta is derived from, and min_val is stored as FP32, so a non-finite - // endpoint could not be stored under any arithmetic. Clamping it here is therefore - // both the storage limit and the guard: it keeps inf out of delta, where 1/inf = 0 - // would turn the per-element scaling into inf * 0 = NaN. Only the endpoints are - // clamped, so the per-element loop stays FP32 and pays nothing. min_val and max_val - // are finite or +/-inf but never NaN here, since finite minus finite cannot be NaN, - // so std::clamp (which would propagate a NaN) is safe on them. + // Centering is an FP32 subtraction of finite floats, so it can reach +/-inf: FLT_MAX + // against mean -FLT_MAX gives 6.8e38. Clamping the endpoints is both the storage limit + // (min is stored as FP32) and the guard quantize() relies on; only the endpoints are + // clamped, so the loop above stays FP32. std::clamp is safe on these two, which are + // finite or +/-inf but never NaN, since finite minus finite cannot be NaN. constexpr float representable = std::numeric_limits::max(); return {std::clamp(min_val, -representable, representable), std::clamp(max_val, -representable, representable)};