22//!
33//! This module provides blazingly fast bitwise Hamming distance calculation
44//! using SIMD intrinsics where available.
5+ //!
6+ //! # Optimizations
7+ //! - Branchless hex parsing with 256-byte lookup table
8+ //! - SIMD vectorized processing (AVX2/SSE4.1/NEON)
9+ //! - Batched horizontal summation to minimize lane reductions
10+ //! - Unsafe bounds elimination in hot loops
11+ //! - Algorithm selection based on input size thresholds
512
613use pyo3:: exceptions:: PyValueError ;
714use pyo3:: prelude:: * ;
@@ -10,6 +17,24 @@ use std::sync::atomic::{AtomicU8, Ordering};
1017/// Lookup table for popcount of 4-bit values (0-15)
1118const LOOKUP : [ u8 ; 16 ] = [ 0 , 1 , 1 , 2 , 1 , 2 , 2 , 3 , 1 , 2 , 2 , 3 , 2 , 3 , 3 , 4 ] ;
1219
20+ /// Branchless hex character to nibble lookup table (256 entries)
21+ /// Invalid characters map to 0xFF for easy detection
22+ const HEX_LOOKUP : [ u8 ; 256 ] = {
23+ let mut table = [ 0xFFu8 ; 256 ] ;
24+ let mut i = 0u8 ;
25+ loop {
26+ table[ i as usize ] = match i {
27+ b'0' ..=b'9' => i - b'0' ,
28+ b'A' ..=b'F' => i - b'A' + 10 ,
29+ b'a' ..=b'f' => i - b'a' + 10 ,
30+ _ => 0xFF ,
31+ } ;
32+ if i == 255 { break ; }
33+ i += 1 ;
34+ }
35+ table
36+ } ;
37+
1338/// Algorithm selection constants
1439const ALGO_CLASSIC : u8 = 0 ;
1540const ALGO_NATIVE : u8 = 1 ;
@@ -20,11 +45,15 @@ const ALGO_AVX2: u8 = 3;
2045#[ cfg( target_arch = "aarch64" ) ]
2146const ALGO_NEON : u8 = 4 ;
2247
48+ /// Thresholds for algorithm selection (tuned for typical CPU cache behavior)
49+ const SCALAR_THRESHOLD : usize = 16 ; // Below this, scalar may beat SIMD
50+ const SSE_THRESHOLD : usize = 64 ; // Use SSE for medium strings
51+
2352/// Current algorithm selection (global state)
2453static CURRENT_ALGO : AtomicU8 = AtomicU8 :: new ( ALGO_NATIVE ) ;
2554
26- /// Classic popcount implementation using bit manipulation
27- #[ inline]
55+ /// Classic popcount implementation using bit manipulation (Wilkes-Wheeler-Gill)
56+ #[ inline( always ) ]
2857fn popcnt64_classic ( mut x : u64 ) -> u64 {
2958 const M1 : u64 = 0x5555555555555555 ;
3059 const M2 : u64 = 0x3333333333333333 ;
@@ -37,70 +66,152 @@ fn popcnt64_classic(mut x: u64) -> u64 {
3766}
3867
3968/// Native popcount implementation using CPU instruction
40- #[ inline]
69+ #[ inline( always ) ]
4170fn popcnt64_native ( x : u64 ) -> u64 {
4271 x. count_ones ( ) as u64
4372}
4473
74+ /// Branchless hex character to nibble conversion using lookup table
75+ /// Returns 0xFF for invalid characters
76+ #[ inline( always) ]
77+ fn hex_char_to_nibble ( c : u8 ) -> u8 {
78+ // SAFETY: c is u8, so always in bounds of 256-element table
79+ unsafe { * HEX_LOOKUP . get_unchecked ( c as usize ) }
80+ }
81+
4582/// Convert a hex character to its numeric value (0-15)
4683/// Returns None if the character is not a valid hex digit
47- #[ inline]
84+ #[ inline( always ) ]
4885fn hex_char_to_val ( c : u8 ) -> Option < u8 > {
49- match c {
50- b'0' ..=b'9' => Some ( c - b'0' ) ,
51- b'A' ..=b'F' => Some ( c - b'A' + 10 ) ,
52- b'a' ..=b'f' => Some ( c - b'a' + 10 ) ,
53- _ => None ,
54- }
86+ let val = hex_char_to_nibble ( c) ;
87+ if val == 0xFF { None } else { Some ( val) }
5588}
5689
5790/// Calculate hamming distance between two hex strings using classic algorithm
91+ /// Optimized with branchless lookup and bounds check elimination
92+ #[ inline( always) ]
5893fn hamming_distance_string_classic ( a : & [ u8 ] , b : & [ u8 ] ) -> Result < u64 , & ' static str > {
94+ let len = a. len ( ) ;
5995 let mut result: u64 = 0 ;
60- for i in 0 ..a. len ( ) {
61- let val1 = hex_char_to_val ( a[ i] ) . ok_or ( "hex string contains invalid char" ) ?;
62- let val2 = hex_char_to_val ( b[ i] ) . ok_or ( "hex string contains invalid char" ) ?;
63- result += LOOKUP [ ( val1 ^ val2) as usize ] as u64 ;
96+ let mut i = 0 ;
97+
98+ // Process 4 hex chars at a time to reduce loop overhead
99+ while i + 4 <= len {
100+ // SAFETY: i + 3 < len verified by loop condition
101+ unsafe {
102+ let val1_0 = hex_char_to_nibble ( * a. get_unchecked ( i) ) ;
103+ let val2_0 = hex_char_to_nibble ( * b. get_unchecked ( i) ) ;
104+ let val1_1 = hex_char_to_nibble ( * a. get_unchecked ( i + 1 ) ) ;
105+ let val2_1 = hex_char_to_nibble ( * b. get_unchecked ( i + 1 ) ) ;
106+ let val1_2 = hex_char_to_nibble ( * a. get_unchecked ( i + 2 ) ) ;
107+ let val2_2 = hex_char_to_nibble ( * b. get_unchecked ( i + 2 ) ) ;
108+ let val1_3 = hex_char_to_nibble ( * a. get_unchecked ( i + 3 ) ) ;
109+ let val2_3 = hex_char_to_nibble ( * b. get_unchecked ( i + 3 ) ) ;
110+
111+ // Check all 8 values for validity (0xFF indicates invalid)
112+ // Use bitwise OR to combine checks - any 0xFF will result in high bit set
113+ let invalid = ( val1_0 | val2_0 | val1_1 | val2_1 | val1_2 | val2_2 | val1_3 | val2_3) & 0xF0 ;
114+ if invalid != 0 {
115+ return Err ( "hex string contains invalid char" ) ;
116+ }
117+
118+ result += * LOOKUP . get_unchecked ( ( val1_0 ^ val2_0) as usize ) as u64
119+ + * LOOKUP . get_unchecked ( ( val1_1 ^ val2_1) as usize ) as u64
120+ + * LOOKUP . get_unchecked ( ( val1_2 ^ val2_2) as usize ) as u64
121+ + * LOOKUP . get_unchecked ( ( val1_3 ^ val2_3) as usize ) as u64 ;
122+ }
123+ i += 4 ;
124+ }
125+
126+ // Handle remaining characters
127+ while i < len {
128+ // SAFETY: i < len verified by loop condition
129+ unsafe {
130+ let val1 = hex_char_to_nibble ( * a. get_unchecked ( i) ) ;
131+ let val2 = hex_char_to_nibble ( * b. get_unchecked ( i) ) ;
132+ if ( val1 | val2) & 0xF0 != 0 {
133+ return Err ( "hex string contains invalid char" ) ;
134+ }
135+ result += * LOOKUP . get_unchecked ( ( val1 ^ val2) as usize ) as u64 ;
136+ }
137+ i += 1 ;
64138 }
139+
65140 Ok ( result)
66141}
67142
68143/// Calculate hamming distance between two byte arrays using classic algorithm
144+ /// Optimized with loop unrolling and bounds check elimination
145+ #[ inline( always) ]
69146fn hamming_distance_bytes_classic ( a : & [ u8 ] , b : & [ u8 ] , max_dist : i64 ) -> u64 {
70- let mut difference: u64 = 0 ;
71147 let length = a. len ( ) ;
72148
73149 if max_dist < 0 {
74- // Full distance calculation
150+ // Full distance calculation - heavily optimized
151+ let mut difference: u64 = 0 ;
75152 let mut i = 0 ;
76- // Process 8 bytes at a time
153+
154+ // Process 32 bytes at a time (4 x 8-byte chunks)
155+ while i + 32 <= length {
156+ // SAFETY: i + 31 < length verified by loop condition
157+ unsafe {
158+ let a0 = u64:: from_ne_bytes ( * ( a. as_ptr ( ) . add ( i) as * const [ u8 ; 8 ] ) ) ;
159+ let b0 = u64:: from_ne_bytes ( * ( b. as_ptr ( ) . add ( i) as * const [ u8 ; 8 ] ) ) ;
160+ let a1 = u64:: from_ne_bytes ( * ( a. as_ptr ( ) . add ( i + 8 ) as * const [ u8 ; 8 ] ) ) ;
161+ let b1 = u64:: from_ne_bytes ( * ( b. as_ptr ( ) . add ( i + 8 ) as * const [ u8 ; 8 ] ) ) ;
162+ let a2 = u64:: from_ne_bytes ( * ( a. as_ptr ( ) . add ( i + 16 ) as * const [ u8 ; 8 ] ) ) ;
163+ let b2 = u64:: from_ne_bytes ( * ( b. as_ptr ( ) . add ( i + 16 ) as * const [ u8 ; 8 ] ) ) ;
164+ let a3 = u64:: from_ne_bytes ( * ( a. as_ptr ( ) . add ( i + 24 ) as * const [ u8 ; 8 ] ) ) ;
165+ let b3 = u64:: from_ne_bytes ( * ( b. as_ptr ( ) . add ( i + 24 ) as * const [ u8 ; 8 ] ) ) ;
166+
167+ difference += popcnt64_classic ( a0 ^ b0)
168+ + popcnt64_classic ( a1 ^ b1)
169+ + popcnt64_classic ( a2 ^ b2)
170+ + popcnt64_classic ( a3 ^ b3) ;
171+ }
172+ i += 32 ;
173+ }
174+
175+ // Process remaining 8-byte chunks
77176 while i + 8 <= length {
78- let a_chunk = u64:: from_ne_bytes ( a[ i..i + 8 ] . try_into ( ) . unwrap ( ) ) ;
79- let b_chunk = u64:: from_ne_bytes ( b[ i..i + 8 ] . try_into ( ) . unwrap ( ) ) ;
80- difference += popcnt64_classic ( a_chunk ^ b_chunk) ;
177+ unsafe {
178+ let a_chunk = u64:: from_ne_bytes ( * ( a. as_ptr ( ) . add ( i) as * const [ u8 ; 8 ] ) ) ;
179+ let b_chunk = u64:: from_ne_bytes ( * ( b. as_ptr ( ) . add ( i) as * const [ u8 ; 8 ] ) ) ;
180+ difference += popcnt64_classic ( a_chunk ^ b_chunk) ;
181+ }
81182 i += 8 ;
82183 }
184+
83185 // Process remaining bytes
84186 while i < length {
85- difference += popcnt64_classic ( ( a[ i] ^ b[ i] ) as u64 ) ;
187+ unsafe {
188+ difference += popcnt64_classic ( ( * a. get_unchecked ( i) ^ * b. get_unchecked ( i) ) as u64 ) ;
189+ }
86190 i += 1 ;
87191 }
88192 difference
89193 } else {
90194 // Early termination if exceeds max_dist
195+ let max_dist_u64 = max_dist as u64 ;
196+ let mut difference: u64 = 0 ;
91197 let mut i = 0 ;
198+
92199 while i + 8 <= length {
93- let a_chunk = u64:: from_ne_bytes ( a[ i..i + 8 ] . try_into ( ) . unwrap ( ) ) ;
94- let b_chunk = u64:: from_ne_bytes ( b[ i..i + 8 ] . try_into ( ) . unwrap ( ) ) ;
95- difference += popcnt64_classic ( a_chunk ^ b_chunk) ;
96- if difference > max_dist as u64 {
200+ unsafe {
201+ let a_chunk = u64:: from_ne_bytes ( * ( a. as_ptr ( ) . add ( i) as * const [ u8 ; 8 ] ) ) ;
202+ let b_chunk = u64:: from_ne_bytes ( * ( b. as_ptr ( ) . add ( i) as * const [ u8 ; 8 ] ) ) ;
203+ difference += popcnt64_classic ( a_chunk ^ b_chunk) ;
204+ }
205+ if difference > max_dist_u64 {
97206 return 0 ;
98207 }
99208 i += 8 ;
100209 }
101210 while i < length {
102- difference += popcnt64_classic ( ( a[ i] ^ b[ i] ) as u64 ) ;
103- if difference > max_dist as u64 {
211+ unsafe {
212+ difference += popcnt64_classic ( ( * a. get_unchecked ( i) ^ * b. get_unchecked ( i) ) as u64 ) ;
213+ }
214+ if difference > max_dist_u64 {
104215 return 0 ;
105216 }
106217 i += 1 ;
@@ -110,37 +221,74 @@ fn hamming_distance_bytes_classic(a: &[u8], b: &[u8], max_dist: i64) -> u64 {
110221}
111222
112223/// Calculate hamming distance between two byte arrays using native popcount
224+ /// Optimized with aggressive loop unrolling and bounds check elimination
225+ #[ inline( always) ]
113226fn hamming_distance_bytes_native ( a : & [ u8 ] , b : & [ u8 ] , max_dist : i64 ) -> u64 {
114- let mut difference: u64 = 0 ;
115227 let length = a. len ( ) ;
116228
117229 if max_dist < 0 {
230+ let mut difference: u64 = 0 ;
118231 let mut i = 0 ;
232+
233+ // Process 32 bytes at a time (4 x 8-byte chunks) to saturate execution units
234+ while i + 32 <= length {
235+ unsafe {
236+ let a0 = u64:: from_ne_bytes ( * ( a. as_ptr ( ) . add ( i) as * const [ u8 ; 8 ] ) ) ;
237+ let b0 = u64:: from_ne_bytes ( * ( b. as_ptr ( ) . add ( i) as * const [ u8 ; 8 ] ) ) ;
238+ let a1 = u64:: from_ne_bytes ( * ( a. as_ptr ( ) . add ( i + 8 ) as * const [ u8 ; 8 ] ) ) ;
239+ let b1 = u64:: from_ne_bytes ( * ( b. as_ptr ( ) . add ( i + 8 ) as * const [ u8 ; 8 ] ) ) ;
240+ let a2 = u64:: from_ne_bytes ( * ( a. as_ptr ( ) . add ( i + 16 ) as * const [ u8 ; 8 ] ) ) ;
241+ let b2 = u64:: from_ne_bytes ( * ( b. as_ptr ( ) . add ( i + 16 ) as * const [ u8 ; 8 ] ) ) ;
242+ let a3 = u64:: from_ne_bytes ( * ( a. as_ptr ( ) . add ( i + 24 ) as * const [ u8 ; 8 ] ) ) ;
243+ let b3 = u64:: from_ne_bytes ( * ( b. as_ptr ( ) . add ( i + 24 ) as * const [ u8 ; 8 ] ) ) ;
244+
245+ difference += popcnt64_native ( a0 ^ b0)
246+ + popcnt64_native ( a1 ^ b1)
247+ + popcnt64_native ( a2 ^ b2)
248+ + popcnt64_native ( a3 ^ b3) ;
249+ }
250+ i += 32 ;
251+ }
252+
253+ // Process remaining 8-byte chunks
119254 while i + 8 <= length {
120- let a_chunk = u64:: from_ne_bytes ( a[ i..i + 8 ] . try_into ( ) . unwrap ( ) ) ;
121- let b_chunk = u64:: from_ne_bytes ( b[ i..i + 8 ] . try_into ( ) . unwrap ( ) ) ;
122- difference += popcnt64_native ( a_chunk ^ b_chunk) ;
255+ unsafe {
256+ let a_chunk = u64:: from_ne_bytes ( * ( a. as_ptr ( ) . add ( i) as * const [ u8 ; 8 ] ) ) ;
257+ let b_chunk = u64:: from_ne_bytes ( * ( b. as_ptr ( ) . add ( i) as * const [ u8 ; 8 ] ) ) ;
258+ difference += popcnt64_native ( a_chunk ^ b_chunk) ;
259+ }
123260 i += 8 ;
124261 }
262+
263+ // Process remaining bytes
125264 while i < length {
126- difference += popcnt64_native ( ( a[ i] ^ b[ i] ) as u64 ) ;
265+ unsafe {
266+ difference += ( * a. get_unchecked ( i) ^ * b. get_unchecked ( i) ) . count_ones ( ) as u64 ;
267+ }
127268 i += 1 ;
128269 }
129270 difference
130271 } else {
272+ let max_dist_u64 = max_dist as u64 ;
273+ let mut difference: u64 = 0 ;
131274 let mut i = 0 ;
275+
132276 while i + 8 <= length {
133- let a_chunk = u64:: from_ne_bytes ( a[ i..i + 8 ] . try_into ( ) . unwrap ( ) ) ;
134- let b_chunk = u64:: from_ne_bytes ( b[ i..i + 8 ] . try_into ( ) . unwrap ( ) ) ;
135- difference += popcnt64_native ( a_chunk ^ b_chunk) ;
136- if difference > max_dist as u64 {
277+ unsafe {
278+ let a_chunk = u64:: from_ne_bytes ( * ( a. as_ptr ( ) . add ( i) as * const [ u8 ; 8 ] ) ) ;
279+ let b_chunk = u64:: from_ne_bytes ( * ( b. as_ptr ( ) . add ( i) as * const [ u8 ; 8 ] ) ) ;
280+ difference += popcnt64_native ( a_chunk ^ b_chunk) ;
281+ }
282+ if difference > max_dist_u64 {
137283 return 0 ;
138284 }
139285 i += 8 ;
140286 }
141287 while i < length {
142- difference += popcnt64_native ( ( a[ i] ^ b[ i] ) as u64 ) ;
143- if difference > max_dist as u64 {
288+ unsafe {
289+ difference += ( * a. get_unchecked ( i) ^ * b. get_unchecked ( i) ) . count_ones ( ) as u64 ;
290+ }
291+ if difference > max_dist_u64 {
144292 return 0 ;
145293 }
146294 i += 1 ;
0 commit comments