Added support for binary (1bit) vectors and SIMD-Optimized Hamming distance

Author: marcobambini

src/distance-avx2.c

@@ -949,6 +949,47 @@ float int8_distance_cosine_avx2 (const void *a, const void *b, int n) {
     return 1.0f - cosine_similarity;
 }
 
+// MARK: - BIT -
+
+// lookup table for popcount of 4-bit values
+static const __m256i popcount_lut = _mm256_setr_epi8(0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4);
+
+static inline __m256i popcount_avx2(__m256i v) {
+    __m256i low_mask = _mm256_set1_epi8(0x0f);
+    __m256i lo = _mm256_and_si256(v, low_mask);
+    __m256i hi = _mm256_and_si256(_mm256_srli_epi16(v, 4), low_mask);
+    __m256i cnt_lo = _mm256_shuffle_epi8(popcount_lut, lo);
+    __m256i cnt_hi = _mm256_shuffle_epi8(popcount_lut, hi);
+    return _mm256_add_epi8(cnt_lo, cnt_hi);
+}
+
+float bit1_distance_hamming_avx2 (const void *v1, const void *v2, int n) {
+    const uint8_t *a = (const uint8_t *)v1;
+    const uint8_t *b = (const uint8_t *)v2;
+    __m256i acc = _mm256_setzero_si256();
+    int i = 0;
+    
+    // Process 32 bytes at a time
+    for (; i + 32 <= n; i += 32) {
+        __m256i va = _mm256_loadu_si256((const __m256i *)(a + i));
+        __m256i vb = _mm256_loadu_si256((const __m256i *)(b + i));
+        __m256i xored = _mm256_xor_si256(va, vb);
+        __m256i popcnt = popcount_avx2(xored);
+        acc = _mm256_add_epi64(acc, _mm256_sad_epu8(popcnt, _mm256_setzero_si256()));
+    }
+    
+    // Horizontal sum
+    __m128i sum128 = _mm_add_epi64(_mm256_extracti128_si256(acc, 0), _mm256_extracti128_si256(acc, 1));
+    int distance = _mm_extract_epi64(sum128, 0) + _mm_extract_epi64(sum128, 1);
+    
+    // Handle remainder with scalar
+    for (; i < n; i++) {
+        distance += __builtin_popcount(a[i] ^ b[i]);
+    }
+    
+    return (float)distance;
+}
+
 #endif
 
 // MARK: -
@@ -985,6 +1026,8 @@ void init_distance_functions_avx2 (void) {
     dispatch_distance_table[VECTOR_DISTANCE_L1][VECTOR_TYPE_U8] = uint8_distance_l1_avx2;
     dispatch_distance_table[VECTOR_DISTANCE_L1][VECTOR_TYPE_I8] = int8_distance_l1_avx2;
 
+    dispatch_distance_table[VECTOR_DISTANCE_HAMMING][VECTOR_TYPE_BIT] = bit1_distance_hamming_avx2;
+    
     distance_backend_name = "AVX2";
 #endif
 }

src/distance-avx512.c

@@ -48,7 +48,7 @@ static inline bool block_has_l2_inf_mismatch_16(const uint16_t* a, const uint16_
     return false;
 }
 
-/* 16×bf16 -> 16×f32: widen to u32, shift <<16, reinterpret as f32 */
+/* 16�bf16 -> 16�f32: widen to u32, shift <<16, reinterpret as f32 */
 static inline __m512 bf16x16_to_f32x16_loadu(const uint16_t* p) {
     // Load 16x u16 (256 bits)
     __m256i v16 = _mm256_loadu_si256((const __m256i*)p);
@@ -846,6 +846,72 @@ float int8_distance_cosine_avx512(const void* a, const void* b, int n) {
     return 1.0f - cosine_similarity;
 }
 
+// MARK: - BIT -
+
+// AVX-512 popcount using lookup table (works on all AVX-512 CPUs)
+static inline __m512i popcount_avx512(__m512i v) {
+    // Lookup table for popcount of 4-bit values
+    const __m512i popcount_lut = _mm512_set_epi8(
+        4, 3, 3, 2, 3, 2, 2, 1, 3, 2, 2, 1, 2, 1, 1, 0,
+        4, 3, 3, 2, 3, 2, 2, 1, 3, 2, 2, 1, 2, 1, 1, 0,
+        4, 3, 3, 2, 3, 2, 2, 1, 3, 2, 2, 1, 2, 1, 1, 0,
+        4, 3, 3, 2, 3, 2, 2, 1, 3, 2, 2, 1, 2, 1, 1, 0
+    );
+    const __m512i low_mask = _mm512_set1_epi8(0x0f);
+
+    __m512i lo = _mm512_and_si512(v, low_mask);
+    __m512i hi = _mm512_and_si512(_mm512_srli_epi16(v, 4), low_mask);
+    __m512i cnt_lo = _mm512_shuffle_epi8(popcount_lut, lo);
+    __m512i cnt_hi = _mm512_shuffle_epi8(popcount_lut, hi);
+    return _mm512_add_epi8(cnt_lo, cnt_hi);
+}
+
+// Hamming distance for 1-bit packed binary vectors
+// n = number of dimensions (bits), not bytes
+static float bit1_distance_hamming_avx512(const void *v1, const void *v2, int n) {
+    const uint8_t *a = (const uint8_t *)v1;
+    const uint8_t *b = (const uint8_t *)v2;
+    int num_bytes = (n + 7) / 8;
+
+    __m512i acc = _mm512_setzero_si512();
+    int i = 0;
+
+    // Process 64 bytes at a time
+    for (; i + 64 <= num_bytes; i += 64) {
+        __m512i va = _mm512_loadu_si512((const __m512i *)(a + i));
+        __m512i vb = _mm512_loadu_si512((const __m512i *)(b + i));
+        __m512i xored = _mm512_xor_si512(va, vb);
+
+#if defined(__AVX512VPOPCNTDQ__)
+        // Native popcount (Ice Lake+)
+        __m512i popcnt = _mm512_popcnt_epi64(xored);
+        acc = _mm512_add_epi64(acc, popcnt);
+#else
+        // Lookup table popcount (Skylake-X compatible)
+        __m512i popcnt = popcount_avx512(xored);
+        // Sum bytes to 64-bit using SAD against zero
+        acc = _mm512_add_epi64(acc, _mm512_sad_epu8(popcnt, _mm512_setzero_si512()));
+#endif
+    }
+
+    // Horizontal sum
+    uint64_t distance = _mm512_reduce_add_epi64(acc);
+
+    // Handle remaining bytes with scalar code
+    for (; i < num_bytes; i++) {
+#if defined(__GNUC__) || defined(__clang__)
+        distance += __builtin_popcount(a[i] ^ b[i]);
+#else
+        uint8_t x = a[i] ^ b[i];
+        x = x - ((x >> 1) & 0x55);
+        x = (x & 0x33) + ((x >> 2) & 0x33);
+        distance += (x + (x >> 4)) & 0x0f;
+#endif
+    }
+
+    return (float)distance;
+}
+
 #endif
 
 // MARK: -
@@ -882,6 +948,8 @@ void init_distance_functions_avx512(void) {
     dispatch_distance_table[VECTOR_DISTANCE_L1][VECTOR_TYPE_U8] = uint8_distance_l1_avx512;
     dispatch_distance_table[VECTOR_DISTANCE_L1][VECTOR_TYPE_I8] = int8_distance_l1_avx512;
 
+    dispatch_distance_table[VECTOR_DISTANCE_HAMMING][VECTOR_TYPE_BIT] = bit1_distance_hamming_avx512;
+
     distance_backend_name = "AVX512";
 #endif
-}
+}

src/distance-avx512.h

@@ -10,6 +10,6 @@
 
 #include <stdio.h>
 
-void init_distance_functions_avx512(void);
+void init_distance_functions_avx512 (void);
 
 #endif

src/distance-cpu.c

@@ -693,6 +693,42 @@ float int8_distance_l1_cpu (const void *v1, const void *v2, int n) {
     return sum;
 }
 
+// MARK: - BIT -
+
+static inline int popcount64(uint64_t x) {
+    #if defined(__GNUC__) || defined(__clang__)
+    return __builtin_popcountll(x);
+    #else
+    // fallback: bit manipulation
+    x = x - ((x >> 1) & 0x5555555555555555ULL);
+    x = (x & 0x3333333333333333ULL) + ((x >> 2) & 0x3333333333333333ULL);
+    x = (x + (x >> 4)) & 0x0f0f0f0f0f0f0f0fULL;
+    return (x * 0x0101010101010101ULL) >> 56;
+    #endif
+}
+
+float bit1_distance_hamming_cpu (const void *v1, const void *v2, int n) {
+    const uint8_t *a = (const uint8_t *)v1;
+    const uint8_t *b = (const uint8_t *)v2;
+    
+    int distance = 0;
+    int i = 0;
+    
+    // process 8 bytes at a time
+    for (; i + 8 <= n; i += 8) {
+        uint64_t xa = *(const uint64_t *)(a + i);
+        uint64_t xb = *(const uint64_t *)(b + i);
+        distance += popcount64(xa ^ xb);
+    }
+    
+    // handle remainder
+    for (; i < n; i++) {
+        distance += popcount64(a[i] ^ b[i]);
+    }
+    
+    return (float)distance;
+}
+
 // MARK: - ENTRYPOINT -
 
 #if defined(__x86_64__) || defined(_M_X64) || defined(__i386__) || defined(_M_IX86)
@@ -845,6 +881,9 @@ void init_cpu_functions (void) {
                 [VECTOR_TYPE_BF16] = bfloat16_distance_l1_cpu,
                 [VECTOR_TYPE_U8]  = uint8_distance_l1_cpu,
                 [VECTOR_TYPE_I8]  = int8_distance_l1_cpu,
+            },
+            [VECTOR_DISTANCE_HAMMING] = {
+                [VECTOR_TYPE_BIT] = bit1_distance_hamming_cpu
             }
     };
 

src/distance-cpu.h

@@ -38,14 +38,16 @@ typedef enum {
     VECTOR_TYPE_F16,
     VECTOR_TYPE_BF16,
     VECTOR_TYPE_U8,
-    VECTOR_TYPE_I8
+    VECTOR_TYPE_I8,
+    VECTOR_TYPE_BIT
 } vector_type;
-#define VECTOR_TYPE_MAX         6
+#define VECTOR_TYPE_MAX         7
 
 typedef enum {
     VECTOR_QUANT_AUTO = 0,
     VECTOR_QUANT_U8BIT = 1,
-    VECTOR_QUANT_S8BIT = 2
+    VECTOR_QUANT_S8BIT = 2,
+    VECTOR_QUANT_1BIT = 3
 } vector_qtype;
 
 typedef enum {
@@ -54,8 +56,9 @@ typedef enum {
     VECTOR_DISTANCE_COSINE,
     VECTOR_DISTANCE_DOT,
     VECTOR_DISTANCE_L1,
+    VECTOR_DISTANCE_HAMMING
 } vector_distance;
-#define VECTOR_DISTANCE_MAX     6
+#define VECTOR_DISTANCE_MAX     7
 
 typedef float (*distance_function_t)(const void *v1, const void *v2, int n);
 

src/distance-neon.c

@@ -1230,6 +1230,38 @@ float int8_distance_l1_neon(const void *v1, const void *v2, int n) {
 
     return (float)final;
 }
+
+// MARK: - BIT -
+
+float bit1_distance_hamming_neon (const void *v1, const void *v2, int n) {
+    const uint8_t *a = (const uint8_t *)v1;
+    const uint8_t *b = (const uint8_t *)v2;
+    uint64x2_t acc = vdupq_n_u64(0);
+    int i = 0;
+    
+    // Process 16 bytes at a time
+    for (; i + 16 <= n; i += 16) {
+        uint8x16_t va = vld1q_u8(a + i);
+        uint8x16_t vb = vld1q_u8(b + i);
+        uint8x16_t xored = veorq_u8(va, vb);
+        
+        // vcntq_u8: popcount per byte
+        uint8x16_t popcnt = vcntq_u8(xored);
+        
+        // Sum bytes to 64-bit accumulators
+        acc = vpadalq_u32(acc, vpaddlq_u16(vpaddlq_u8(popcnt)));
+    }
+    
+    int distance = (int)(vgetq_lane_u64(acc, 0) + vgetq_lane_u64(acc, 1));
+    
+    // Handle remainder
+    for (; i < n; i++) {
+        distance += __builtin_popcount(a[i] ^ b[i]);
+    }
+    
+    return (float)distance;
+}
+
 #endif
 
 // MARK: -
@@ -1266,6 +1298,8 @@ void init_distance_functions_neon (void) {
     dispatch_distance_table[VECTOR_DISTANCE_L1][VECTOR_TYPE_U8] = uint8_distance_l1_neon;
     dispatch_distance_table[VECTOR_DISTANCE_L1][VECTOR_TYPE_I8] = int8_distance_l1_neon;
 
+    dispatch_distance_table[VECTOR_DISTANCE_HAMMING][VECTOR_TYPE_BIT] = bit1_distance_hamming_neon;
+    
     distance_backend_name = "NEON";
 #endif
 }

src/distance-sse2.c

@@ -1006,6 +1006,67 @@ float int8_distance_cosine_sse2 (const void *v1, const void *v2, int n) {
     return 1.0f - cosine_sim;
 }
 
+// MARK: - BIT -
+
+static inline __m128i popcount_sse2 (__m128i v) {
+    // Classic parallel bit count algorithm vectorized for SSE2
+    
+    const __m128i mask1 = _mm_set1_epi8(0x55);  // 01010101
+    const __m128i mask2 = _mm_set1_epi8(0x33);  // 00110011
+    const __m128i mask4 = _mm_set1_epi8(0x0f);  // 00001111
+    
+    // x = x - ((x >> 1) & 0x55555555)
+    __m128i t = _mm_and_si128(_mm_srli_epi16(v, 1), mask1);
+    v = _mm_sub_epi8(v, t);
+    
+    // x = (x & 0x33333333) + ((x >> 2) & 0x33333333)
+    t = _mm_and_si128(_mm_srli_epi16(v, 2), mask2);
+    v = _mm_add_epi8(_mm_and_si128(v, mask2), t);
+    
+    // x = (x + (x >> 4)) & 0x0f0f0f0f
+    t = _mm_srli_epi16(v, 4);
+    v = _mm_and_si128(_mm_add_epi8(v, t), mask4);
+    
+    // Now each byte contains popcount for that byte (0-8)
+    return v;
+}
+
+float bit1_distance_hamming_sse2 (const void *v1, const void *v2, int n) {
+    const uint8_t *a = (const uint8_t *)v1;
+    const uint8_t *b = (const uint8_t *)v2;
+    __m128i acc = _mm_setzero_si128();
+    int i = 0;
+    
+    // Process 16 bytes at a time
+    for (; i + 16 <= n; i += 16) {
+        __m128i va = _mm_loadu_si128((const __m128i *)(a + i));
+        __m128i vb = _mm_loadu_si128((const __m128i *)(b + i));
+        __m128i xored = _mm_xor_si128(va, vb);
+        __m128i popcnt = popcount_sse2(xored);
+        
+        // Sum bytes using SAD (sum of absolute differences against zero)
+        // This sums all 16 bytes into two 64-bit values
+        acc = _mm_add_epi64(acc, _mm_sad_epu8(popcnt, _mm_setzero_si128()));
+    }
+    
+    // Horizontal sum of the two 64-bit accumulators
+    int distance = _mm_cvtsi128_si64(acc) + _mm_cvtsi128_si64(_mm_srli_si128(acc, 8));
+    
+    // Handle remainder with scalar code
+    for (; i < n; i++) {
+        #if defined(__GNUC__) || defined(__clang__)
+        distance += __builtin_popcount(a[i] ^ b[i]);
+        #else
+        uint8_t x = a[i] ^ b[i];
+        x = x - ((x >> 1) & 0x55);
+        x = (x & 0x33) + ((x >> 2) & 0x33);
+        distance += (x + (x >> 4)) & 0x0f;
+        #endif
+    }
+    
+    return (float)distance;
+}
+
 #endif
 
 // MARK: -
@@ -1042,6 +1103,8 @@ void init_distance_functions_sse2 (void) {
     dispatch_distance_table[VECTOR_DISTANCE_L1][VECTOR_TYPE_U8] = uint8_distance_l1_sse2;
     dispatch_distance_table[VECTOR_DISTANCE_L1][VECTOR_TYPE_I8] = int8_distance_l1_sse2;
 
+    dispatch_distance_table[VECTOR_DISTANCE_HAMMING][VECTOR_TYPE_BIT] = bit1_distance_hamming_sse2;
+    
     distance_backend_name = "SSE2";
 #endif
 }