Add: u8 kernels & FMA

Author: ashvardanian

README.md

@@ -91,15 +91,16 @@ You can learn more about the technical implementation details in the following b
 ## Benchmarks
 
 For reference, we use 1536-dimensional vectors, like the embeddings produced by the OpenAI Ada API.
-Comparing the serial code throughput produced by GCC 12 to hand-optimized kernels in SimSIMD, we see the following single-core improvements:
+Comparing the serial code throughput produced by GCC 12 to hand-optimized kernels in SimSIMD, we see the following single-core improvements for the two most common vector-vector similarity metrics - the Cosine similarity and the Euclidean distance:
 
-| Type   |                       Apple M2 Pro |                          AMD Genoa |                     AWS Graviton 4 |
-| :----- | ---------------------------------: | ---------------------------------: | ---------------------------------: |
-| `f64`  | 18.5 → 28.8 GB/s <br/>      + 56 % | 21.9 → 41.4 GB/s <br/>      + 89 % | 20.7 → 41.3 GB/s <br/>      + 99 % |
-| `f32`  |  9.2 → 29.6 GB/s <br/>     + 221 % | 10.9 → 95.8 GB/s <br/>     + 779 % |  4.9 → 41.9 GB/s <br/>     + 755 % |
-| `f16`  |  4.6 → 14.6 GB/s <br/>     + 217 % | 3.1 → 108.4 GB/s <br/>   + 3,397 % |  5.4 → 39.3 GB/s <br/>     + 627 % |
-| `bf16` |  4.6 → 26.3 GB/s <br/>     + 472 % | 0.8 → 59.5 GB/s <br/>     +7,437 % |  2.5 → 29.9 GB/s <br/>   + 1,096 % |
-| `i8`   | 25.8 → 47.1 GB/s <br/>      + 83 % | 33.1 → 65.3 GB/s <br/>      + 97 % | 35.2 → 43.5 GB/s <br/>      + 24 % |
+| Type   |                  Apple M2 Pro |            Intel Sapphire Rapids |                  AWS Graviton 4 |
+| :----- | ----------------------------: | -------------------------------: | ------------------------------: |
+| `f64`  | 18.5 → 28.8 GB/s <br/> + 56 % |    21.9 → 41.4 GB/s <br/> + 89 % |   20.7 → 41.3 GB/s <br/> + 99 % |
+| `f32`  | 9.2 → 29.6 GB/s <br/> + 221 % |   10.9 → 95.8 GB/s <br/> + 779 % |   4.9 → 41.9 GB/s <br/> + 755 % |
+| `f16`  | 4.6 → 14.6 GB/s <br/> + 217 % | 3.1 → 108.4 GB/s <br/> + 3,397 % |   5.4 → 39.3 GB/s <br/> + 627 % |
+| `bf16` | 4.6 → 26.3 GB/s <br/> + 472 % |   0.8 → 59.5 GB/s <br/> +7,437 % | 2.5 → 29.9 GB/s <br/> + 1,096 % |
+| `i8`   | 25.8 → 47.1 GB/s <br/> + 83 % |    33.1 → 65.3 GB/s <br/> + 97 % |   35.2 → 43.5 GB/s <br/> + 24 % |
+| `u8`   |                               |   32.5 → 66.5 GB/s <br/> + 105 % |                                 |
 
 Similar speedups are often observed even when compared to BLAS and LAPACK libraries underlying most numerical computing libraries, including NumPy and SciPy in Python.
 Broader benchmarking results:
@@ -112,7 +113,7 @@ Broader benchmarking results:
 
 The package is intended to replace the usage of `numpy.inner`, `numpy.dot`, and `scipy.spatial.distance`.
 Aside from drastic performance improvements, SimSIMD significantly improves accuracy in mixed precision setups.
-NumPy and SciPy, processing `i8` or `f16` vectors, will use the same types for accumulators, while SimSIMD can combine `i8` enumeration, `i16` multiplication, and `i32` accumulation to avoid overflows entirely.
+NumPy and SciPy, processing `i8`, `u8` or `f16` vectors, will use the same types for accumulators, while SimSIMD can combine `i8` enumeration, `i16` multiplication, and `i32` accumulation to avoid overflows entirely.
 The same applies to processing `f16` and `bf16` values with `f32` precision.
 
 ### Installation

c/lib.c

@@ -108,6 +108,8 @@ extern "C" {
     }
 
 // Dot products
+SIMSIMD_DECLARATION_DENSE(dot, i8, i8)
+SIMSIMD_DECLARATION_DENSE(dot, u8, u8)
 SIMSIMD_DECLARATION_DENSE(dot, f16, f16)
 SIMSIMD_DECLARATION_DENSE(dot, bf16, bf16)
 SIMSIMD_DECLARATION_DENSE(dot, f32, f32)
@@ -123,16 +125,19 @@ SIMSIMD_DECLARATION_DENSE(vdot, f64c, f64)
 
 // Spatial distances
 SIMSIMD_DECLARATION_DENSE(cos, i8, i8)
+SIMSIMD_DECLARATION_DENSE(cos, u8, u8)
 SIMSIMD_DECLARATION_DENSE(cos, f16, f16)
 SIMSIMD_DECLARATION_DENSE(cos, bf16, bf16)
 SIMSIMD_DECLARATION_DENSE(cos, f32, f32)
 SIMSIMD_DECLARATION_DENSE(cos, f64, f64)
 SIMSIMD_DECLARATION_DENSE(l2sq, i8, i8)
+SIMSIMD_DECLARATION_DENSE(l2sq, u8, u8)
 SIMSIMD_DECLARATION_DENSE(l2sq, f16, f16)
 SIMSIMD_DECLARATION_DENSE(l2sq, bf16, bf16)
 SIMSIMD_DECLARATION_DENSE(l2sq, f32, f32)
 SIMSIMD_DECLARATION_DENSE(l2sq, f64, f64)
 SIMSIMD_DECLARATION_DENSE(l2, i8, i8)
+SIMSIMD_DECLARATION_DENSE(l2, u8, u8)
 SIMSIMD_DECLARATION_DENSE(l2, f16, f16)
 SIMSIMD_DECLARATION_DENSE(l2, bf16, bf16)
 SIMSIMD_DECLARATION_DENSE(l2, f32, f32)
@@ -199,10 +204,13 @@ SIMSIMD_DYNAMIC simsimd_capability_t simsimd_capabilities(void) {
     void* dummy = 0;
 
     // Dense:
+    simsimd_dot_i8((simsimd_i8_t*)dummy, (simsimd_i8_t*)dummy, 0, dummy_results);
+    simsimd_dot_u8((simsimd_u8_t*)dummy, (simsimd_u8_t*)dummy, 0, dummy_results);
     simsimd_dot_f16((simsimd_f16_t*)dummy, (simsimd_f16_t*)dummy, 0, dummy_results);
     simsimd_dot_bf16((simsimd_bf16_t*)dummy, (simsimd_bf16_t*)dummy, 0, dummy_results);
     simsimd_dot_f32((simsimd_f32_t*)dummy, (simsimd_f32_t*)dummy, 0, dummy_results);
     simsimd_dot_f64((simsimd_f64_t*)dummy, (simsimd_f64_t*)dummy, 0, dummy_results);
+
     simsimd_dot_f16c((simsimd_f16_t*)dummy, (simsimd_f16_t*)dummy, 0, dummy_results);
     simsimd_dot_bf16c((simsimd_bf16_t*)dummy, (simsimd_bf16_t*)dummy, 0, dummy_results);
     simsimd_dot_f32c((simsimd_f32_t*)dummy, (simsimd_f32_t*)dummy, 0, dummy_results);
@@ -211,23 +219,32 @@ SIMSIMD_DYNAMIC simsimd_capability_t simsimd_capabilities(void) {
     simsimd_vdot_bf16c((simsimd_bf16_t*)dummy, (simsimd_bf16_t*)dummy, 0, dummy_results);
     simsimd_vdot_f32c((simsimd_f32_t*)dummy, (simsimd_f32_t*)dummy, 0, dummy_results);
     simsimd_vdot_f64c((simsimd_f64_t*)dummy, (simsimd_f64_t*)dummy, 0, dummy_results);
+
     simsimd_cos_i8((simsimd_i8_t*)dummy, (simsimd_i8_t*)dummy, 0, dummy_results);
+    simsimd_cos_u8((simsimd_u8_t*)dummy, (simsimd_u8_t*)dummy, 0, dummy_results);
     simsimd_cos_f16((simsimd_f16_t*)dummy, (simsimd_f16_t*)dummy, 0, dummy_results);
     simsimd_cos_bf16((simsimd_bf16_t*)dummy, (simsimd_bf16_t*)dummy, 0, dummy_results);
     simsimd_cos_f32((simsimd_f32_t*)dummy, (simsimd_f32_t*)dummy, 0, dummy_results);
     simsimd_cos_f64((simsimd_f64_t*)dummy, (simsimd_f64_t*)dummy, 0, dummy_results);
+
     simsimd_l2sq_i8((simsimd_i8_t*)dummy, (simsimd_i8_t*)dummy, 0, dummy_results);
+    simsimd_l2sq_u8((simsimd_u8_t*)dummy, (simsimd_u8_t*)dummy, 0, dummy_results);
     simsimd_l2sq_f16((simsimd_f16_t*)dummy, (simsimd_f16_t*)dummy, 0, dummy_results);
     simsimd_l2sq_bf16((simsimd_bf16_t*)dummy, (simsimd_bf16_t*)dummy, 0, dummy_results);
     simsimd_l2sq_f32((simsimd_f32_t*)dummy, (simsimd_f32_t*)dummy, 0, dummy_results);
     simsimd_l2sq_f64((simsimd_f64_t*)dummy, (simsimd_f64_t*)dummy, 0, dummy_results);
+
     simsimd_l2_i8((simsimd_i8_t*)dummy, (simsimd_i8_t*)dummy, 0, dummy_results);
+    simsimd_l2_i8((simsimd_i8_t*)dummy, (simsimd_i8_t*)dummy, 0, dummy_results);
+    simsimd_l2_u8((simsimd_u8_t*)dummy, (simsimd_u8_t*)dummy, 0, dummy_results);
     simsimd_l2_f16((simsimd_f16_t*)dummy, (simsimd_f16_t*)dummy, 0, dummy_results);
     simsimd_l2_bf16((simsimd_bf16_t*)dummy, (simsimd_bf16_t*)dummy, 0, dummy_results);
     simsimd_l2_f32((simsimd_f32_t*)dummy, (simsimd_f32_t*)dummy, 0, dummy_results);
     simsimd_l2_f64((simsimd_f64_t*)dummy, (simsimd_f64_t*)dummy, 0, dummy_results);
+
     simsimd_hamming_b8((simsimd_b8_t*)dummy, (simsimd_b8_t*)dummy, 0, dummy_results);
     simsimd_jaccard_b8((simsimd_b8_t*)dummy, (simsimd_b8_t*)dummy, 0, dummy_results);
+
     simsimd_kl_f16((simsimd_f16_t*)dummy, (simsimd_f16_t*)dummy, 0, dummy_results);
     simsimd_kl_bf16((simsimd_bf16_t*)dummy, (simsimd_bf16_t*)dummy, 0, dummy_results);
     simsimd_kl_f32((simsimd_f32_t*)dummy, (simsimd_f32_t*)dummy, 0, dummy_results);

include/simsimd/binary.h

@@ -319,7 +319,20 @@ SIMSIMD_PUBLIC void simsimd_jaccard_b8_ice(simsimd_b8_t const* a, simsimd_b8_t c
                                            simsimd_distance_t* result) {
 
     simsimd_size_t intersection = 0, union_ = 0;
-    // It's harder to squeeze out performance from tiny representations, so we unroll the loops for binary metrics.
+    //? On such vectors we can clearly see that the CPU struggles to perform this many parallel
+    //? population counts, because the throughput of Jaccard and Hamming in this case starts to differ.
+    //? One optimization, aside from Harley-Seal transforms can be using "shuffles" for nibble-popcount
+    //? lookups, to utilize other ports on the CPU.
+    //? https://github.com/ashvardanian/SimSIMD/pull/138
+    //
+    //  - `_mm512_popcnt_epi64` maps to `VPOPCNTQ (ZMM, K, ZMM)`:
+    //      - On Ice Lake: 3 cycles latency, ports: 1*p5
+    //      - On Genoa: 2 cycles latency, ports: 1*FP01
+    //  - `_mm512_shuffle_epi8` maps to `VPSHUFB (ZMM, ZMM, ZMM)`:
+    //      - On Ice Lake: 1 cycles latency, ports: 1*p5
+    //      - On Genoa: 2 cycles latency, ports: 1*FP12
+    //
+    //  It's harder to squeeze out performance from tiny representations, so we unroll the loops for binary metrics.
     if (n_words <= 64) { // Up to 512 bits.
         __mmask64 mask = (__mmask64)_bzhi_u64(0xFFFFFFFFFFFFFFFF, n_words);
         __m512i a_vec = _mm512_maskz_loadu_epi8(mask, a);
@@ -341,20 +354,6 @@ SIMSIMD_PUBLIC void simsimd_jaccard_b8_ice(simsimd_b8_t const* a, simsimd_b8_t c
         intersection = _mm512_reduce_add_epi64(_mm512_add_epi64(and2_count_vec, and1_count_vec));
         union_ = _mm512_reduce_add_epi64(_mm512_add_epi64(or2_count_vec, or1_count_vec));
     } else if (n_words <= 196) { // Up to 1568 bits.
-        // TODO: On such vectors we can clearly see that the CPU struggles to perform this many parallel
-        // population counts, because the throughput of Jaccard and Hamming in this case starts to differ.
-        // One optimization, aside from Harley-Seal transforms can be using "shuffles" for nibble-popcount
-        // lookups, to utilize other ports on the CPU.
-        // https://github.com/ashvardanian/SimSIMD/pull/138
-        //
-        // On Ice Lake:
-        //  - `VPOPCNTQ (ZMM, K, ZMM)` can only execute on port 5, which is a bottleneck.
-        //  - `VPSHUFB (ZMM, ZMM, ZMM)` can only run on the same port 5 as well!
-        // On Zen4:
-        //  - `VPOPCNTQ (ZMM, K, ZMM)` can run on ports: 0, 1.
-        //  - `VPSHUFB (ZMM, ZMM, ZMM)` can run on ports: 1, 2.
-        // https://uops.info/table.html?search=VPOPCNTQ%20(ZMM%2C%20K%2C%20ZMM)&cb_lat=on&cb_tp=on&cb_uops=on&cb_ports=on&cb_SKX=on&cb_ICL=on&cb_TGL=on&cb_measurements=on&cb_doc=on&cb_avx512=on
-        // https://uops.info/table.html?search=VPSHUFB%20(ZMM%2C%20ZMM%2C%20ZMM)&cb_lat=on&cb_tp=on&cb_uops=on&cb_ports=on&cb_SKX=on&cb_ICL=on&cb_TGL=on&cb_measurements=on&cb_doc=on&cb_avx512=on
         __mmask64 mask = (__mmask64)_bzhi_u64(0xFFFFFFFFFFFFFFFF, n_words - 128);
         __m512i a1_vec = _mm512_loadu_epi8(a);
         __m512i b1_vec = _mm512_loadu_epi8(b);