Speed up batch scans and Python calls (#53)

* perf(array): stop best scan at exact match — 99.5% on 512x16

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Copilot-Session: e7b02886-a2c5-43cd-bb97-bfa037d3b625

* perf(array): bound Rayon fan-out — 84.8% on 100000x128 best

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Copilot-Session: e7b02886-a2c5-43cd-bb97-bfa037d3b625

* perf(python): bypass PyBuffer for bytes — 64.2% on within-distance

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Copilot-Session: e7b02886-a2c5-43cd-bb97-bfa037d3b625

* perf(python): keep small buffers attached — 22.7% on bytearray best

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Copilot-Session: e7b02886-a2c5-43cd-bb97-bfa037d3b625

* perf(rust): inline hex entrypoint — 21.8% on 16-char

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Copilot-Session: e7b02886-a2c5-43cd-bb97-bfa037d3b625

* perf(neon): check first hex block eagerly — 49.9% on tight 1024-char

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Copilot-Session: e7b02886-a2c5-43cd-bb97-bfa037d3b625

* perf(neon): scale eager blocks with cutoff — 47.3% on random max100

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Copilot-Session: e7b02886-a2c5-43cd-bb97-bfa037d3b625

* perf(aarch64): prefer native byte kernel — 24.0% on 64-byte Rust

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Copilot-Session: e7b02886-a2c5-43cd-bb97-bfa037d3b625

* docs: refresh M4 Max performance benchmarks

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Copilot-Session: e7b02886-a2c5-43cd-bb97-bfa037d3b625

---------

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>

Author: mrecachinas

README.rst

@@ -182,70 +182,81 @@ immutable that is a very slow operation. Use a ``bytearray`` instead, and cast i
 Benchmark
 ---------
 
-All benchmarks on Apple M-series (ARM64) with hexhamming v3.0.0, ``rustc`` 1.85, Python 3.14.
+All benchmarks were run on an Apple M4 Max (ARM64, 16 logical cores, 64 GiB)
+with hexhamming v3.0.0, ``rustc`` 1.96.1, and Python 3.14.6. Values are the
+median of the means from three independent runs.
 
 Raw Rust (no Python overhead)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 These numbers show the pure computation time using Rust's ``criterion`` benchmarks
 (``cargo bench --no-default-features``), with no Python/PyO3 overhead.
 
-===========================================  ===========
-Name                                           Mean (ns)
-===========================================  ===========
-hex_string (NEON) [16 chars]                        2.0
-hex_string (NEON) [64 chars]                        5.8
-hex_string (NEON) [128 chars]                      11.2
-hex_string (NEON) [254 chars]                      21.2
-bytes (NEON) [8 bytes]                              1.8
-bytes (NEON) [32 bytes]                             2.8
-bytes (NEON) [64 bytes]                             2.6
-bytes (NEON) [127 bytes]                            6.6
-bytes_within_dist [127 bytes]                       1.6
-array first [512×16, at start]                      2.3
-array first [512×16, at end]                      671.4
-array best [512×16]                               925.8
-array all [512×16]                                799.4
-===========================================  ===========
-
-Larger array workloads cross the parallel (Rayon) threshold; see the Python
-table below for representative end-to-end numbers.
+================================================  ===========
+Name                                              Mean (ns)
+================================================  ===========
+hex_string (NEON) [16 chars]                           2.4
+hex_string (NEON) [64 chars]                           8.3
+hex_string (NEON) [128 chars]                         16.2
+hex_string (NEON) [254 chars]                         30.2
+bytes (native) [8 bytes]                               1.7
+bytes (native) [32 bytes]                              2.4
+bytes (native) [64 bytes]                              3.2
+bytes (native) [127 bytes]                             8.4
+bytes_within_dist [127 bytes]                          2.4
+array first [512×16, at start]                         6.6
+array first [512×16, at end]                       1,397.0
+array best [512×16, exact at start]                     8.3
+array best [512×16, exact at end]                   1,599.2
+array all [512×16]                                  1,610.1
+array best [16384×64, match at mid]                71,121.0
+array all [16384×64, match at mid]                 79,365.0
+array best [100000×128, parallel]                  50,996.0
+array all [100000×128, parallel]                  144,800.0
+================================================  ===========
+
+On AArch64, LLVM's auto-vectorized native byte loop is faster than the
+hand-written NEON byte kernel for these sizes, while hexadecimal strings still
+use the packed NEON implementation. Large array workloads use four balanced
+Rayon jobs to avoid oversubscribing the memory-bound scan.
 
 Python API (via PyO3)
 ~~~~~~~~~~~~~~~~~~~~~
 
-These numbers include Python function call overhead (~45 ns) using ``pytest-benchmark``.
+These numbers include Python wrapper and function-call overhead using
+``pytest-benchmark``.
 
 ======================================================  ===========
 Name                                                      Mean (ns)
 ======================================================  ===========
-hamming_distance_string [3 chars, same]                       84.6
-hamming_distance_string [3 chars, diff]                       83.6
-hamming_distance_string [64 chars, diff]                      91.3
-hamming_distance_string [1024 chars, diff]                   207.4
-hamming_distance_bytes [3 bytes, same]                       127.3
-hamming_distance_bytes [3 bytes, diff]                        96.2
-hamming_distance_bytes [64 bytes, diff]                      169.9
-hamming_distance_bytes [1024 bytes, diff]                    175.2
-check_hexstrings_within_dist [1000 chars]                    221.3
-check_bytes_within_dist [16 bytes]                           146.2
-check_bytes_within_dist [64 bytes]                           107.7
-check_bytes_within_dist [127 bytes]                          100.6
-first_within_dist [512×16, at start]                          98.7
-first_within_dist [512×16, mid]                              570.8
-first_within_dist [512×16, at end]                         1,030.1
-first_within_dist [16384×64, at start]                        99.6
-first_within_dist [16384×64, mid]                         20,838.8
-first_within_dist [16384×64, at end]                      41,489.5
-best_within_dist [512×16, at start]                        1,609.8
-best_within_dist [512×16, at end]                          1,116.4
-best_within_dist [16384×64, mid]                          46,826.1
-all_within_dist [512×16, at start]                         1,342.9
-all_within_dist [512×16, at end]                           1,365.9
-all_within_dist [16384×64, mid]                           48,067.2
+hamming_distance_string [3 chars, same]                       56.3
+hamming_distance_string [3 chars, diff]                      105.8
+hamming_distance_string [64 chars, diff]                      60.0
+hamming_distance_string [1024 chars, diff]                   177.1
+hamming_distance_bytes [3 bytes, same]                        51.7
+hamming_distance_bytes [3 bytes, diff]                        51.8
+hamming_distance_bytes [64 bytes, diff]                       51.9
+hamming_distance_bytes [1024 bytes, diff]                     68.9
+check_hexstrings_within_dist [1000 chars]                     56.4
+check_bytes_within_dist [16 bytes]                            52.5
+check_bytes_within_dist [64 bytes]                            51.8
+check_bytes_within_dist [127 bytes]                           52.8
+first_within_dist [512×16, at start]                          58.5
+first_within_dist [512×16, mid]                              771.4
+first_within_dist [512×16, at end]                         1,475.1
+first_within_dist [16384×64, at start]                       160.3
+first_within_dist [16384×64, mid]                         23,031.5
+first_within_dist [16384×64, at end]                      45,801.2
+best_within_dist [512×16, at start]                           75.3
+best_within_dist [512×16, at end]                          1,703.5
+best_within_dist [16384×64, mid]                          93,212.8
+all_within_dist [512×16, at start]                         1,735.5
+all_within_dist [512×16, at end]                           1,747.3
+all_within_dist [16384×64, mid]                           93,056.3
 ======================================================  ===========
 
-For small inputs, Python call overhead dominates (~45 ns). For large inputs
+For small inputs, Python call and wrapper overhead dominates (roughly 40–55 ns
+on this machine). For large inputs
 (1024+ chars, 16384-element arrays), computation dominates and Python overhead
 is negligible. Array APIs transparently parallelize with Rayon once the input
 exceeds ~64 KiB; the ``first`` variant additionally short-circuits on the first

benches/hamming_bench.rs

@@ -47,7 +47,7 @@ fn bench_bytes_by_algo(c: &mut Criterion) {
     let algos: &[&str] = if cfg!(target_arch = "x86_64") {
         &["classic", "sse", "avx2", "avx512"]
     } else if cfg!(target_arch = "aarch64") {
-        &["classic", "neon"]
+        &["classic", "native", "neon"]
     } else {
         &["classic", "native"]
     };

src/api.rs

@@ -9,7 +9,27 @@ use rayon::prelude::*;
 use std::sync::atomic::Ordering;
 
 /// Minimum total byte size of big_array before we use rayon parallel paths.
-const PAR_THRESHOLD_BYTES: usize = 64 * 1024;
+const PAR_THRESHOLD_BYTES: usize = 256 * 1024;
+/// Keep byte-array scans to a small number of coarse jobs. More workers spend
+/// more time scheduling these very small per-record calculations than running
+/// them on current many-core CPUs.
+const PAR_JOBS: usize = 4;
+
+#[inline]
+fn partition_element_ranges(num_elements: usize) -> [(usize, usize); PAR_JOBS] {
+    let base = num_elements / PAR_JOBS;
+    let remainder = num_elements % PAR_JOBS;
+    let mut ranges = [(0, 0); PAR_JOBS];
+    let mut start = 0;
+
+    for (job, range) in ranges.iter_mut().enumerate() {
+        let end = start + base + usize::from(job < remainder);
+        *range = (start, end);
+        start = end;
+    }
+
+    ranges
+}
 
 /// Calculate the bitwise hamming distance between two equal-length hex strings.
 ///
@@ -23,6 +43,7 @@ const PAR_THRESHOLD_BYTES: usize = 64 * 1024;
 /// let dist = hexhamming::hex_hamming_distance("deadbeef", "00000000").unwrap();
 /// assert_eq!(dist, 24);
 /// ```
+#[inline]
 pub fn hex_hamming_distance(a: &str, b: &str) -> Result<u64, &'static str> {
     if a.len() != b.len() {
         return Err("strings are NOT the same length");
@@ -121,43 +142,27 @@ pub fn bytes_array_best_within_dist(
         return Ok(serial_best_within_dist(big_array, small_array, max_dist));
     }
     let elem_size = small_array.len();
-    Ok(big_array
-        .par_chunks_exact(elem_size)
-        .enumerate()
-        .fold(
-            || None::<(u64, usize)>,
-            |acc, (i, chunk)| {
-                let threshold = acc
-                    .map(|(d, _)| (d as i64).saturating_sub(1))
-                    .unwrap_or(max_dist);
-                let d = hamming_distance_bytes_dispatch(chunk, small_array, threshold);
-                if d == u64::MAX {
-                    return acc;
-                }
-                match acc {
-                    None => Some((d, i)),
-                    Some((best_d, _)) if d < best_d => Some((d, i)),
-                    _ => acc,
-                }
-            },
-        )
-        .reduce(
-            || None,
-            |a, b| match (a, b) {
-                (None, x) | (x, None) => x,
-                (Some(x), Some(y)) => {
-                    if x.0 < y.0 {
-                        Some(x)
-                    } else if y.0 < x.0 {
-                        Some(y)
-                    } else if x.1 < y.1 {
-                        Some(x)
-                    } else {
-                        Some(y)
-                    }
-                }
-            },
-        ))
+    let num_elements = big_array.len() / elem_size;
+    let ranges = partition_element_ranges(num_elements);
+
+    Ok(ranges
+        .par_iter()
+        .with_max_len(1)
+        .map(|&(start, end)| {
+            let chunk = &big_array[start * elem_size..end * elem_size];
+            serial_best_within_dist(chunk, small_array, max_dist)
+                .map(|(distance, index)| (distance, index + start))
+        })
+        .reduce(|| None, merge_best))
+}
+
+#[inline]
+fn merge_best(a: Option<(u64, usize)>, b: Option<(u64, usize)>) -> Option<(u64, usize)> {
+    match (a, b) {
+        (None, x) | (x, None) => x,
+        (Some(x), Some(y)) if x.0 < y.0 || (x.0 == y.0 && x.1 <= y.1) => Some(x),
+        (Some(_), Some(y)) => Some(y),
+    }
 }
 
 #[inline]
@@ -180,6 +185,9 @@ fn serial_best_within_dist(
         }
         if best.is_none() || d < best.unwrap().0 {
             best = Some((d, i));
+            if d == 0 {
+                return best;
+            }
         }
     }
     best
@@ -203,19 +211,27 @@ pub fn bytes_array_all_within_dist(
         return Ok(serial_all_within_dist(big_array, small_array, max_dist));
     }
     let elem_size = small_array.len();
-    let mut results: Vec<(u64, usize)> = big_array
-        .par_chunks_exact(elem_size)
-        .enumerate()
-        .filter_map(|(i, chunk)| {
-            let d = hamming_distance_bytes_dispatch(chunk, small_array, max_dist);
-            if d == u64::MAX {
-                None
-            } else {
-                Some((d, i))
-            }
+    let num_elements = big_array.len() / elem_size;
+    let ranges = partition_element_ranges(num_elements);
+    let per_job: Vec<Vec<(u64, usize)>> = ranges
+        .par_iter()
+        .with_max_len(1)
+        .map(|&(start, end)| {
+            let chunk = &big_array[start * elem_size..end * elem_size];
+            serial_all_within_dist(chunk, small_array, max_dist)
+                .into_iter()
+                .map(|(distance, index)| (distance, index + start))
+                .collect()
         })
         .collect();
-    results.sort_unstable_by_key(|&(_, idx)| idx);
+
+    // The indexed parallel iterator preserves range order, and each serial
+    // result is already ordered, so flattening preserves ascending indices.
+    let result_count = per_job.iter().map(Vec::len).sum();
+    let mut results = Vec::with_capacity(result_count);
+    for job_results in per_job {
+        results.extend(job_results);
+    }
     Ok(results)
 }
 
@@ -423,11 +439,25 @@ mod tests {
         big
     }
 
+    #[test]
+    fn parallel_ranges_cover_elements_once() {
+        for num_elements in [1, 3, 4, 5, 7, 8, 9, 100, 100_001] {
+            let ranges = partition_element_ranges(num_elements);
+            assert_eq!(ranges[0].0, 0);
+            assert_eq!(ranges[PAR_JOBS - 1].1, num_elements);
+            for pair in ranges.windows(2) {
+                assert_eq!(pair[0].1, pair[1].0);
+            }
+            let lengths: Vec<usize> = ranges.iter().map(|&(start, end)| end - start).collect();
+            assert!(lengths.iter().max().unwrap() - lengths.iter().min().unwrap() <= 1);
+        }
+    }
+
     #[test]
     fn first_within_dist_small_batch() {
         // Below PAR_THRESHOLD_BYTES → serial path
         let elem_size = 4;
-        let n = 100; // 400 bytes < 64KB
+        let n = 100; // 400 bytes < parallel threshold
         let needle = vec![0x00u8; elem_size];
         let big = make_batch(elem_size, n, 0xFF, &[50], &needle);
         assert_eq!(

src/neon_simd.rs

@@ -481,10 +481,40 @@ pub unsafe fn hamming_distance_string_neon_pack_with_max(
     let mut difference: u64 = 0;
     let mut bad_acc = zero;
 
-    // Process 32 hex chars at a time, batching cross-lane reductions and the
-    // threshold check. Each packed byte popcount is ≤8, so BATCH iterations are
-    // safe in u8 lanes (BATCH*8 < 256). Early exit is granular to one batch.
-    const BATCH: usize = 16;
+    // Check enough initial blocks to cross a tight threshold under typical
+    // two-bits-per-nibble data. Larger thresholds go straight to batching.
+    let eager_blocks = if max_dist < 256 {
+        (max_dist as usize / 64) + 1
+    } else {
+        0
+    };
+    let mut eager = 0;
+    while eager < eager_blocks && i + 32 <= length {
+        let (packed, bad) = pack32_xor_neon(
+            a.as_ptr().add(i),
+            b.as_ptr().add(i),
+            case_mask,
+            ascii_0,
+            seven,
+            nine,
+            ten,
+            fifteen_u,
+        );
+        bad_acc = vorrq_u8(bad_acc, bad);
+        if vmaxvq_u8(bad_acc) != 0 {
+            return Err("hex string contains invalid char");
+        }
+        difference += vaddlvq_u8(vcntq_u8(packed)) as u64;
+        if difference > max_dist {
+            return Ok(u64::MAX);
+        }
+        i += 32;
+        eager += 1;
+    }
+
+    // Each packed byte popcount is at most 8, so 31 iterations fit safely in
+    // u8 lanes (31 * 8 = 248) before one cross-lane reduction.
+    const BATCH: usize = 31;
     while i + 32 <= length {
         let mut acc = zero;
         let mut n = 0;