Skip to main content

lev/
lib.rs

1//! An extremely fast Levenshtein distance and similarity ratio.
2//!
3//! * [`distance`] – the Levenshtein edit distance between two strings.
4//! * [`ratio`]    – the indel similarity score in `[0.0, 1.0]`.
5//!
6//! ```
7//! assert_eq!(lev::distance("kitten", "sitting"), 3);
8//! assert_eq!(lev::ratio("abc", "abc"), 1.0);
9//! ```
10//!
11//! Characters are counted as Unicode scalar values (`char`s), so an accented
12//! letter or an emoji counts as one character.
13//!
14//! Both functions use bit-parallel algorithms (Myers 1999, Hyyrö 2003)
15//! after stripping common prefixes and suffixes, and a diagonal band for
16//! long, similar strings.  On x86-64 the hot path is compiled a second time
17//! for x86-64-v2 and picked at run time.
18//!
19//! The same code backs the [`lev-rs`](https://pypi.org/project/lev-rs/)
20//! Python package.  Its bindings are behind the `python` feature, which
21//! only the Python build enables; Python 3.10 through 3.15 are supported,
22//! including free-threaded builds (PEP 703).
23
24// # Algorithmic Architecture
25//
26// See <https://en.wikipedia.org/wiki/Levenshtein_distance> for the distance
27// definition.  Five performance layers are composed to minimize constant factors
28// across all string length and similarity regimes:
29//
30// 1. **Identity short-circuit** (Python only) – if both arguments point to the
31//    identical Python string object (`s1.is(s2)`), return `0` (or `1.0` for
32//    `ratio`) immediately without examining string contents.
33// 2. **Common-affix stripping** – shared leading and trailing code units are
34//    stripped before the main computation using vectorized 16-element/32-byte
35//    block comparisons followed by scalar word scans (`strip_affix`, `strip_affix_mixed`).
36// 3. **Zero-copy CPython buffer access** – Python stores strings in one of
37//    three compact internal encodings (UCS-1 / UCS-2 / UCS-4); see
38//    [PEP 393](https://peps.python.org/pep-0393/).  Kind, ascii flag, data
39//    pointer, and length are extracted from the `PyASCIIObject` header into a
40//    [`UniView`] in a single pass without copying; subsequent dispatch operates
41//    exclusively on local slices.  The Rust API borrows pure-ASCII `&str`s as
42//    they are and transcodes any other string into the narrowest of the three
43//    encodings that holds it (see [`Units`]).
44//    - *UCS-1* (`u8`, ≤ U+00FF): peq is a flat `[u64; 128]` (pure ASCII) or
45//      `[u64; 256]` (Latin-1) stack array — O(1) direct-index lookup.
46//    - *UCS-2* (`u16`) and *UCS-4* (`u32`): peq is a 128-slot stack-allocated
47//      open-addressing hash table with Fibonacci hashing — O(1) amortized
48//      lookup at ≤ 50% load.  Mixed-kind pairs compare code units in place;
49//      no UTF-8 transcoding or allocations, preserving lone surrogates.
50// 4. **Bit-parallel algorithm (Hyyrö/Myers) & mbleven fast path** – O(⌈m/w⌉ · n)
51//    time with w = 64 (Hyyrö 2003, Myers 1999):
52//    - *Tiny patterns* (`m ≤ 16`): table-free execution; the match word is
53//      computed per text character using SWAR for byte patterns up to length 8,
54//      128-bit SIMD (NEON on AArch64, SSE2 on x86-64) for byte patterns 9..=16,
55//      or unrolled scalar scans for wider code units.
56//    - *Small-distance fast path* (`ub ≤ 3`, `m ≤ 512`): for strings differing by
57//      very few edits, a SIMD Hamming scan (`small_ub`) establishes an upper bound
58//      and gates into `mbleven` (Hyyrö 2018), evaluating candidate branch patterns
59//      in O(ub · n) time without building peq tables or allocating memory.
60//    - *Single-word* (`m ≤ 64`): one 64-bit word covers the entire pattern.
61//    - *Multi-word* (`64 < m ≤ 512`): ⌈m/64⌉ words with carry propagation.
62//      UCS-1 uses a `[[u64; W]; 128 or 256]` stack peq (O(1) direct lookup, no heap).
63//      UCS-2/4 use stack-allocated open-addressing hash tables.
64//    - *Large strings* (`m > 512`): one loop over all words per column; on
65//      AArch64, strips of eight words sweep the text one after another
66//      with the const-width kernel instead, handing the per-column carry
67//      and boundary bits upwards, so each strip's state stays in registers.
68//      The LCS full matrix runs in strips of four words on x86-64 too.
69//    - Steps with target-specific forms live in [`kernels`], one module per
70//      target.  On x86-64 the dispatchers pick a form at runtime by CPU
71//      feature detection; on AArch64 the form is fixed at compile time by
72//      `cfg` (Apple or Neoverse tuning).
73// 5. **Ukkonen diagonal banding (`m > 512`) & LCS ratio** – similar long strings
74//    are resolved by a narrow diagonal band over Myers' blocked recurrence (O(n · d/w)
75//    instead of O(n · m/w)).  Band widths double geometrically until certifying the
76//    exact result, with early bailout on dissimilar pairs.
77//
78// [`ratio`] reuses steps 1–3 verbatim and swaps step 4 for the bit-parallel
79// LCS recurrence of Crochemore et al. (2001) / Hyyrö (2004), from which the
80// indel distance follows as `m + n - 2 · LLCS`; indel distance is invariant
81// under common-affix stripping for the same reason Levenshtein distance is.
82// Step 5 has an LCS counterpart that restricts the recurrence to a band of
83// diagonals and certifies its result (see `LargeCtx::lcs_banded`).  mbleven
84// stays Levenshtein-only: against the three-op LCS step it does not pay.
85//
86// The crate builds for each target's baseline ISA.  On x86-64 the whole
87// pipeline is compiled a second time with the x86-64-v2 features, and the
88// module picks one of the two copies when it loads (see [`Isa`]).
89
90use std::marker::PhantomData;
91use std::mem::MaybeUninit;
92use std::os::raw::c_uint;
93
94mod kernels;
95#[cfg(feature = "python")]
96mod python;
97
98// ---------------------------------------------------------------------------
99// Public API
100// ---------------------------------------------------------------------------
101
102/// Levenshtein edit distance between two strings.
103///
104/// The distance is the minimum number of single-character insertions,
105/// deletions, or substitutions required to transform `s1` into `s2`.
106///
107/// Characters are counted as Unicode scalar values (`char`s), so an accented
108/// letter or an emoji counts as one character.
109///
110/// # Examples
111///
112/// ```
113/// assert_eq!(lev::distance("kitten", "sitting"), 3);
114/// assert_eq!(lev::distance("flaw", "lawn"), 2);
115/// assert_eq!(lev::distance("résumé", "resume"), 2);
116/// assert_eq!(lev::distance("日本語", "日本"), 1);
117/// ```
118#[must_use]
119pub fn distance(s1: &str, s2: &str) -> usize {
120    let (u1, u2) = (Units::new(s1), Units::new(s2));
121    let (v1, v2) = (u1.view(), u2.view());
122    #[cfg(target_arch = "x86_64")]
123    if kernels::x86_64::v2::detected() {
124        // SAFETY: the CPU has the features of `V2`, and the views borrow `u1` and `u2`.
125        return unsafe { kernels::x86_64::v2::distance_views(&v1, &v2) };
126    }
127    // SAFETY: the views borrow `u1` and `u2`.
128    unsafe { compute::<Lev<Baseline>>(&v1, &v2) }
129}
130
131/// Calculate how similar two strings are, as a score from `0.0` to `1.0`.
132///
133/// A score of `1.0` means the strings are identical, and `0.0` means they
134/// have no character in common.  In between, the score is the share of
135/// characters the two strings have in common, known as the indel similarity
136/// ratio:
137///
138/// ```text
139/// ratio(s1, s2) = 2 · LCS(s1, s2) / (|s1| + |s2|)
140/// ```
141///
142/// where `LCS` is the length of the longest common subsequence and `|s|` is
143/// the number of characters in `s`.  Two empty strings score `1.0`.
144///
145/// # Examples
146///
147/// ```
148/// assert_eq!(lev::ratio("kitten", "sitting"), 0.6153846153846154);
149/// assert_eq!(lev::ratio("abc", "xyz"), 0.0);
150/// assert_eq!(lev::ratio("", ""), 1.0);
151/// assert_eq!(lev::ratio("résumé", "resume"), 0.6666666666666667);
152/// ```
153#[must_use]
154pub fn ratio(s1: &str, s2: &str) -> f64 {
155    let (u1, u2) = (Units::new(s1), Units::new(s2));
156    let (v1, v2) = (u1.view(), u2.view());
157    let total = v1.len + v2.len;
158    if total == 0 {
159        return 1.0;
160    }
161    #[cfg(target_arch = "x86_64")]
162    if kernels::x86_64::v2::detected() {
163        // SAFETY: the CPU has the features of `V2`, and the views borrow `u1` and `u2`.
164        let indel = unsafe { kernels::x86_64::v2::indel_views(&v1, &v2) };
165        return 1.0 - indel as f64 / total as f64;
166    }
167    // SAFETY: the views borrow `u1` and `u2`.
168    1.0 - unsafe { compute::<Lcs<Baseline>>(&v1, &v2) } as f64 / total as f64
169}
170
171/// A `&str` as code units in the narrowest of CPython's three string
172/// encodings that holds it, so Rust strings take the same paths through
173/// [`compute`] as Python strings with the same characters.  Pure-ASCII
174/// strings are borrowed as they are; others are transcoded once.
175enum Units<'a> {
176    Ascii(&'a [u8]),
177    Latin1(Vec<u8>),
178    Ucs2(Vec<u16>),
179    Ucs4(Vec<u32>),
180}
181
182impl<'a> Units<'a> {
183    fn new(s: &'a str) -> Self {
184        if s.is_ascii() {
185            return Units::Ascii(s.as_bytes());
186        }
187        // `collect` from `Chars` would start at a quarter of the UTF-8 length
188        // and grow; the UTF-8 length bounds the char count, so one
189        // allocation suffices.
190        fn transcode<T>(s: &str, unit: impl Fn(char) -> T) -> Vec<T> {
191            let mut units = Vec::with_capacity(s.len());
192            units.extend(s.chars().map(unit));
193            units
194        }
195        match s.chars().max().map_or(0, u32::from) {
196            0..=0xFF => Units::Latin1(transcode(s, |c| c as u8)),
197            0x100..=0xFFFF => Units::Ucs2(transcode(s, |c| c as u16)),
198            _ => Units::Ucs4(transcode(s, u32::from)),
199        }
200    }
201
202    fn view(&self) -> UniView {
203        let (kind, ascii, data, len) = match self {
204            Units::Ascii(b) => (KIND_1BYTE, true, b.as_ptr(), b.len()),
205            Units::Latin1(b) => (KIND_1BYTE, false, b.as_ptr(), b.len()),
206            Units::Ucs2(u) => (KIND_2BYTE, false, u.as_ptr().cast::<u8>(), u.len()),
207            Units::Ucs4(u) => (KIND_4BYTE, false, u.as_ptr().cast::<u8>(), u.len()),
208        };
209        UniView {
210            kind,
211            ascii,
212            data,
213            len,
214        }
215    }
216}
217
218// ---------------------------------------------------------------------------
219// Code units of the three CPython string kinds (UCS-1 / UCS-2 / UCS-4).
220// ---------------------------------------------------------------------------
221
222/// A CPython code unit (`u8`, `u16` or `u32`).  Mixed-kind pairs compare
223/// code units widened to `u64`, or to `u32` where blocks of them are
224/// compared as vectors.
225trait CodeUnit: Copy + Eq + Into<u64> + Into<u32> {
226    #[inline(always)]
227    fn as_u64(self) -> u64 {
228        Into::<u64>::into(self)
229    }
230
231    #[inline(always)]
232    fn as_u32(self) -> u32 {
233        Into::<u32>::into(self)
234    }
235}
236impl CodeUnit for u8 {}
237impl CodeUnit for u16 {}
238impl CodeUnit for u32 {}
239
240/// Fibonacci hash slot: maps a 64-bit key into `[0, mask]` where `mask = 2^k - 1`.
241///
242/// Shifts by `64 - k` to extract the top k bits of the product.  The high bits
243/// of a multiplicative hash have the best avalanche properties — they mix
244/// contributions from all input bits via carry propagation — so this gives
245/// better distribution than extracting middle or low bits.  The result is
246/// already in `[0, mask]`, so no masking step is needed.
247///
248/// See <https://en.wikipedia.org/wiki/Hash_function#Fibonacci_hashing>.
249#[inline(always)]
250fn hslot(key: u64, shift: u32) -> usize {
251    (key.wrapping_mul(0x9e37_79b9_7f4a_7c15_u64) >> shift) as usize
252}
253
254// ---------------------------------------------------------------------------
255// CPython internal-buffer accessors
256// ---------------------------------------------------------------------------
257
258/// [`UniView::kind`] of UCS-1 (`u8`) code units, as in CPython's `PyUnicode_Kind`.
259const KIND_1BYTE: c_uint = 1;
260/// [`UniView::kind`] of UCS-2 (`u16`) code units.
261const KIND_2BYTE: c_uint = 2;
262/// [`UniView::kind`] of UCS-4 (`u32`) code units.
263const KIND_4BYTE: c_uint = 4;
264
265/// Snapshot of a string's code units, laid out like a Python string's
266/// internals.  For a Python string all four fields are derived from a
267/// single PyASCIIObject header read; subsequent dispatch branches operate
268/// on locals only — no further FFI calls.
269struct UniView {
270    kind: c_uint,
271    ascii: bool,
272    data: *const u8,
273    len: usize,
274}
275
276/// Reinterpret a [`UniView`]'s data as a `u8` slice.
277///
278/// # Safety
279///
280/// `v.data` must point to `v.len` valid, initialized `u8` code units.
281#[inline(always)]
282unsafe fn as_u8(v: &UniView) -> &[u8] {
283    // SAFETY: caller guarantees v.data points to v.len valid u8 elements.
284    unsafe { std::slice::from_raw_parts(v.data, v.len) }
285}
286
287/// Reinterpret a [`UniView`]'s data as a `u16` slice.
288///
289/// # Safety
290///
291/// `v.data` must be properly aligned for `u16` and point to `v.len` valid, initialized `u16` code units.
292#[inline(always)]
293unsafe fn as_u16(v: &UniView) -> &[u16] {
294    // SAFETY: caller guarantees v.data is aligned and points to v.len valid u16 elements.
295    unsafe { std::slice::from_raw_parts(v.data.cast::<u16>(), v.len) }
296}
297
298/// Reinterpret a [`UniView`]'s data as a `u32` slice.
299///
300/// # Safety
301///
302/// `v.data` must be properly aligned for `u32` and point to `v.len` valid, initialized `u32` code units.
303#[inline(always)]
304unsafe fn as_u32(v: &UniView) -> &[u32] {
305    // SAFETY: caller guarantees v.data is aligned and points to v.len valid u32 elements.
306    unsafe { std::slice::from_raw_parts(v.data.cast::<u32>(), v.len) }
307}
308
309// ---------------------------------------------------------------------------
310// Copies of the hot path, one per instruction set
311// ---------------------------------------------------------------------------
312
313/// One copy of the hot path, compiled for one instruction set, which
314/// [`Kernel::Isa`] carries through the pipeline.
315///
316/// [`Baseline`] is the copy for the target's baseline ISA, which every CPU
317/// of the target runs, and the only one off x86-64.  On x86-64 the crate
318/// builds for the x86-64-v1 ISA, and [`kernels::x86_64::v2::V2`] adds the
319/// x86-64-v2 features; [`_lev`] installs its entry points when the CPU has
320/// them.  A crate-wide `-C target-cpu=x86-64-v2` would get the features
321/// too, but lets LLVM use them anywhere, so the wheels would crash with
322/// SIGILL on older CPUs.  Detecting them per kernel instead would leave the
323/// code between the kernels (affix scans, peq tables, the final popcounts)
324/// on x86-64-v1.
325///
326/// Code inlined into a function is compiled with that function's target
327/// features, so a copy needs its own form of each function the pipeline
328/// keeps out of line.  [`hyrro_64_tiny`], [`hyrro_multiword_direct`],
329/// [`hyrro_multiword_hash`] and [`LargeCtx::banded`] call their copy's form
330/// of their body; [`mbleven`] is called through [`mbleven_in`], and the
331/// copies with target features inline [`LargeCtx::run`] by hand.  These
332/// functions test [`Isa::X86_64_V2`], a constant, rather than call methods
333/// of the copy: rustc inlines MIR only into calls it can resolve before
334/// monomorphization, and calls through trait methods changed the AArch64
335/// code of the baseline copy.  As it is, that code is the same per function
336/// as before there were copies.
337///
338/// A pipeline of a copy other than [`Baseline`] is only reached from that
339/// copy's entry points, so the CPU has its features.
340trait Isa {
341    /// Whether the copy has the x86-64-v2 features: its out-of-line
342    /// functions are those of [`kernels::x86_64::v2`], and the forms of
343    /// [`kernels::x86_64::bmi1`], which need them besides BMI1, may run.
344    #[cfg_attr(not(target_arch = "x86_64"), allow(dead_code))]
345    const X86_64_V2: bool;
346}
347
348/// The copy of the hot path for the target's baseline ISA (see [`Isa`]).
349struct Baseline;
350
351impl Isa for Baseline {
352    const X86_64_V2: bool = false;
353}
354
355// ---------------------------------------------------------------------------
356// Top-level dispatch: read both views' kinds, branch once, run the right path.
357// ---------------------------------------------------------------------------
358
359/// Dispatch on `(kind, ascii)` tuples and return the `K` edit distance directly.
360///
361/// # Safety
362///
363/// Both views must point to valid initialized code-unit buffers corresponding
364/// to their respective `kind` flags.
365#[inline(always)]
366unsafe fn compute<K: Kernel>(v1: &UniView, v2: &UniView) -> usize {
367    use self::{KIND_1BYTE as K1, KIND_2BYTE as K2, KIND_4BYTE as K4};
368
369    // 1. Normalize: v1 is the shorter pattern string (m), v2 is the text (n).
370    let (v1, v2) = if v1.len <= v2.len { (v1, v2) } else { (v2, v1) };
371
372    if v1.len == 0 {
373        return v2.len;
374    }
375
376    // 2. Dispatch based on kinds.
377    match (v1.kind, v2.kind) {
378        (K1, K1) => {
379            // SAFETY: kind is 1-byte, so data is valid &[u8].
380            let (b1, b2) = unsafe { (as_u8(v1), as_u8(v2)) };
381            if v1.ascii && v2.ascii {
382                compute_u8::<K, true>(b1, b2)
383            } else {
384                compute_u8::<K, false>(b1, b2)
385            }
386        }
387        // SAFETY: kind is 2-byte, so data is valid &[u16].
388        (K2, K2) => compute_sorted::<K, _>(unsafe { as_u16(v1) }, unsafe { as_u16(v2) }),
389        // SAFETY: kind is 4-byte, so data is valid &[u32].
390        (K4, K4) => compute_sorted::<K, _>(unsafe { as_u32(v1) }, unsafe { as_u32(v2) }),
391        // Mixed kinds compare code units in place, without a widened copy.
392        // SAFETY: kinds are 1- and 2-byte, so data are valid &[u8] and &[u16].
393        (K1, K2) => compute_sorted_mixed::<K, _, _>(unsafe { as_u8(v1) }, unsafe { as_u16(v2) }),
394        // SAFETY: kinds are 1- and 4-byte, so data are valid &[u8] and &[u32].
395        (K1, K4) => compute_sorted_mixed::<K, _, _>(unsafe { as_u8(v1) }, unsafe { as_u32(v2) }),
396        // SAFETY: kinds are 2- and 1-byte, so data are valid &[u16] and &[u8].
397        (K2, K1) => compute_sorted_mixed::<K, _, _>(unsafe { as_u16(v1) }, unsafe { as_u8(v2) }),
398        // SAFETY: kinds are 2- and 4-byte, so data are valid &[u16] and &[u32].
399        (K2, K4) => compute_sorted_mixed::<K, _, _>(unsafe { as_u16(v1) }, unsafe { as_u32(v2) }),
400        // SAFETY: kinds are 4- and 1-byte, so data are valid &[u32] and &[u8].
401        (K4, K1) => compute_sorted_mixed::<K, _, _>(unsafe { as_u32(v1) }, unsafe { as_u8(v2) }),
402        // SAFETY: kinds are 4- and 2-byte, so data are valid &[u32] and &[u16].
403        (K4, K2) => compute_sorted_mixed::<K, _, _>(unsafe { as_u32(v1) }, unsafe { as_u16(v2) }),
404        _ => unreachable!("unknown string kind"),
405    }
406}
407
408// ---------------------------------------------------------------------------
409// Affix stripping
410// ---------------------------------------------------------------------------
411
412/// Code units compared per branch by the mixed-kind scans.
413const MIXED_BLOCK: usize = 16;
414
415/// Number of positions where two blocks of different code-unit types
416/// differ.  No early exit, so LLVM widens and compares the block as vectors.
417#[inline(always)]
418fn block_mismatches<T1: CodeUnit, T2: CodeUnit>(
419    x: &[T1; MIXED_BLOCK],
420    y: &[T2; MIXED_BLOCK],
421) -> usize {
422    x.iter()
423        .zip(y)
424        .map(|(p, q)| (p.as_u32() != q.as_u32()) as usize)
425        .sum()
426}
427
428/// Strip common prefix and suffix from two slices of different types: whole
429/// blocks first, then code unit by code unit from the first differing block.
430#[inline(always)]
431fn strip_affix_mixed<'a, T1: CodeUnit, T2: CodeUnit>(
432    a: &'a [T1],
433    b: &'a [T2],
434) -> (&'a [T1], &'a [T2]) {
435    let n = a.len().min(b.len());
436    let mut prefix = 0;
437    while prefix + MIXED_BLOCK <= n {
438        let Some(blk_a) = a[prefix..].first_chunk::<MIXED_BLOCK>() else {
439            break;
440        };
441        let Some(blk_b) = b[prefix..].first_chunk::<MIXED_BLOCK>() else {
442            break;
443        };
444        if block_mismatches(blk_a, blk_b) != 0 {
445            break;
446        }
447        prefix += MIXED_BLOCK;
448    }
449    while prefix < n && a[prefix].as_u64() == b[prefix].as_u64() {
450        prefix += 1;
451    }
452    if prefix == a.len() {
453        return (&a[..0], &b[prefix..]);
454    }
455
456    let a = &a[prefix..];
457    let b = &b[prefix..];
458
459    let (n, la, lb) = (n - prefix, a.len(), b.len());
460    let mut suffix = 0;
461    while suffix + MIXED_BLOCK <= n {
462        let Some(blk_a) = a[..la - suffix].last_chunk::<MIXED_BLOCK>() else {
463            break;
464        };
465        let Some(blk_b) = b[..lb - suffix].last_chunk::<MIXED_BLOCK>() else {
466            break;
467        };
468        if block_mismatches(blk_a, blk_b) != 0 {
469            break;
470        }
471        suffix += MIXED_BLOCK;
472    }
473    while suffix < n && a[la - suffix - 1].as_u64() == b[lb - suffix - 1].as_u64() {
474        suffix += 1;
475    }
476
477    (&a[..la - suffix], &b[..lb - suffix])
478}
479
480/// Reinterpret a code-unit slice as raw bytes (sound: u8/u16/u32 are plain data).
481#[inline(always)]
482pub(crate) fn as_bytes<T: CodeUnit>(s: &[T]) -> &[u8] {
483    // SAFETY: `s` is a valid slice of `CodeUnit` primitives (`u8`, `u16`, or `u32`),
484    // which are plain data with no uninitialized bytes or padding. The total byte
485    // length `std::mem::size_of_val(s)` correctly represents the backing memory.
486    unsafe { std::slice::from_raw_parts(s.as_ptr().cast::<u8>(), std::mem::size_of_val(s)) }
487}
488
489/// Bytes compared per branch by the affix scans: the XORs of four words are
490/// ORed into one test, which LLVM lowers to 16-byte vector compares.
491const AFFIX_BLOCK: usize = 32;
492
493/// Whether two `AFFIX_BLOCK`-byte blocks differ, without an early exit.
494#[inline(always)]
495fn blocks_differ(x: &[u8; AFFIX_BLOCK], y: &[u8; AFFIX_BLOCK]) -> bool {
496    let (wx, _) = x.as_chunks::<8>();
497    let (wy, _) = y.as_chunks::<8>();
498    wx.iter().zip(wy).fold(0u64, |acc, (p, q)| {
499        acc | (u64::from_ne_bytes(*p) ^ u64::from_ne_bytes(*q))
500    }) != 0
501}
502
503/// Length in bytes of the common prefix of `a` and `b`: whole blocks first,
504/// then 8 bytes at a time from the first differing block on.
505#[inline(always)]
506fn common_prefix_bytes(a: &[u8], b: &[u8]) -> usize {
507    let n = a.len().min(b.len());
508    let mut i = 0;
509    while i + AFFIX_BLOCK <= n {
510        let Some(blk_a) = a[i..].first_chunk::<AFFIX_BLOCK>() else {
511            break;
512        };
513        let Some(blk_b) = b[i..].first_chunk::<AFFIX_BLOCK>() else {
514            break;
515        };
516        if blocks_differ(blk_a, blk_b) {
517            break;
518        }
519        i += AFFIX_BLOCK;
520    }
521    while i + 8 <= n {
522        let Some(chunk_a) = a[i..].first_chunk::<8>() else {
523            break;
524        };
525        let Some(chunk_b) = b[i..].first_chunk::<8>() else {
526            break;
527        };
528        let diff = u64::from_le_bytes(*chunk_a) ^ u64::from_le_bytes(*chunk_b);
529        if diff != 0 {
530            return i + (diff.trailing_zeros() / 8) as usize;
531        }
532        i += 8;
533    }
534    while i < n && a[i] == b[i] {
535        i += 1;
536    }
537    i
538}
539
540/// Length in bytes of the common suffix of `a` and `b`, scanned like
541/// [`common_prefix_bytes`] from the end.
542#[inline(always)]
543fn common_suffix_bytes(a: &[u8], b: &[u8]) -> usize {
544    let n = a.len().min(b.len());
545    let (la, lb) = (a.len(), b.len());
546    let mut i = 0;
547    while i + AFFIX_BLOCK <= n {
548        let Some(blk_a) = a[..la - i].last_chunk::<AFFIX_BLOCK>() else {
549            break;
550        };
551        let Some(blk_b) = b[..lb - i].last_chunk::<AFFIX_BLOCK>() else {
552            break;
553        };
554        if blocks_differ(blk_a, blk_b) {
555            break;
556        }
557        i += AFFIX_BLOCK;
558    }
559    while i + 8 <= n {
560        let Some(chunk_a) = a[..la - i].last_chunk::<8>() else {
561            break;
562        };
563        let Some(chunk_b) = b[..lb - i].last_chunk::<8>() else {
564            break;
565        };
566        let diff = u64::from_le_bytes(*chunk_a) ^ u64::from_le_bytes(*chunk_b);
567        if diff != 0 {
568            return i + (diff.leading_zeros() / 8) as usize;
569        }
570        i += 8;
571    }
572    while i < n && a[la - i - 1] == b[lb - i - 1] {
573        i += 1;
574    }
575    i
576}
577
578/// Strip common prefix and suffix from two same-kind slices.
579///
580/// Code-unit equality is byte equality for equal-width slices, so the scan
581/// runs on the raw bytes 8 at a time; matched byte counts are floor-divided
582/// by the element size, which discards any partial element at a difference
583/// boundary.
584#[inline(always)]
585fn strip_affix<'a, T: CodeUnit>(a: &'a [T], b: &'a [T]) -> (&'a [T], &'a [T]) {
586    let size = std::mem::size_of::<T>();
587    let prefix = common_prefix_bytes(as_bytes(a), as_bytes(b)) / size;
588    if prefix == a.len() {
589        return (&a[..0], &b[prefix..]);
590    }
591
592    let a = &a[prefix..];
593    let b = &b[prefix..];
594
595    let suffix = common_suffix_bytes(as_bytes(a), as_bytes(b)) / size;
596
597    (&a[..a.len() - suffix], &b[..b.len() - suffix])
598}
599
600// ---------------------------------------------------------------------------
601// UCS-1 pipeline
602// ---------------------------------------------------------------------------
603
604/// Patterns at or below this length skip peq-table construction entirely; the
605/// pattern-match word is computed per text char by vector comparisons or a
606/// fully unrolled linear scan.  Zero-initializing the peq tables (1-2 KiB)
607/// dominates the runtime for tiny inputs, so trading one table load for ≤ `TINY_M`
608/// SIMD or compare+or ops wins.
609const TINY_M: usize = 16;
610
611/// Tiny-pattern kernel: no peq table, `pm` computed per text char directly
612/// from the pattern.
613///
614/// Single-byte patterns up to 16 characters use 128-bit SIMD comparisons
615/// (NEON on AArch64, SSE2 on x86-64), avoiding any peq table setup. Wider
616/// patterns dispatch on a const length so the scan unrolls into `m` compares.
617///
618/// Kept out of line: inlining these loops into the dispatcher perturbs
619/// register allocation of the single-word paths sharing that function.
620/// Each copy of the hot path compiles its own (see [`Isa`]).
621#[inline(never)]
622fn hyrro_64_tiny<K: Kernel, T1: CodeUnit, T2: CodeUnit>(pattern: &[T1], text: &[T2]) -> usize {
623    #[cfg(target_arch = "x86_64")]
624    if K::Isa::X86_64_V2 {
625        // SAFETY: the CPU runs `K::Isa` (see `Isa`).
626        return unsafe { kernels::x86_64::v2::hyrro_64_tiny::<K, T1, T2>(pattern, text) };
627    }
628    hyrro_64_tiny_body::<K, T1, T2>(pattern, text)
629}
630
631/// The body of [`hyrro_64_tiny`].
632#[inline(always)]
633fn hyrro_64_tiny_body<K: Kernel, T1: CodeUnit, T2: CodeUnit>(pattern: &[T1], text: &[T2]) -> usize {
634    debug_assert!((1..=TINY_M).contains(&pattern.len()));
635    if std::mem::size_of::<T1>() == 1 {
636        if pattern.len() <= 8 {
637            return tiny_swar::<K, _, _>(pattern, text);
638        }
639        #[cfg(target_arch = "aarch64")]
640        {
641            // SAFETY: AArch64 baseline guarantees NEON support, `pattern` code units
642            // are single-byte (`size_of::<T1>() == 1`), and `pattern.len() <= TINY_M`
643            // (caller contract, asserted above).
644            return unsafe { tiny_u8_simd::<K, T1, T2>(pattern, text) };
645        }
646        #[cfg(target_arch = "x86_64")]
647        {
648            // SAFETY: x86_64 baseline guarantees SSE2 support, `pattern` code units
649            // are single-byte (`size_of::<T1>() == 1`), and `pattern.len() <= TINY_M`
650            // (caller contract, asserted above).
651            return unsafe { tiny_u8_simd::<K, T1, T2>(pattern, text) };
652        }
653    }
654    match pattern.len() {
655        1 => tiny_n::<K, 1, _, _>(pattern, text),
656        2 => tiny_n::<K, 2, _, _>(pattern, text),
657        3 => tiny_n::<K, 3, _, _>(pattern, text),
658        4 => tiny_n::<K, 4, _, _>(pattern, text),
659        5 => tiny_n::<K, 5, _, _>(pattern, text),
660        6 => tiny_n::<K, 6, _, _>(pattern, text),
661        7 => tiny_n::<K, 7, _, _>(pattern, text),
662        8 => tiny_n::<K, 8, _, _>(pattern, text),
663        9 => tiny_n::<K, 9, _, _>(pattern, text),
664        10 => tiny_n::<K, 10, _, _>(pattern, text),
665        11 => tiny_n::<K, 11, _, _>(pattern, text),
666        12 => tiny_n::<K, 12, _, _>(pattern, text),
667        13 => tiny_n::<K, 13, _, _>(pattern, text),
668        14 => tiny_n::<K, 14, _, _>(pattern, text),
669        15 => tiny_n::<K, 15, _, _>(pattern, text),
670        16 => tiny_n::<K, 16, _, _>(pattern, text),
671        _ => unreachable!("tiny path requires 1..=TINY_M pattern chars"),
672    }
673}
674
675/// [`hyrro_64_tiny`] for single-byte patterns using 128-bit NEON vector comparisons.
676///
677/// # Safety
678///
679/// Requires `size_of::<T1>() == 1` and `pattern.len() <= TINY_M`.
680/// Target must support AArch64 NEON vector instructions.
681#[cfg(target_arch = "aarch64")]
682#[inline(always)]
683unsafe fn tiny_u8_simd<K: Kernel, T1: CodeUnit, T2: CodeUnit>(
684    pattern: &[T1],
685    text: &[T2],
686) -> usize {
687    use std::arch::aarch64::*;
688    let m = pattern.len();
689    debug_assert!(std::mem::size_of::<T1>() == 1 && (1..=TINY_M).contains(&m));
690    let mut pat_buf = [0u8; 16];
691    pat_buf[..m].copy_from_slice(as_bytes(pattern));
692    let pat_vec = vld1q_u8(pat_buf.as_ptr());
693    let valid = if m == 64 { u64::MAX } else { (1u64 << m) - 1 };
694
695    const IDX_BYTES: [u8; 16] = [0, 8, 1, 9, 2, 10, 3, 11, 4, 12, 5, 13, 6, 14, 7, 15];
696    const MBITS_U16: [u16; 8] = [
697        0x0101, 0x0202, 0x0404, 0x0808, 0x1010, 0x2020, 0x4040, 0x8080,
698    ];
699    let idx = vld1q_u8(IDX_BYTES.as_ptr());
700    let mbits = vld1q_u16(MBITS_U16.as_ptr());
701
702    K::word(
703        m,
704        text.len(),
705        text.iter().map(|&c| {
706            let c_val = c.as_u64();
707            if c_val > 0xff {
708                return 0;
709            }
710            let c_vec = vdupq_n_u8(c_val as u8);
711            let eq = vceqq_u8(pat_vec, c_vec);
712            let perm = vqtbl1q_u8(eq, idx);
713            let anded = vandq_u16(vreinterpretq_u16_u8(perm), mbits);
714            let mask = vaddvq_u16(anded) as u64;
715            mask & valid
716        }),
717    )
718}
719
720/// [`hyrro_64_tiny`] for single-byte patterns using 128-bit SSE2 vector comparisons.
721///
722/// # Safety
723///
724/// Requires `size_of::<T1>() == 1` and `pattern.len() <= TINY_M`.
725/// Target must support x86-64 SSE2 vector instructions.
726#[cfg(target_arch = "x86_64")]
727#[inline(always)]
728unsafe fn tiny_u8_simd<K: Kernel, T1: CodeUnit, T2: CodeUnit>(
729    pattern: &[T1],
730    text: &[T2],
731) -> usize {
732    use std::arch::x86_64::*;
733    let m = pattern.len();
734    debug_assert!(std::mem::size_of::<T1>() == 1 && (1..=TINY_M).contains(&m));
735    let mut pat_buf = [0u8; 16];
736    pat_buf[..m].copy_from_slice(as_bytes(pattern));
737    let pat_vec = _mm_loadu_si128(pat_buf.as_ptr() as *const __m128i);
738    let valid = if m == 64 { u64::MAX } else { (1u64 << m) - 1 };
739
740    K::word(
741        m,
742        text.len(),
743        text.iter().map(|&c| {
744            let c_val = c.as_u64();
745            if c_val > 0xff {
746                return 0;
747            }
748            let c_vec = _mm_set1_epi8(c_val as i8);
749            let eq = _mm_cmpeq_epi8(pat_vec, c_vec);
750            let mask = _mm_movemask_epi8(eq) as u64;
751            mask & valid
752        }),
753    )
754}
755
756/// [`hyrro_64_tiny`] for a multi-byte pattern of exactly `N` code units.
757#[inline(always)]
758fn tiny_n<K: Kernel, const N: usize, T1: CodeUnit, T2: CodeUnit>(
759    pattern: &[T1],
760    text: &[T2],
761) -> usize {
762    let pattern: [u64; N] = std::array::from_fn(|i| pattern[i].as_u64());
763    K::word(
764        N,
765        text.len(),
766        text.iter().map(|&c| {
767            let c = c.as_u64();
768            let mut pm = 0u64;
769            for (i, &p) in pattern.iter().enumerate() {
770                pm |= ((p == c) as u64) << i;
771            }
772            pm
773        }),
774    )
775}
776
777/// [`hyrro_64_tiny`] for single-byte patterns: the pattern is packed one byte
778/// per lane into a u64 and matched against each text char with SWAR, ~10
779/// scalar ops per char.
780#[inline(always)]
781fn tiny_swar<K: Kernel, T1: CodeUnit, T2: CodeUnit>(pattern: &[T1], text: &[T2]) -> usize {
782    const LO7: u64 = 0x7f7f_7f7f_7f7f_7f7f;
783    const HI: u64 = 0x8080_8080_8080_8080;
784    const ONES: u64 = 0x0101_0101_0101_0101;
785    // Multiplying the lane LSBs (bits 8i) by this moves bit 8i to 56 + i with
786    // no overlapping partial products, i.e. no carries.
787    const GATHER: u64 = 0x0102_0408_1020_4080;
788    let m = pattern.len();
789    debug_assert!(std::mem::size_of::<T1>() == 1 && (1..=8).contains(&m));
790
791    let packed = pattern
792        .iter()
793        .rev()
794        .fold(0u64, |acc, &p| (acc << 8) | p.as_u64());
795    // Lanes >= m hold zero bytes and would match NUL; mask them off.
796    let valid = u64::MAX >> (64 - m);
797    K::word(
798        m,
799        text.len(),
800        text.iter().map(|&c| {
801            let c = c.as_u64();
802            if c > 0xff {
803                return 0;
804            }
805            // Exact zero-byte test: the high bit of each lane of `hit` is set
806            // iff that lane of `x` is zero (no cross-lane borrows).
807            let x = packed ^ c.wrapping_mul(ONES);
808            let hit = !(((x & LO7) + LO7) | x) & HI;
809            ((hit >> 7).wrapping_mul(GATHER) >> 56) & valid
810        }),
811    )
812}
813
814/// Strip affixes and run kernel `K`.  `ASCII = true` enables the 128-entry peq
815/// fast path.
816#[inline(always)]
817fn compute_u8<K: Kernel, const ASCII: bool>(a: &[u8], b: &[u8]) -> usize {
818    let (a, b) = strip_affix(a, b);
819    if a.is_empty() {
820        return b.len();
821    }
822    if a.len() <= TINY_M {
823        return hyrro_64_tiny::<K, _, _>(a, b);
824    }
825    if K::MBLEVEN {
826        if let Some(ub) = small_ub::<K::Isa, _>(a, b) {
827            return mbleven_in::<K::Isa, _, _>(a, b, ub);
828        }
829    }
830    if a.len() <= 64 {
831        if ASCII {
832            hyrro_64_u8::<K, 128>(a, b)
833        } else {
834            hyrro_64_u8::<K, 256>(a, b)
835        }
836    } else {
837        if ASCII {
838            hyrro_multiword_direct::<K, 128, _, _>(a, b)
839        } else {
840            hyrro_multiword_direct::<K, 256, _, _>(a, b)
841        }
842    }
843}
844
845// ---------------------------------------------------------------------------
846// Small-distance fast path (mbleven)
847// ---------------------------------------------------------------------------
848
849/// Largest internally-derived upper bound handled by [`mbleven`].
850pub(crate) const MBLEVEN_MAX: usize = 3;
851
852/// Longest pattern eligible for the mbleven gate.  Bounds the gate's
853/// hamming scan; beyond this the banded kernels already handle
854/// near-identical strings well.
855pub(crate) const MBLEVEN_MAX_LEN: usize = 512;
856
857/// Gate for the mbleven fast path (same-kind slices): `Some(ub)` when
858/// `length difference + aligned mismatches <= MBLEVEN_MAX`.
859///
860/// Adapted from SIMD string equality & mismatch scanning techniques (such as
861/// simdjson / rapidfuzz vector classification): compares chunks in 32-byte
862/// AVX2 blocks on x86-64 when available, or 16-byte SIMD vector blocks on targets
863/// supporting NEON or SSE2, falling back to 8-byte scalar blocks. Only blocks
864/// containing differences inspect the mismatch count, and aborts as soon as
865/// `ub > MBLEVEN_MAX`.
866#[inline(always)]
867fn small_ub<X: Isa, T: CodeUnit>(short: &[T], long: &[T]) -> Option<usize> {
868    #[cfg(target_arch = "x86_64")]
869    if X::X86_64_V2 {
870        // SAFETY: the CPU runs `X` (see `Isa`).
871        return unsafe { kernels::x86_64::v2::small_ub::<T>(short, long) };
872    }
873    small_ub_body::<T>(short, long)
874}
875
876/// The body of [`small_ub`] for baseline ISAs.
877#[inline(always)]
878pub(crate) fn small_ub_body<T: CodeUnit>(short: &[T], long: &[T]) -> Option<usize> {
879    debug_assert!(short.len() <= long.len());
880    let mut ub = long.len() - short.len();
881    if ub > MBLEVEN_MAX || short.len() > MBLEVEN_MAX_LEN {
882        return None;
883    }
884    let size = std::mem::size_of::<T>();
885    let a = as_bytes(short);
886    let b = &as_bytes(long)[..a.len()];
887    let n = a.len();
888    let mut i = 0;
889
890    #[cfg(target_arch = "aarch64")]
891    // SAFETY: `i + 16 <= n`, so reading 16 bytes starting at `a.as_ptr().add(i)` and
892    // `b.as_ptr().add(i)` is within the bounds of `a` and `b`. NEON intrinsics are part
893    // of the baseline AArch64 ISA.
894    unsafe {
895        use std::arch::aarch64::*;
896        while i + 16 <= n {
897            let va = vld1q_u8(a.as_ptr().add(i));
898            let vb = vld1q_u8(b.as_ptr().add(i));
899            let eq = vceqq_u8(va, vb);
900            let neq = vmvnq_u8(eq);
901            if vmaxvq_u8(neq) != 0 {
902                let diff_lanes = match size {
903                    1 => vaddlvq_u8(vshrq_n_u8(neq, 7)) as usize,
904                    2 => {
905                        let neq16 = vmvnq_u16(vceqq_u16(
906                            vreinterpretq_u16_u8(va),
907                            vreinterpretq_u16_u8(vb),
908                        ));
909                        vaddlvq_u16(vshrq_n_u16(neq16, 15)) as usize
910                    }
911                    4 => {
912                        let neq32 = vmvnq_u32(vceqq_u32(
913                            vreinterpretq_u32_u8(va),
914                            vreinterpretq_u32_u8(vb),
915                        ));
916                        vaddlvq_u32(vshrq_n_u32(neq32, 31)) as usize
917                    }
918                    _ => unreachable!(),
919                };
920                ub += diff_lanes;
921                if ub > MBLEVEN_MAX {
922                    return None;
923                }
924            }
925            i += 16;
926        }
927    }
928
929    #[cfg(target_arch = "x86_64")]
930    // SAFETY: `i + 16 <= n`, so unaligned 16-byte reads from `a.as_ptr().add(i)` and
931    // `b.as_ptr().add(i)` are within bounds. SSE2 is part of the baseline x86-64 ISA.
932    unsafe {
933        use std::arch::x86_64::*;
934        while i + 16 <= n {
935            let va = _mm_loadu_si128(a.as_ptr().add(i) as *const __m128i);
936            let vb = _mm_loadu_si128(b.as_ptr().add(i) as *const __m128i);
937            match size {
938                1 => {
939                    let eq = _mm_cmpeq_epi8(va, vb);
940                    let mask = _mm_movemask_epi8(eq) as u32;
941                    let neq_mask = (!mask) & 0xFFFF;
942                    if neq_mask != 0 {
943                        ub += neq_mask.count_ones() as usize;
944                        if ub > MBLEVEN_MAX {
945                            return None;
946                        }
947                    }
948                }
949                2 => {
950                    let eq = _mm_cmpeq_epi16(va, vb);
951                    let mask = _mm_movemask_epi8(eq) as u32;
952                    let neq_mask = (!mask) & 0xFFFF;
953                    if neq_mask != 0 {
954                        ub += (neq_mask.count_ones() / 2) as usize;
955                        if ub > MBLEVEN_MAX {
956                            return None;
957                        }
958                    }
959                }
960                4 => {
961                    let eq = _mm_cmpeq_epi32(va, vb);
962                    let mask = _mm_movemask_epi8(eq) as u32;
963                    let neq_mask = (!mask) & 0xFFFF;
964                    if neq_mask != 0 {
965                        ub += (neq_mask.count_ones() / 4) as usize;
966                        if ub > MBLEVEN_MAX {
967                            return None;
968                        }
969                    }
970                }
971                _ => unreachable!(),
972            }
973            i += 16;
974        }
975    }
976
977    let lane_bits = 8 * size;
978    let lane_mask = if lane_bits >= 64 {
979        u64::MAX
980    } else {
981        (1u64 << lane_bits) - 1
982    };
983
984    while i + 8 <= n {
985        let Some(ca) = a[i..].first_chunk::<8>() else {
986            break;
987        };
988        let Some(cb) = b[i..].first_chunk::<8>() else {
989            break;
990        };
991        let mut d = u64::from_le_bytes(*ca) ^ u64::from_le_bytes(*cb);
992        while d != 0 {
993            ub += 1;
994            if ub > MBLEVEN_MAX {
995                return None;
996            }
997            let lane = d.trailing_zeros() as usize / lane_bits;
998            d &= !(lane_mask << (lane * lane_bits));
999        }
1000        i += 8;
1001    }
1002    for e in i / size..short.len() {
1003        if short[e].as_u64() != long[e].as_u64() {
1004            ub += 1;
1005            if ub > MBLEVEN_MAX {
1006                return None;
1007            }
1008        }
1009    }
1010    Some(ub)
1011}
1012
1013/// [`small_ub`] for mixed-kind slices (no byte-comparable representation):
1014/// mismatches are counted a block at a time, aborting between blocks.
1015#[inline(always)]
1016fn small_ub_mixed<T1: CodeUnit, T2: CodeUnit>(short: &[T1], long: &[T2]) -> Option<usize> {
1017    debug_assert!(short.len() <= long.len());
1018    let mut ub = long.len() - short.len();
1019    if ub > MBLEVEN_MAX || short.len() > MBLEVEN_MAX_LEN {
1020        return None;
1021    }
1022    let blocks = short.len() / MIXED_BLOCK * MIXED_BLOCK;
1023    for i in (0..blocks).step_by(MIXED_BLOCK) {
1024        // `i + MIXED_BLOCK <= blocks <= short.len() <= long.len()`, so both chunks exist;
1025        // a silent `break` here would skip blocks and undercount `ub`.
1026        let Some(blk_s) = short[i..].first_chunk::<MIXED_BLOCK>() else {
1027            unreachable!("block lies within short");
1028        };
1029        let Some(blk_l) = long[i..].first_chunk::<MIXED_BLOCK>() else {
1030            unreachable!("block lies within long");
1031        };
1032        ub += block_mismatches(blk_s, blk_l);
1033        if ub > MBLEVEN_MAX {
1034            return None;
1035        }
1036    }
1037    for (x, y) in short[blocks..].iter().zip(long[blocks..].iter()) {
1038        if x.as_u64() != y.as_u64() {
1039            ub += 1;
1040            if ub > MBLEVEN_MAX {
1041                return None;
1042            }
1043        }
1044    }
1045    Some(ub)
1046}
1047
1048/// [`mbleven`] in the copy `X` of the hot path: the function itself, or on
1049/// x86-64 its form in [`kernels::x86_64::v2`].  Unlike the other functions
1050/// that stay out of line, `mbleven` is left to LLVM to inline or not: an
1051/// `#[inline(never)]` wrapper that picks the copy changed its AArch64 code.
1052#[inline(always)]
1053fn mbleven_in<X: Isa, T1: CodeUnit, T2: CodeUnit>(short: &[T1], long: &[T2], ub: usize) -> usize {
1054    #[cfg(target_arch = "x86_64")]
1055    if X::X86_64_V2 {
1056        // SAFETY: the CPU runs `X` (see `Isa`).
1057        return unsafe { kernels::x86_64::v2::mbleven::<X, T1, T2>(short, long, ub) };
1058    }
1059    mbleven::<X, T1, T2>(short, long, ub)
1060}
1061
1062/// Common prefix of two code-unit slices.  Same-kind slices compare words via
1063/// [`common_prefix_bytes`]; mixed-kind slices compare code units linearly.
1064#[inline(always)]
1065fn common_prefix<T1: CodeUnit, T2: CodeUnit>(a: &[T1], b: &[T2]) -> usize {
1066    if std::mem::size_of::<T1>() == std::mem::size_of::<T2>() {
1067        let size = std::mem::size_of::<T1>();
1068        common_prefix_bytes(as_bytes(a), as_bytes(b)) / size
1069    } else {
1070        let n = a.len().min(b.len());
1071        let mut i = 0;
1072        while i < n && a[i].as_u64() == b[i].as_u64() {
1073            i += 1;
1074        }
1075        i
1076    }
1077}
1078
1079/// Exact distance when an upper bound `ub` in `1..=3` is known (mbleven,
1080/// Hyyrö 2018 formulation as used by rapidfuzz): try every candidate edit-op
1081/// sequence of cost <= `ub` and take the cheapest that aligns the strings.
1082///
1083/// Each sequence packs ops two bits per op from the LSB: `0b01` advances only
1084/// `long` (deletion from the longer string), `0b10` advances only `short`
1085/// (insertion), `0b11` advances both (substitution).  Identical matching runs
1086/// are jumped via [`common_prefix`].  O(ub * n) worst case, no peq table, no
1087/// allocation.  Instantiated per copy `X` of the hot path and called through
1088/// [`mbleven_in`].
1089fn mbleven<X: Isa, T1: CodeUnit, T2: CodeUnit>(short: &[T1], long: &[T2], ub: usize) -> usize {
1090    if ub == 0 {
1091        return 0;
1092    }
1093    debug_assert!((1..=MBLEVEN_MAX).contains(&ub));
1094    let diff = long.len() - short.len();
1095    debug_assert!(diff <= ub);
1096
1097    // Rows indexed by (ub, diff): (1,0),(1,1),(2,0),(2,1),(2,2),(3,0)..(3,3).
1098    static SEQS: [&[u64]; 9] = [
1099        &[0x03],
1100        &[0x01],
1101        &[0x0F, 0x09, 0x06],
1102        &[0x0D, 0x07],
1103        &[0x05],
1104        &[0x3F, 0x27, 0x2D, 0x39, 0x36, 0x1E, 0x1B],
1105        &[0x3D, 0x37, 0x1F, 0x25, 0x19, 0x16],
1106        &[0x35, 0x1D, 0x17],
1107        &[0x15],
1108    ];
1109    let row = (ub - 1) * (ub + 2) / 2 + diff;
1110
1111    let mut best = usize::MAX;
1112    for &seq in SEQS[row] {
1113        let mut ops = seq;
1114        let (mut sp, mut lp, mut cost) = (0usize, 0usize, 0usize);
1115        loop {
1116            let matched = common_prefix(&short[sp..], &long[lp..]);
1117            sp += matched;
1118            lp += matched;
1119            if sp == short.len() || lp == long.len() {
1120                break;
1121            }
1122            cost += 1;
1123            if ops == 0 {
1124                break;
1125            }
1126            lp += (ops & 1) as usize;
1127            sp += ((ops >> 1) & 1) as usize;
1128            ops >>= 2;
1129        }
1130        cost += (short.len() - sp) + (long.len() - lp);
1131        best = best.min(cost);
1132        if best == diff {
1133            return diff;
1134        }
1135    }
1136    best
1137}
1138
1139// ---------------------------------------------------------------------------
1140// Generic pipelines
1141// ---------------------------------------------------------------------------
1142
1143/// Strip affixes and run kernel `K` for same-kind strings (non-UCS-1).
1144#[inline(always)]
1145fn compute_sorted<K: Kernel, T: CodeUnit>(a: &[T], b: &[T]) -> usize {
1146    let (a, b) = strip_affix(a, b);
1147    if a.is_empty() {
1148        return b.len();
1149    }
1150    if a.len() <= TINY_M {
1151        return hyrro_64_tiny::<K, _, _>(a, b);
1152    }
1153    if K::MBLEVEN {
1154        if let Some(ub) = small_ub::<K::Isa, _>(a, b) {
1155            return mbleven_in::<K::Isa, _, _>(a, b, ub);
1156        }
1157    }
1158    if a.len() <= 64 {
1159        hyrro_64_hash::<K, _, _>(a, b)
1160    } else {
1161        hyrro_multiword_hash::<K, _, _>(a, b)
1162    }
1163}
1164
1165/// Strip affixes and run kernel `K` for mixed-kind strings.  A pair with a
1166/// UCS-1 side can only match below U+0100, so multi-word patterns take the
1167/// direct-table kernel instead of hashing every text char; single-word ones
1168/// keep the hash, whose 1 KiB table is cheaper to set up than 257 rows.
1169#[inline(always)]
1170fn compute_sorted_mixed<K: Kernel, T1: CodeUnit, T2: CodeUnit>(a: &[T1], b: &[T2]) -> usize {
1171    let (a, b) = strip_affix_mixed(a, b);
1172    if a.is_empty() {
1173        return b.len();
1174    }
1175    if a.len() <= TINY_M {
1176        return hyrro_64_tiny::<K, _, _>(a, b);
1177    }
1178    if K::MBLEVEN {
1179        if let Some(ub) = small_ub_mixed(a, b) {
1180            return mbleven_in::<K::Isa, _, _>(a, b, ub);
1181        }
1182    }
1183    if a.len() <= 64 {
1184        hyrro_64_hash::<K, _, _>(a, b)
1185    } else if std::mem::size_of::<T1>() == 1 || std::mem::size_of::<T2>() == 1 {
1186        hyrro_multiword_direct::<K, 257, _, _>(a, b)
1187    } else {
1188        hyrro_multiword_hash::<K, _, _>(a, b)
1189    }
1190}
1191
1192// ---------------------------------------------------------------------------
1193// Hyyrö single-word variants
1194// ---------------------------------------------------------------------------
1195
1196/// Row of code unit `c` in the direct peq table of [`hyrro_multiword_direct`].  UCS-1 units
1197/// index it as they are (the caller guarantees they are below `ROWS`).
1198/// Wider units occur only on the wide side of a mixed-kind pair with a
1199/// UCS-1 side, where `ROWS` is 257: units above U+00FF can never match
1200/// there, so they share row 256, which the text side never fills.
1201#[inline(always)]
1202fn direct_row<T: CodeUnit, const ROWS: usize>(c: T) -> usize {
1203    if std::mem::size_of::<T>() == 1 {
1204        c.as_u64() as usize
1205    } else {
1206        c.as_u64().min(ROWS as u64 - 1) as usize
1207    }
1208}
1209
1210/// Single-word variant for byte patterns (UCS-1).
1211#[inline(always)]
1212fn hyrro_64_u8<K: Kernel, const SLOTS: usize>(pattern: &[u8], text: &[u8]) -> usize {
1213    debug_assert!((1..=64).contains(&pattern.len()));
1214    let mut peq = [0u64; SLOTS];
1215    for (i, &c) in pattern.iter().enumerate() {
1216        // SAFETY: caller guarantees c < SLOTS (128 or 256).
1217        unsafe {
1218            *peq.get_unchecked_mut(c as usize) |= 1u64 << i;
1219        }
1220    }
1221    K::word(
1222        pattern.len(),
1223        text.len(),
1224        text.iter().map(|&c| {
1225            // SAFETY: When SLOTS == 128, input is verified ASCII (c < 128);
1226            // when SLOTS == 256, any u8 is < 256.
1227            unsafe { *peq.get_unchecked(c as usize) }
1228        }),
1229    )
1230}
1231
1232/// Single-word variant for UCS-2/4 and mixed kinds, with the peq table in a
1233/// stack-allocated open-addressing hash of `SLOTS` entries.
1234#[inline(always)]
1235fn hyrro_64_hash_slots<K: Kernel, const SLOTS: usize, T1: CodeUnit, T2: CodeUnit>(
1236    pattern: &[T1],
1237    text: &[T2],
1238) -> usize {
1239    let m = pattern.len();
1240    // At most half the slots fill, so every probe chain ends at an empty slot
1241    // within a few steps; a full table would make a lookup of an absent key
1242    // loop forever.
1243    debug_assert!((1..=64).contains(&m) && 2 * m <= SLOTS);
1244    const { assert!(SLOTS.is_power_of_two()) };
1245
1246    let mask = SLOTS - 1;
1247    let shift = 64 - mask.count_ones();
1248    // A slot's key is written before its value becomes non-zero and read only
1249    // after, so `vals` alone marks occupancy and `keys` needs no fill.
1250    let mut keys = [const { MaybeUninit::<T1>::uninit() }; SLOTS];
1251    let mut vals = [0u64; SLOTS];
1252
1253    for (i, &c) in pattern.iter().enumerate() {
1254        let mut slot = hslot(c.as_u64(), shift);
1255        loop {
1256            if vals[slot] == 0 {
1257                keys[slot].write(c);
1258                vals[slot] = 1u64 << i;
1259                break;
1260            }
1261            // SAFETY: vals[slot] != 0, so keys[slot] was written.
1262            if unsafe { keys[slot].assume_init() } == c {
1263                vals[slot] |= 1u64 << i;
1264                break;
1265            }
1266            slot = (slot + 1) & mask;
1267        }
1268    }
1269
1270    K::word(
1271        m,
1272        text.len(),
1273        text.iter().map(|&c| {
1274            let c = c.as_u64();
1275            let mut slot = hslot(c, shift);
1276            loop {
1277                // SAFETY: `slot < SLOTS`: `hslot` with `shift = 64 - log2(SLOTS)` keeps
1278                // only the top `log2(SLOTS)` bits, and each probe step is masked with `MASK`.
1279                let v = unsafe { *vals.get_unchecked(slot) };
1280                if v == 0 {
1281                    return 0;
1282                }
1283                // SAFETY: `v != 0`, meaning `keys[slot]` was initialized during table setup,
1284                // and `slot < SLOTS` as established above.
1285                if unsafe { keys.get_unchecked(slot).assume_init().as_u64() } == c {
1286                    return v;
1287                }
1288                slot = (slot + 1) & mask;
1289            }
1290        }),
1291    )
1292}
1293
1294/// Single-word variant for UCS-2/4 and mixed kinds, with the peq table in a
1295/// stack-allocated open-addressing hash. Adapts table size between 64 and 128
1296/// slots based on pattern length to cut stack zeroing overhead when `m <= 32`.
1297#[inline(always)]
1298fn hyrro_64_hash<K: Kernel, T1: CodeUnit, T2: CodeUnit>(pattern: &[T1], text: &[T2]) -> usize {
1299    let m = pattern.len();
1300    debug_assert!((1..=64).contains(&m));
1301
1302    if m <= 32 {
1303        hyrro_64_hash_slots::<K, 64, T1, T2>(pattern, text)
1304    } else {
1305        hyrro_64_hash_slots::<K, 128, T1, T2>(pattern, text)
1306    }
1307}
1308
1309// ---------------------------------------------------------------------------
1310// Hyyrö's bit-parallel inner loop
1311// ---------------------------------------------------------------------------
1312
1313/// Core bit-parallel loop for a pattern of `m <= 64` rows (Myers 1999 in
1314/// Hyyrö's global form: the top boundary shifts in `+1`).
1315///
1316/// Bits above row `m - 1` are never masked off: every step is bitwise, a left
1317/// shift, or an addition, so information only flows towards higher bits and
1318/// garbage above the pattern never reaches the pattern rows.
1319///
1320/// The step is arranged for a short loop-carried chain: `vp & vn == 0` keeps
1321/// `vn` out of the addition, and `(sum ^ vp) | vp == sum | vp` lets `ph` skip
1322/// Myers' `xh`.  No score is tracked per column: `vp`/`vn` hold the vertical
1323/// deltas of the current column and `D(0, j) = j`, so the distance after the
1324/// last of the `n` columns is `n + popcount(vp) - popcount(vn)` over the
1325/// pattern rows.
1326#[inline(always)]
1327fn hyrro_inner<I: Iterator<Item = u64>>(m: usize, n: usize, pm_iter: I) -> usize {
1328    debug_assert!((1..=64).contains(&m));
1329    let (mut vp, mut vn) = (!0u64, 0u64);
1330    for eq in pm_iter {
1331        (vp, vn) = kernels::myers_step(eq, vp, vn);
1332    }
1333    let rows = u64::MAX >> (64 - m);
1334    n + (vp & rows).count_ones() as usize - (vn & rows).count_ones() as usize
1335}
1336
1337// ---------------------------------------------------------------------------
1338// Bit-parallel LCS (Crochemore et al. 2001 / Hyyrö 2004)
1339// ---------------------------------------------------------------------------
1340//
1341// The state `s` is one bit per pattern row: a *zero* at row `i` marks a column
1342// position where the LCS length increases, so `LLCS = popcount(!s)` over the
1343// pattern's bits.  Starting from all-ones (column 0 has no increments) the
1344// column step is
1345//
1346//     u = s & pm;  s = (s + u) | (s - u)
1347//
1348// where `s - u == s & !pm` because `u` is a submask of `s`, so the subtraction
1349// never borrows.  Bits above the pattern length are never set in `pm`, so the
1350// `| (s - u)` term always restores them to one; a carry from the addition can
1351// disturb them only transiently and they never contribute to the popcount.
1352// Over several words the addition is one big-integer addition, so the carry
1353// is threaded from the low word upwards; `s - u` stays word-local.
1354
1355/// One column of the LCS recurrence over the words `s`, with `pm(k)` giving
1356/// the pattern-match bits of word `k`; returns the outgoing carry.
1357#[inline(always)]
1358fn lcs_step(s: &mut [u64], mut carry: bool, pm: impl Fn(usize) -> u64) -> bool {
1359    for (k, sk) in s.iter_mut().enumerate() {
1360        let u = *sk & pm(k);
1361        let (sum, nc) = sk.carrying_add(u, carry);
1362        carry = nc;
1363        *sk = sum | (*sk - u);
1364    }
1365    carry
1366}
1367
1368/// LCS length encoded in the state words `s`.
1369#[inline(always)]
1370fn lcs_len(s: &[u64]) -> usize {
1371    s.iter().map(|&x| (!x).count_ones() as usize).sum()
1372}
1373
1374/// Single-word LCS loop; returns the LCS length.
1375#[inline(always)]
1376fn lcs_inner_64<I: Iterator<Item = u64>>(pm_iter: I) -> usize {
1377    let mut s = [!0u64];
1378    for pm in pm_iter {
1379        lcs_step(&mut s, false, |_| pm);
1380    }
1381    lcs_len(&s)
1382}
1383
1384/// Multi-word LCS over `W` words; returns their share of the LCS length.
1385/// With `carries`, the kernel is one strip of a taller pattern (see
1386/// [`kernels::lcs_full`]): it reads each column's incoming carry from there
1387/// and writes the outgoing one back.
1388#[inline(always)]
1389fn lcs_multiword_kernel<const W: usize, I: Iterator<Item = [u64; W]>>(
1390    pm_iter: I,
1391    mut carries: Option<&mut [u8]>,
1392) -> usize {
1393    let mut s = [!0u64; W];
1394    for (j, pm_row) in pm_iter.enumerate() {
1395        let carry_in = carries.as_ref().is_some_and(|c| c[j] != 0);
1396        let carry = lcs_step(&mut s, carry_in, |k| pm_row[k]);
1397        if let Some(c) = carries.as_deref_mut() {
1398            c[j] = carry as u8;
1399        }
1400    }
1401    lcs_len(&s)
1402}
1403
1404// ---------------------------------------------------------------------------
1405// Kernel selection: Levenshtein vs. indel (LCS) distance
1406// ---------------------------------------------------------------------------
1407
1408/// Chooses the bit-parallel recurrence run by the shared peq-building
1409/// pipelines.  Every method returns an *edit distance*, so the pipelines,
1410/// their affix stripping, and their empty-pattern shortcut are identical for
1411/// both metrics (indel distance is affix-invariant for the same reason
1412/// Levenshtein distance is).
1413trait Kernel {
1414    /// The copy of the hot path the pipeline runs in.
1415    type Isa: Isa;
1416
1417    /// Whether the mbleven fast path applies; it is Levenshtein-specific.
1418    const MBLEVEN: bool;
1419
1420    /// Single-word kernel (`m <= 64`).
1421    fn word<I: Iterator<Item = u64>>(m: usize, n: usize, pm_iter: I) -> usize;
1422
1423    /// Multi-word kernel (`64 < m <= 64 * MAX_W`), unrolled over
1424    /// `W = ceil(m / 64)`.
1425    fn words<const W: usize, I: Iterator<Item = [u64; W]>>(m: usize, n: usize, pm_iter: I)
1426        -> usize;
1427
1428    /// Heap-backed kernel for very long patterns (`w > MAX_W`).
1429    fn large<F: Fn(usize) -> usize>(ctx: &LargeCtx<'_, F>) -> usize;
1430}
1431
1432/// Levenshtein distance (substitution, insertion, and deletion each cost 1),
1433/// in the copy `X` of the hot path.
1434struct Lev<X = Baseline>(PhantomData<X>);
1435
1436impl<X: Isa> Kernel for Lev<X> {
1437    type Isa = X;
1438
1439    const MBLEVEN: bool = true;
1440
1441    #[inline(always)]
1442    fn word<I: Iterator<Item = u64>>(m: usize, n: usize, pm_iter: I) -> usize {
1443        hyrro_inner(m, n, pm_iter)
1444    }
1445
1446    #[inline(always)]
1447    fn words<const W: usize, I: Iterator<Item = [u64; W]>>(
1448        m: usize,
1449        n: usize,
1450        pm_iter: I,
1451    ) -> usize {
1452        (n as isize + kernels::multiword::<X, W, I>(m, pm_iter, None)) as usize
1453    }
1454
1455    #[inline(always)]
1456    fn large<F: Fn(usize) -> usize>(ctx: &LargeCtx<'_, F>) -> usize {
1457        // In a copy with target features, LLVM keeps `run` (which has none)
1458        // out of line, so its body is inlined by hand.
1459        #[cfg(target_arch = "x86_64")]
1460        if X::X86_64_V2 {
1461            return ctx.run_body::<X>();
1462        }
1463        ctx.run::<X>()
1464    }
1465}
1466
1467/// Indel distance: insertions and deletions cost 1, substitution is not an
1468/// operation (it decomposes into one of each, cost 2).  Derived from the LCS
1469/// length as `m + n - 2 * LLCS`.  In the copy `X` of the hot path.
1470struct Lcs<X = Baseline>(PhantomData<X>);
1471
1472impl<X: Isa> Kernel for Lcs<X> {
1473    type Isa = X;
1474
1475    const MBLEVEN: bool = false;
1476
1477    #[inline(always)]
1478    fn word<I: Iterator<Item = u64>>(m: usize, n: usize, pm_iter: I) -> usize {
1479        m + n - 2 * lcs_inner_64(pm_iter)
1480    }
1481
1482    #[inline(always)]
1483    fn words<const W: usize, I: Iterator<Item = [u64; W]>>(
1484        m: usize,
1485        n: usize,
1486        pm_iter: I,
1487    ) -> usize {
1488        m + n - 2 * lcs_multiword_kernel::<W, I>(pm_iter, None)
1489    }
1490
1491    #[inline(always)]
1492    fn large<F: Fn(usize) -> usize>(ctx: &LargeCtx<'_, F>) -> usize {
1493        ctx.m + ctx.n - 2 * ctx.lcs_run::<X>()
1494    }
1495}
1496
1497// ---------------------------------------------------------------------------
1498// Multi-word Hyyrö (m > 64)
1499// ---------------------------------------------------------------------------
1500
1501/// Widest pattern, in 64-bit words, handled by the const-width kernels with
1502/// stack peq tables; longer patterns go through [`LargeCtx`].
1503const MAX_W: usize = 8;
1504
1505/// Multi-word entry point with a direct peq table of `ROWS` rows: UCS-1
1506/// pairs (128 rows for pure ASCII, 256 for Latin-1) and mixed-kind pairs
1507/// with a UCS-1 side (257 rows, see [`direct_row`]).  Each copy of the hot
1508/// path compiles its own (see [`Isa`]).
1509#[inline(never)]
1510fn hyrro_multiword_direct<K: Kernel, const ROWS: usize, T1: CodeUnit, T2: CodeUnit>(
1511    short: &[T1],
1512    long: &[T2],
1513) -> usize {
1514    #[cfg(target_arch = "x86_64")]
1515    if K::Isa::X86_64_V2 {
1516        // SAFETY: the CPU runs `K::Isa` (see `Isa`).
1517        return unsafe {
1518            kernels::x86_64::v2::hyrro_multiword_direct::<K, ROWS, T1, T2>(short, long)
1519        };
1520    }
1521    hyrro_multiword_direct_body::<K, ROWS, T1, T2>(short, long)
1522}
1523
1524/// The body of [`hyrro_multiword_direct`].
1525#[inline(always)]
1526fn hyrro_multiword_direct_body<K: Kernel, const ROWS: usize, T1: CodeUnit, T2: CodeUnit>(
1527    short: &[T1],
1528    long: &[T2],
1529) -> usize {
1530    debug_assert!(short.len() > 64);
1531    let m = short.len();
1532    let w = m.div_ceil(64);
1533
1534    macro_rules! run {
1535        ($W:literal) => {{
1536            // An array-of-arrays literal is zeroed row by row in a loop;
1537            // `zeroed` lowers to one `bzero` call.
1538            // SAFETY: all-zero bytes are a valid `[[u64; W]; ROWS]`.
1539            let mut peq: [[u64; $W]; ROWS] = unsafe { std::mem::zeroed() };
1540            for (i, &c) in short.iter().enumerate() {
1541                // SAFETY: direct_row < ROWS and i / 64 < $W.
1542                unsafe {
1543                    *peq.get_unchecked_mut(direct_row::<_, ROWS>(c))
1544                        .get_unchecked_mut(i / 64) |= 1u64 << (i % 64);
1545                }
1546            }
1547            K::words::<$W, _>(
1548                m,
1549                long.len(),
1550                // SAFETY: direct_row < ROWS.
1551                long.iter()
1552                    .map(|&c| unsafe { *peq.get_unchecked(direct_row::<_, ROWS>(c)) }),
1553            )
1554        }};
1555    }
1556    match w {
1557        2 => run!(2),
1558        3 => run!(3),
1559        4 => run!(4),
1560        5 => run!(5),
1561        6 => run!(6),
1562        7 => run!(7),
1563        8 => run!(8),
1564        _ => {
1565            let mut peq = vec![0u64; ROWS * w];
1566            for (i, &c) in short.iter().enumerate() {
1567                peq[direct_row::<_, ROWS>(c) * w + i / 64] |= 1u64 << (i % 64);
1568            }
1569            K::large(&LargeCtx {
1570                m,
1571                n: long.len(),
1572                w,
1573                data: &peq,
1574                // SAFETY: the kernels only ask for j < n.
1575                base_of: |j: usize| direct_row::<_, ROWS>(unsafe { *long.get_unchecked(j) }) * w,
1576            })
1577        }
1578    }
1579}
1580
1581// ---------------------------------------------------------------------------
1582// Banded multi-word kernels (w > MAX_W)
1583// ---------------------------------------------------------------------------
1584
1585/// Band parameter for the next pass after a band of parameter `t` gave up at
1586/// text column `j` of `n`.  The band's bound on the distance grows about
1587/// linearly with `j`, so extrapolating it to `n` (with a margin) estimates
1588/// the distance; at least double `t` so that the passes grow geometrically.
1589#[inline(always)]
1590fn next_band(t: usize, j: usize, n: usize) -> usize {
1591    (t * n / j * 5 / 4).max(2 * t)
1592}
1593
1594/// Number of 64-row blocks a Levenshtein band of half-width `t` can touch in
1595/// one column.
1596#[inline(always)]
1597fn band_blocks(t: usize) -> usize {
1598    2 * t / 64 + 2
1599}
1600
1601/// Number of words an LCS band of parameter `t` (diagonals `[-t, slack + t]`)
1602/// can touch in one column.
1603#[inline(always)]
1604fn lcs_band_words(slack: usize, t: usize) -> usize {
1605    (slack + 2 * t) / 64 + 2
1606}
1607
1608/// Shared context for the very-long-pattern kernels (w > MAX_W): `data` is
1609/// the flat peq table with stride `w`, and `base_of(j)` is the offset of
1610/// text char `j`'s row in it.
1611struct LargeCtx<'a, F> {
1612    m: usize,
1613    n: usize,
1614    w: usize,
1615    data: &'a [u64],
1616    base_of: F,
1617}
1618
1619/// One block of the Myers column step in [`LargeCtx::banded_body`] for match
1620/// mask `eq`, incoming horizontal delta `hin` and vertical deltas `pv`/`mv`:
1621/// returns the new `(pv, mv)` and the block's horizontal deltas `(ph, mh)`.
1622#[inline(always)]
1623fn banded_step(eq: u64, hin: i32, pv: u64, mv: u64) -> ((u64, u64), (u64, u64)) {
1624    let xv = eq | mv;
1625    let eq_in = eq | ((hin < 0) as u64);
1626    let xh = ((eq_in & pv).wrapping_add(pv) ^ pv) | eq_in;
1627    let ph = mv | !(xh | pv);
1628    let mh = pv & xh;
1629    let ph_s = (ph << 1) | ((hin > 0) as u64);
1630    let mh_s = (mh << 1) | ((hin < 0) as u64);
1631    ((mh_s | !(xv | ph_s), ph_s & xv), (ph, mh))
1632}
1633
1634impl<F: Fn(usize) -> usize> LargeCtx<'_, F> {
1635    /// Copy of `self` that reads row offsets from `bases`, one per text char.
1636    /// Resolving rows once pays off when the band passes sweep the text
1637    /// several times; `bases` costs a word per text char, so only build it
1638    /// for those.
1639    fn with_bases<'b>(&'b self, bases: &'b [usize]) -> LargeCtx<'b, impl Fn(usize) -> usize + 'b> {
1640        debug_assert_eq!(bases.len(), self.n);
1641        LargeCtx {
1642            m: self.m,
1643            n: self.n,
1644            w: self.w,
1645            data: self.data,
1646            // SAFETY: the kernels only ask for j < n = bases.len().
1647            base_of: move |j: usize| unsafe { *bases.get_unchecked(j) },
1648        }
1649    }
1650
1651    /// Calls `f` with the row offset of each text char in order.  The
1652    /// offsets are resolved a chunk at a time into a stack buffer, so the
1653    /// lookups (hash probes for UCS-2/4) stay out of the kernel's serial
1654    /// loop without a buffer the size of the text.
1655    #[inline(always)]
1656    fn for_each_base(&self, mut f: impl FnMut(usize)) {
1657        const CHUNK: usize = 512;
1658        let mut buf = [0usize; CHUNK];
1659        for start in (0..self.n).step_by(CHUNK) {
1660            let chunk = &mut buf[..(self.n - start).min(CHUNK)];
1661            for (j, b) in (start..).zip(chunk.iter_mut()) {
1662                *b = (self.base_of)(j);
1663            }
1664            chunk.iter().for_each(|&base| f(base));
1665        }
1666    }
1667
1668    /// Peq word `k` of the row starting at `base`.
1669    #[inline(always)]
1670    fn peq(&self, base: usize, k: usize) -> u64 {
1671        debug_assert!(k < self.w && base + k < self.data.len());
1672        // SAFETY: every base starts a length-w row of `data`, and k < w.
1673        unsafe { *self.data.get_unchecked(base + k) }
1674    }
1675
1676    /// Ukkonen band of half-width `t` over Myers' blocked recurrence: blocks
1677    /// exchange only the horizontal delta `hin`/`hout` (Edlib's formulation),
1678    /// so per column only blocks intersecting rows `[j - t, j + t]` run.
1679    ///
1680    /// Out-of-band boundaries are pessimistic (`hin = +1` below a dropped
1681    /// block; `+1`-per-row extension above entering blocks), so every computed
1682    /// cell is `>=` its true value, while cells on any optimal path are exact
1683    /// whenever the true distance is `<= t` (such paths never leave the band).
1684    /// Returns `Ok(distance)` when the result proves itself (`score <= t`),
1685    /// else `Err(j)` with the text column where the band gave up.
1686    ///
1687    /// Kept out of line: inlined into the large entry functions, the block
1688    /// loop reloads the `pv`, `mv` and peq pointers from the stack on x86-64
1689    /// and loses the branch around the score update on AArch64.  Each copy
1690    /// `X` of the hot path compiles its own (see [`Isa`]).
1691    #[inline(never)]
1692    fn banded<X: Isa>(&self, t: usize, pv: &mut [u64], mv: &mut [u64]) -> Result<usize, usize> {
1693        #[cfg(target_arch = "x86_64")]
1694        if X::X86_64_V2 {
1695            // SAFETY: the CPU runs `X` (see `Isa`).
1696            return unsafe { kernels::x86_64::v2::banded(self, t, pv, mv) };
1697        }
1698        self.banded_body(t, pv, mv)
1699    }
1700
1701    /// The body of [`Self::banded`].
1702    #[inline(always)]
1703    fn banded_body(&self, t: usize, pv: &mut [u64], mv: &mut [u64]) -> Result<usize, usize> {
1704        let (m, n, w) = (self.m, self.n, self.w);
1705        debug_assert!(n >= m);
1706        debug_assert!(t >= n - m);
1707        let (pv, mv) = (&mut pv[..w], &mut mv[..w]);
1708        // Lowest row covered when `k` is the last live block.
1709        let bottom = |k: usize| (64 * (k + 1)).min(m);
1710        let mut last = (t / 64).min(w - 1);
1711        let mut first = 0usize;
1712        pv[..=last].fill(!0);
1713        mv[..=last].fill(0);
1714        let mut score = bottom(last) as isize;
1715        let mut top_bit = (bottom(last) - 1) & 63;
1716
1717        for j in 1..=n {
1718            let base = (self.base_of)(j - 1);
1719            let needed = ((j + t - 1) / 64).min(w - 1);
1720            if needed > last {
1721                pv[last + 1..=needed].fill(!0);
1722                mv[last + 1..=needed].fill(0);
1723                score += (bottom(needed) - bottom(last)) as isize;
1724                last = needed;
1725                top_bit = (bottom(last) - 1) & 63;
1726            }
1727            if j > t + 1 {
1728                first = ((j - t - 1) / 64).min(last);
1729            }
1730            let mut hin: i32 = 1;
1731            debug_assert!(first <= last && last < w);
1732            for k in first..last {
1733                let eq = self.peq(base, k);
1734                // SAFETY: k < last < w = pv.len() = mv.len().
1735                let (pv_k, mv_k) = unsafe { (*pv.get_unchecked(k), *mv.get_unchecked(k)) };
1736                let ((pv_n, mv_n), (ph, mh)) = banded_step(eq, hin, pv_k, mv_k);
1737                hin = ((ph >> 63) & 1) as i32 - ((mh >> 63) & 1) as i32;
1738                // SAFETY: as for the reads above.
1739                unsafe {
1740                    *pv.get_unchecked_mut(k) = pv_n;
1741                    *mv.get_unchecked_mut(k) = mv_n;
1742                }
1743            }
1744            // The final iteration `k = last` is peeled: only it updates
1745            // `score` at `top_bit`, and its outgoing carry is not needed.
1746            let eq = self.peq(base, last);
1747            // SAFETY: last < w = pv.len() = mv.len().
1748            let (pv_l, mv_l) = unsafe { (*pv.get_unchecked(last), *mv.get_unchecked(last)) };
1749            let ((pv_n, mv_n), (ph, mh)) = banded_step(eq, hin, pv_l, mv_l);
1750            score += ((ph >> top_bit) & 1) as isize - ((mh >> top_bit) & 1) as isize;
1751            // SAFETY: as for the reads above.
1752            unsafe {
1753                *pv.get_unchecked_mut(last) = pv_n;
1754                *mv.get_unchecked_mut(last) = mv_n;
1755            }
1756            // `score` never underestimates D[bottom(last)][j], and finishing
1757            // from there takes at least `score - excess` edits (D never
1758            // decreases along diagonal steps); once that exceeds `t` this band
1759            // can no longer prove a result, so bail out and let the caller
1760            // widen or fall back.
1761            let excess = (n - j) as isize - (m - bottom(last)) as isize;
1762            if score - excess.max(0) > t as isize {
1763                return Err(j);
1764            }
1765        }
1766        debug_assert_eq!(last, w - 1);
1767        if score as usize <= t {
1768            Ok(score as usize)
1769        } else {
1770            Err(n)
1771        }
1772    }
1773
1774    /// Peq rows restricted to words `w0..w0 + W`, one per text char.
1775    #[inline(always)]
1776    fn rows<const W: usize>(&self, w0: usize) -> impl Iterator<Item = [u64; W]> + '_ {
1777        debug_assert!(w0 + W <= self.w);
1778        (0..self.n).map(move |j| {
1779            let base = (self.base_of)(j);
1780            std::array::from_fn(|k| self.peq(base, w0 + k))
1781        })
1782    }
1783
1784    /// Bands of growing width before the full matrix.  Each pass aborts as
1785    /// soon as it can no longer prove a result, so a miss on dissimilar
1786    /// strings costs a few percent of the full kernel.  The bands run in
1787    /// the copy `X` of the hot path.
1788    fn run<X: Isa>(&self) -> usize {
1789        self.run_body::<X>()
1790    }
1791
1792    /// The body of [`Self::run`], which the copies other than [`Baseline`]
1793    /// inline (see [`Lev::large`]).
1794    #[inline(always)]
1795    fn run_body<X: Isa>(&self) -> usize {
1796        let mut t = (self.n - self.m + 32).max(64);
1797        if band_blocks(t) * 3 > self.w {
1798            return kernels::full(self);
1799        }
1800        let bases: Vec<usize> = (0..self.n).map(&self.base_of).collect();
1801        let ctx = self.with_bases(&bases);
1802        let (mut pv, mut mv) = (vec![0u64; self.w], vec![0u64; self.w]);
1803        while band_blocks(t) * 3 <= self.w {
1804            match ctx.banded::<X>(t, &mut pv, &mut mv) {
1805                Ok(d) => return d,
1806                Err(j) => t = next_band(t, j, self.n),
1807            }
1808        }
1809        kernels::full(&ctx)
1810    }
1811
1812    /// LCS over the matches on diagonals `j - i` in `[-t, n - m + t]`, rows
1813    /// widened to whole words, computed by running only the words that
1814    /// intersect that band in each column.  Returns `Ok(lcs)` when the band
1815    /// certifies it, else `Err(j)` with the text column where it gave up.
1816    ///
1817    /// Skipping a word is the same as zeroing its pattern-match bits for that
1818    /// column: words below the band are still all ones, and all-ones words
1819    /// without matches only pass the carry through; words above it, without
1820    /// matches and with no carry from row 0, stay unchanged and emit no
1821    /// carry.  The band only moves down, so `l` is the exact LCS over a
1822    /// superset of the in-band matches, hence `l <= LLCS`.  An alignment using
1823    /// any match off the band costs more than `n - m + 2t` indels, so `l` is
1824    /// exact whenever `m + n - 2l <= n - m + 2t`.
1825    fn lcs_banded(&self, t: usize, s: &mut [u64]) -> Result<usize, usize> {
1826        let (m, n, w) = (self.m, self.n, self.w);
1827        debug_assert!(n >= m);
1828        let slack = n - m;
1829        let s = &mut s[..w];
1830        s.fill(!0);
1831        for j in 1..=n {
1832            let base = (self.base_of)(j - 1);
1833            // Rows of the band in column j: `[j - 1 - slack - t, j - 1 + t]`.
1834            let first = j.saturating_sub(slack + t + 1) / 64;
1835            let last = ((j + t - 1) / 64).min(w - 1);
1836            lcs_step(&mut s[first..=last], false, |k| self.peq(base, first + k));
1837            if j % 32 == 0 && self.lcs_band_hopeless(j, t, s) {
1838                return Err(j);
1839            }
1840        }
1841        let l = lcs_len(s);
1842        if m + n - 2 * l <= slack + 2 * t {
1843            Ok(l)
1844        } else {
1845            Err(n)
1846        }
1847    }
1848
1849    /// Whether the band pass can no longer certify (see [`Self::lcs_banded`]).
1850    ///
1851    /// Every alignment crosses column j at some row r and then costs at least
1852    /// `D(r, j) + |(m - r) - (n - j)|`, where `D(r, j) = r + j - 2 * P(r)` and
1853    /// `P(r)` counts LCS increments in rows < r.  As P grows by at most one
1854    /// per row, the minimum over r sits at `r = j - slack`, where the bound
1855    /// is `2j - slack - 2 * P(j - slack)`; for `j < slack` it sits at row 0
1856    /// and is just `slack`, which never exceeds the limit.
1857    fn lcs_band_hopeless(&self, j: usize, t: usize, s: &[u64]) -> bool {
1858        let slack = self.n - self.m;
1859        let Some(c) = j.checked_sub(slack) else {
1860            return false;
1861        };
1862        let (full, part) = (c / 64, c % 64);
1863        let mut p = lcs_len(&s[..full]);
1864        if part != 0 {
1865            p += (!s[full] & ((1u64 << part) - 1)).count_ones() as usize;
1866        }
1867        // p <= c, so this cannot underflow.
1868        2 * j - slack - 2 * p > slack + 2 * t
1869    }
1870
1871    /// LCS length via bands of growing width, like [`Self::run`].
1872    fn lcs_run<X: Isa>(&self) -> usize {
1873        let slack = self.n - self.m;
1874        let mut t = 64;
1875        if lcs_band_words(slack, t) * 3 > self.w {
1876            return kernels::lcs_full::<X, _>(self);
1877        }
1878        let bases: Vec<usize> = (0..self.n).map(&self.base_of).collect();
1879        let ctx = self.with_bases(&bases);
1880        let mut s = vec![0u64; self.w];
1881        while lcs_band_words(slack, t) * 3 <= self.w {
1882            match ctx.lcs_banded(t, &mut s) {
1883                Ok(l) => return l,
1884                Err(j) => t = next_band(t, j, self.n),
1885            }
1886        }
1887        kernels::lcs_full::<X, _>(&ctx)
1888    }
1889}
1890
1891/// `len` copies of `value`, in `stack` when they fit and in `heap` otherwise;
1892/// only the used prefix of `stack` is initialized.
1893#[inline(always)]
1894fn scratch<'a, const N: usize>(
1895    stack: &'a mut [MaybeUninit<u64>; N],
1896    heap: &'a mut Vec<u64>,
1897    len: usize,
1898    value: u64,
1899) -> &'a mut [u64] {
1900    if len <= N {
1901        let buf = &mut stack[..len];
1902        for x in buf.iter_mut() {
1903            x.write(value);
1904        }
1905        // SAFETY: All `len` elements in `buf` were initialized above with `x.write(value)`.
1906        unsafe { std::slice::from_raw_parts_mut(buf.as_mut_ptr().cast::<u64>(), len) }
1907    } else {
1908        *heap = vec![value; len];
1909        heap
1910    }
1911}
1912
1913/// Multi-word entry point for UCS-2/4 and mixed kinds: the peq rows are
1914/// indexed through an open-addressing hash of the pattern alphabet.  Each
1915/// copy of the hot path compiles its own (see [`Isa`]).
1916#[inline(never)]
1917fn hyrro_multiword_hash<K: Kernel, T1: CodeUnit, T2: CodeUnit>(short: &[T1], long: &[T2]) -> usize {
1918    #[cfg(target_arch = "x86_64")]
1919    if K::Isa::X86_64_V2 {
1920        // SAFETY: the CPU runs `K::Isa` (see `Isa`).
1921        return unsafe { kernels::x86_64::v2::hyrro_multiword_hash::<K, T1, T2>(short, long) };
1922    }
1923    hyrro_multiword_hash_body::<K, T1, T2>(short, long)
1924}
1925
1926/// The body of [`hyrro_multiword_hash`].
1927#[inline(always)]
1928fn hyrro_multiword_hash_body<K: Kernel, T1: CodeUnit, T2: CodeUnit>(
1929    short: &[T1],
1930    long: &[T2],
1931) -> usize {
1932    debug_assert!(short.len() > 64);
1933    let m = short.len();
1934    let w = m.div_ceil(64);
1935
1936    // Dense-index the pattern alphabet in first-seen order with a single
1937    // open-addressing pass; each live entry packs (dense index << 32) | key.
1938    // Code units fit in 32 bits and the dense index is bounded by the string
1939    // length, so no live entry can equal the EMPTY sentinel.
1940    const EMPTY: u64 = u64::MAX;
1941    const STACK_M: usize = 64 * MAX_W;
1942    let hash_size = (m * 2).next_power_of_two();
1943    let hash_mask = hash_size - 1;
1944    let hshift = 64 - hash_mask.count_ones();
1945    // Both tables fit on the stack for every const-width kernel; two heap
1946    // round trips would cost as much as the kernel on ~100 chars.
1947    let (mut hash_stack, mut hash_heap) = ([MaybeUninit::uninit(); 2 * STACK_M], Vec::new());
1948    let hash = scratch(&mut hash_stack, &mut hash_heap, hash_size, EMPTY);
1949    let mut ki_stack = [0u16; STACK_M];
1950    let mut ki_heap = Vec::new();
1951    let kis: &mut [u16] = if m <= STACK_M {
1952        &mut ki_stack[..m]
1953    } else {
1954        ki_heap.resize(m, 0);
1955        &mut ki_heap[..]
1956    };
1957    let mut n_keys = 0u64;
1958    for (i, &c) in short.iter().enumerate() {
1959        let key = c.as_u64();
1960        debug_assert!(
1961            key <= u64::from(u32::MAX),
1962            "CodeUnit value must fit in 32 bits"
1963        );
1964        let mut slot = hslot(key, hshift);
1965        loop {
1966            let entry = hash[slot];
1967            if entry == EMPTY {
1968                let ki = n_keys;
1969                hash[slot] = ki << 32 | key;
1970                kis[i] = ki as u16;
1971                n_keys += 1;
1972                break;
1973            }
1974            if entry & 0xFFFF_FFFF == key {
1975                kis[i] = (entry >> 32) as u16;
1976                break;
1977            }
1978            slot = (slot + 1) & hash_mask;
1979        }
1980    }
1981    let n_keys = n_keys as usize;
1982    let hash = &*hash;
1983
1984    // Dense row of a code unit; `n_keys` for units absent from the pattern.
1985    let row_of = |key: u64| -> usize {
1986        let mut slot = hslot(key, hshift);
1987        loop {
1988            // SAFETY: `slot < hash_size == hash.len()`: `hslot` with
1989            // `hshift = 64 - log2(hash_size)` keeps only the top `log2(hash_size)` bits,
1990            // and each probe step is masked with `hash_mask = hash_size - 1`.
1991            let entry = unsafe { *hash.get_unchecked(slot) };
1992            if entry == EMPTY {
1993                return n_keys;
1994            }
1995            if entry & 0xFFFF_FFFF == key {
1996                return (entry >> 32) as usize;
1997            }
1998            slot = (slot + 1) & hash_mask;
1999        }
2000    };
2001
2002    // Rows sized by the actual alphabet; the trailing row stays all-zero and
2003    // serves text chars that do not occur in the pattern.
2004    let (mut data_stack, mut data_heap) =
2005        ([MaybeUninit::uninit(); (STACK_M + 1) * MAX_W], Vec::new());
2006    let data = scratch(&mut data_stack, &mut data_heap, (n_keys + 1) * w, 0);
2007    for (i, &ki) in kis.iter().enumerate() {
2008        let ki = ki as usize;
2009        debug_assert!(ki < n_keys, "pattern key must already be inserted");
2010        data[ki * w + i / 64] |= 1u64 << (i % 64);
2011    }
2012
2013    macro_rules! run {
2014        ($W:literal) => {{
2015            let (rows, _) = data.as_chunks::<$W>();
2016            K::words::<$W, _>(
2017                m,
2018                long.len(),
2019                // SAFETY: row_of returns <= n_keys < rows.len().
2020                long.iter()
2021                    .map(|&c| unsafe { *rows.get_unchecked(row_of(c.as_u64())) }),
2022            )
2023        }};
2024    }
2025    match w {
2026        2 => run!(2),
2027        3 => run!(3),
2028        4 => run!(4),
2029        5 => run!(5),
2030        6 => run!(6),
2031        7 => run!(7),
2032        8 => run!(8),
2033        _ => {
2034            K::large(&LargeCtx {
2035                m,
2036                n: long.len(),
2037                w,
2038                data,
2039                // SAFETY: the kernels only ask for j < n.
2040                base_of: |j: usize| row_of(unsafe { long.get_unchecked(j) }.as_u64()) * w,
2041            })
2042        }
2043    }
2044}
2045
2046// ---------------------------------------------------------------------------
2047// Test helpers (bypass the PyO3 layer)
2048// ---------------------------------------------------------------------------
2049
2050#[cfg(test)]
2051fn metric<K: Kernel>(s1: &str, s2: &str) -> usize {
2052    // Normalize: s1 is shorter pattern, s2 is longer text.
2053    let (s1, s2) = if s1.chars().count() <= s2.chars().count() {
2054        (s1, s2)
2055    } else {
2056        (s2, s1)
2057    };
2058    if s1.is_ascii() && s2.is_ascii() {
2059        compute_u8::<K, true>(s1.as_bytes(), s2.as_bytes())
2060    } else {
2061        let a: Vec<u32> = s1.chars().map(|c| c as u32).collect();
2062        let b: Vec<u32> = s2.chars().map(|c| c as u32).collect();
2063        compute_sorted::<K, _>(&a, &b)
2064    }
2065}
2066
2067#[cfg(test)]
2068fn levenshtein(s1: &str, s2: &str) -> usize {
2069    metric::<Lev>(s1, s2)
2070}
2071
2072#[cfg(test)]
2073fn indel(s1: &str, s2: &str) -> usize {
2074    metric::<Lcs>(s1, s2)
2075}
2076
2077// ---------------------------------------------------------------------------
2078// Tests
2079// ---------------------------------------------------------------------------
2080
2081#[cfg(test)]
2082mod tests {
2083    use super::*;
2084
2085    /// Reference O(m·n) DP — small, obviously correct, used as oracle.
2086    pub(crate) fn naive(a: &[char], b: &[char]) -> usize {
2087        let (m, n) = (a.len(), b.len());
2088        let mut dp = vec![vec![0usize; n + 1]; m + 1];
2089        for (i, row) in dp.iter_mut().enumerate() {
2090            row[0] = i;
2091        }
2092        for (j, val) in dp[0].iter_mut().enumerate() {
2093            *val = j;
2094        }
2095        for i in 1..=m {
2096            for j in 1..=n {
2097                let cost = (a[i - 1] != b[j - 1]) as usize;
2098                dp[i][j] = (dp[i - 1][j] + 1)
2099                    .min(dp[i][j - 1] + 1)
2100                    .min(dp[i - 1][j - 1] + cost);
2101            }
2102        }
2103        dp[m][n]
2104    }
2105
2106    fn check(a: &str, b: &str, expected: usize) {
2107        assert_eq!(levenshtein(a, b), expected, "({a:?}, {b:?})");
2108        assert_eq!(levenshtein(b, a), expected, "({b:?}, {a:?}) symmetry");
2109    }
2110
2111    /// Exhaustive mbleven validation: every pair over a binary alphabet with
2112    /// lengths (9..=11, 9..=11) whose gate bound qualifies must match the
2113    /// oracle exactly.
2114    #[test]
2115    fn mbleven_matches_oracle_exhaustive() {
2116        let strings = |len: usize| -> Vec<Vec<u8>> {
2117            (0u32..1 << len)
2118                .map(|bits| (0..len).map(|i| b'a' + ((bits >> i) & 1) as u8).collect())
2119                .collect()
2120        };
2121        let mut hits = 0usize;
2122        for la in 9..=11usize {
2123            for lb in la..=11usize {
2124                for a in strings(la) {
2125                    for b in strings(lb) {
2126                        if let Some(ub) = small_ub::<Baseline, _>(&a, &b) {
2127                            if ub == 0 {
2128                                continue;
2129                            }
2130                            let ac: Vec<char> = a.iter().map(|&c| c as char).collect();
2131                            let bc: Vec<char> = b.iter().map(|&c| c as char).collect();
2132                            assert_eq!(
2133                                mbleven::<Baseline, _, _>(&a, &b, ub),
2134                                naive(&ac, &bc),
2135                                "({a:?}, {b:?}, ub={ub})"
2136                            );
2137                            hits += 1;
2138                        }
2139                    }
2140                }
2141            }
2142        }
2143        assert!(hits > 100_000, "gate should trigger often here: {hits}");
2144    }
2145
2146    /// Integration: random small-edit pairs (which exercise the small_ub gate
2147    /// through compute) agree with the oracle across lengths and alphabets.
2148    #[test]
2149    fn small_edit_pairs_match_oracle() {
2150        let mut state = 0x243F_6A88_85A3_08D3u64; // xorshift, deterministic
2151        let mut rng = move || {
2152            state ^= state << 13;
2153            state ^= state >> 7;
2154            state ^= state << 17;
2155            state
2156        };
2157        let alphabets: [&[char]; 3] = [
2158            &['a', 'b', 'c', 'd', 'e', 'f'],
2159            &['é', 'ü', 'ø', 'å'],
2160            &['日', '本', '語', 'あ', 'い'],
2161        ];
2162        for _ in 0..4000 {
2163            let ab = alphabets[(rng() % 3) as usize];
2164            let n = 9 + (rng() % 120) as usize;
2165            let a: Vec<char> = (0..n)
2166                .map(|_| ab[(rng() % ab.len() as u64) as usize])
2167                .collect();
2168            let mut b = a.clone();
2169            for _ in 0..(rng() % 4) {
2170                match rng() % 3 {
2171                    0 if !b.is_empty() => {
2172                        let i = (rng() % b.len() as u64) as usize;
2173                        b[i] = ab[(rng() % ab.len() as u64) as usize];
2174                    }
2175                    1 if !b.is_empty() => {
2176                        b.remove((rng() % b.len() as u64) as usize);
2177                    }
2178                    _ => {
2179                        let i = (rng() % (b.len() as u64 + 1)) as usize;
2180                        b.insert(i, ab[(rng() % ab.len() as u64) as usize]);
2181                    }
2182                }
2183            }
2184            let sa: String = a.iter().collect();
2185            let sb: String = b.iter().collect();
2186            assert_eq!(levenshtein(&sa, &sb), naive(&a, &b), "({sa:?}, {sb:?})");
2187        }
2188    }
2189
2190    #[test]
2191    fn empty_strings() {
2192        check("", "", 0);
2193        check("abc", "", 3);
2194        check("", "abc", 3);
2195    }
2196
2197    #[test]
2198    fn equal_strings() {
2199        check("hello", "hello", 0);
2200        check("a", "a", 0);
2201        check("日本語", "日本語", 0);
2202    }
2203
2204    #[test]
2205    fn classic_pairs() {
2206        check("kitten", "sitting", 3);
2207        check("saturday", "sunday", 3);
2208        check("flaw", "lawn", 2);
2209        check("gumbo", "gambol", 2);
2210        check("intention", "execution", 5);
2211    }
2212
2213    #[test]
2214    fn single_edits() {
2215        check("a", "b", 1);
2216        check("a", "ab", 1);
2217        check("ab", "a", 1);
2218        check("abc", "abxc", 1);
2219        check("abc", "axc", 1);
2220    }
2221
2222    #[test]
2223    fn fully_disjoint() {
2224        check("aaaa", "bbbb", 4);
2225        check("abcde", "fghij", 5);
2226    }
2227
2228    #[test]
2229    fn unicode_inputs() {
2230        check("résumé", "resume", 2);
2231        check("café", "cafe", 1);
2232        check("日本語", "日本", 1);
2233        check("🦀🐍", "🐍🦀", 2);
2234        check("naïve", "naive", 1);
2235    }
2236
2237    #[test]
2238    fn boundary_64() {
2239        let a64: String = "a".repeat(64);
2240        let a65: String = "a".repeat(65);
2241        check(&a64, &a64, 0);
2242        check(&a64, &a65, 1);
2243        let mut shifted = String::from("b");
2244        shifted.push_str(&"a".repeat(63));
2245        check(&a64, &shifted, 1);
2246    }
2247
2248    #[test]
2249    fn long_inputs_multiword() {
2250        let a: String = "abc".repeat(40); // 120 chars
2251        let mut b = a.clone();
2252        b.insert(0, 'x');
2253        check(&a, &b, 1);
2254        let c = format!("{a}xyz");
2255        check(&a, &c, 3);
2256    }
2257
2258    /// Directly exercise `hyrro_multiword_direct` (m > 64) against the naive DP.
2259    #[test]
2260    fn multiword_bytes_matches_oracle() {
2261        let oracle = |a: &[u8], b: &[u8]| -> usize {
2262            let ac: Vec<char> = a.iter().map(|&c| c as char).collect();
2263            let bc: Vec<char> = b.iter().map(|&c| c as char).collect();
2264            naive(&ac, &bc)
2265        };
2266        // ASCII bytes (< 128): verify correctness.
2267        let check_ascii = |a: &[u8], b: &[u8]| {
2268            let (s, l) = if a.len() <= b.len() { (a, b) } else { (b, a) };
2269            assert!(s.len() > 64 && s.iter().all(|&c| c < 128) && l.iter().all(|&c| c < 128));
2270            assert_eq!(hyrro_multiword_direct::<Lev, 256, _, _>(s, l), oracle(s, l));
2271        };
2272        // Latin-1 bytes: any u8 value allowed.
2273        let check_latin1 = |a: &[u8], b: &[u8]| {
2274            let (s, l) = if a.len() <= b.len() { (a, b) } else { (b, a) };
2275            assert!(s.len() > 64);
2276            assert_eq!(hyrro_multiword_direct::<Lev, 256, _, _>(s, l), oracle(s, l));
2277        };
2278
2279        // Benchmark-like: two long ASCII strings that diverge after a shared prefix.
2280        let s1 = b"Lets pretend Marshall Mathers never picked up a pen".repeat(8);
2281        let s2 = b"Lets pretend things woulda been no different".repeat(8);
2282        let (sh, lo) = if s1.len() <= s2.len() {
2283            (&s1[..], &s2[..])
2284        } else {
2285            (&s2[..], &s1[..])
2286        };
2287        assert_eq!(
2288            hyrro_multiword_direct::<Lev, 256, _, _>(sh, lo),
2289            oracle(sh, lo)
2290        );
2291
2292        // Fully disjoint long ASCII strings.
2293        check_ascii(&b"a".repeat(100), &b"b".repeat(100));
2294
2295        // Insertions at different positions (cycling through lowercase ASCII).
2296        let base: Vec<u8> = (0u8..80).map(|i| b'a' + i % 26).collect();
2297        let mut ins_front = vec![b'z'];
2298        ins_front.extend_from_slice(&base);
2299        check_ascii(&base, &ins_front);
2300
2301        let mut ins_mid = base[..40].to_vec();
2302        ins_mid.push(b'z');
2303        ins_mid.extend_from_slice(&base[40..]);
2304        check_ascii(&base, &ins_mid);
2305
2306        // Boundary: exactly 65 chars.
2307        check_ascii(&b"x".repeat(65), &b"y".repeat(65));
2308
2309        // Latin-1: bytes spanning full 0..=255 range, 2-word boundary (128 chars).
2310        let p128: Vec<u8> = (0u8..128).collect();
2311        let q128: Vec<u8> = (1u8..=128).collect(); // 128 is a Latin-1 byte
2312        check_latin1(&p128, &q128);
2313
2314        // Latin-1 disjoint.
2315        check_latin1(&[200u8; 80], &[201u8; 80]);
2316    }
2317
2318    /// Mixed-kind pairs with a UCS-1 side take the direct-table kernel with
2319    /// wide units folded into row 256; check both orientations, including a
2320    /// wide pattern whose units above U+00FF land in that row.
2321    #[test]
2322    fn mixed_direct_matches_oracle() {
2323        let mut state = 0x0123_4567_89AB_CDEF_u64;
2324        let narrow = ['a', 'b', 'é', '\u{ff}', ' '];
2325        let wide = ['a', 'é', '\u{100}', '\u{2019}', '\u{ffff}'];
2326        let mut gen = |alpha: &[char], len: usize| -> Vec<char> {
2327            (0..len)
2328                .map(|_| alpha[lcg(&mut state) % alpha.len()])
2329                .collect()
2330        };
2331        for &(m, n) in &[(65usize, 65usize), (100, 130), (300, 300), (520, 700)] {
2332            // UCS-1 pattern, UCS-2 text.
2333            let (a, b) = (gen(&narrow, m), gen(&wide, n));
2334            let a8: Vec<u8> = a.iter().map(|&c| c as u8).collect();
2335            let b16: Vec<u16> = b.iter().map(|&c| c as u16).collect();
2336            assert_eq!(compute_sorted_mixed::<Lev, _, _>(&a8, &b16), naive(&a, &b));
2337            assert_eq!(
2338                compute_sorted_mixed::<Lcs, _, _>(&a8, &b16),
2339                naive_indel(&a, &b)
2340            );
2341            // UCS-2 pattern, UCS-1 text.
2342            let (a, b) = (gen(&wide, m), gen(&narrow, n));
2343            let a16: Vec<u16> = a.iter().map(|&c| c as u16).collect();
2344            let b8: Vec<u8> = b.iter().map(|&c| c as u8).collect();
2345            assert_eq!(compute_sorted_mixed::<Lev, _, _>(&a16, &b8), naive(&a, &b));
2346            assert_eq!(
2347                compute_sorted_mixed::<Lcs, _, _>(&a16, &b8),
2348                naive_indel(&a, &b)
2349            );
2350        }
2351    }
2352
2353    /// Directly exercise `hyrro_multiword_hash` (m > 64) against the naive DP.
2354    #[test]
2355    fn multiword_sorted_matches_oracle() {
2356        let check_u32 = |a: &[u32], b: &[u32]| {
2357            let (s, l) = if a.len() <= b.len() { (a, b) } else { (b, a) };
2358            assert!(
2359                s.len() > 64,
2360                "test case must be long enough to hit multiword"
2361            );
2362            let got = hyrro_multiword_hash::<Lev, _, _>(s, l);
2363            let ac: Vec<char> = s
2364                .iter()
2365                .map(|&c| char::from_u32(c).unwrap_or('?'))
2366                .collect();
2367            let bc: Vec<char> = l
2368                .iter()
2369                .map(|&c| char::from_u32(c).unwrap_or('?'))
2370                .collect();
2371            let exp = naive(&ac, &bc);
2372            assert_eq!(got, exp, "u32 mismatch ({} vs {} chars)", s.len(), l.len());
2373        };
2374
2375        // Emoji run: each emoji is one code point — build 80-char runs.
2376        let emoji_a: Vec<u32> = [0x1f980u32, 0x1f40d, 0x1f389, 0x1f38a, 0x1f388]
2377            .iter()
2378            .copied()
2379            .cycle()
2380            .take(80)
2381            .collect();
2382        let emoji_b: Vec<u32> = [0x1f40du32, 0x1f980, 0x1f389, 0x1f38a, 0x1f388]
2383            .iter()
2384            .copied()
2385            .cycle()
2386            .take(80)
2387            .collect();
2388        check_u32(&emoji_a, &emoji_b);
2389
2390        // CJK run: 80 characters from the Japanese test string.
2391        let cjk_a: Vec<u32> = "日本語のテスト文字列"
2392            .chars()
2393            .map(|c| c as u32)
2394            .cycle()
2395            .take(80)
2396            .collect();
2397        let cjk_b: Vec<u32> = "日本語のテスツ文字列"
2398            .chars()
2399            .map(|c| c as u32)
2400            .cycle()
2401            .take(80)
2402            .collect();
2403        check_u32(&cjk_a, &cjk_b);
2404
2405        // Exactly 2 words (128 elements).
2406        let a128: Vec<u32> = (0u32..128).collect();
2407        let b128: Vec<u32> = (1u32..=128).collect();
2408        check_u32(&a128, &b128);
2409
2410        // Fully disjoint.
2411        let all_a: Vec<u32> = vec![1u32; 80];
2412        let all_b: Vec<u32> = vec![2u32; 80];
2413        check_u32(&all_a, &all_b);
2414    }
2415
2416    #[test]
2417    fn affix_stripping_does_not_change_result() {
2418        check("xxx_kitten_yyy", "xxx_sitting_yyy", 3);
2419        check("prefix-foo", "prefix-bar", 3);
2420        check("foo-suffix", "bar-suffix", 3);
2421    }
2422
2423    #[test]
2424    fn strip_affix_unaligned_lengths() {
2425        // One string a prefix of the other, lengths crossing the 8-byte chunk boundary.
2426        for n in 1..=17 {
2427            let a = "a".repeat(n);
2428            let b = "a".repeat(n + 1);
2429            check(&a, &b, 1);
2430            let c = format!("{a}b");
2431            check(&c, &a, 1);
2432        }
2433        // Prefix and suffix scans meeting in the middle.
2434        check("abcdefgh_XY_abcdefgh", "abcdefgh_YX_abcdefgh", 2);
2435        check("abcdefghij", "abcdefghij", 0);
2436    }
2437
2438    /// The block scans must find a difference at every offset inside, at
2439    /// and past a block boundary, from either end.
2440    #[test]
2441    fn affix_scan_finds_every_mismatch() {
2442        for n in [31usize, 32, 33, 63, 64, 65, 100] {
2443            let a: Vec<u8> = (0..n).map(|i| i as u8).collect();
2444            assert_eq!(common_prefix_bytes(&a, &a), n);
2445            assert_eq!(common_suffix_bytes(&a, &a), n);
2446            for pos in 0..n {
2447                let mut b = a.clone();
2448                b[pos] ^= 0x80;
2449                assert_eq!(common_prefix_bytes(&a, &b), pos, "n={n}");
2450                assert_eq!(common_suffix_bytes(&a, &b), n - 1 - pos, "n={n}");
2451            }
2452        }
2453    }
2454
2455    #[test]
2456    fn strip_affix_discards_partial_elements() {
2457        // u16 units where raw bytes match past an element boundary: the
2458        // half-matched element must not be stripped.
2459        let a16: Vec<u16> = vec![0x0101, 0x0102, 0x0201];
2460        let b16: Vec<u16> = vec![0x0101, 0x0202, 0x0201];
2461        let (sa, sb) = strip_affix(&a16, &b16);
2462        assert_eq!(sa, &[0x0102]);
2463        assert_eq!(sb, &[0x0202]);
2464        assert_eq!(compute_sorted::<Lev, _>(&a16, &b16), 1);
2465    }
2466
2467    #[test]
2468    fn tiny_patterns_match_oracle() {
2469        // Deterministic LCG; small alphabet forces shared affixes and hits the
2470        // tiny path (m ≤ TINY_M after stripping) as well as empty-after-strip.
2471        let mut state = 0x243F_6A88_85A3_08D3_u64;
2472        let mut rng = move || {
2473            state = state
2474                .wrapping_mul(6364136223846793005)
2475                .wrapping_add(1442695040888963407);
2476            (state >> 33) as usize
2477        };
2478        for _ in 0..4000 {
2479            let la = rng() % 25;
2480            let lb = rng() % 25;
2481            let a: String = (0..la)
2482                .map(|_| (b'a' + (rng() % 4) as u8) as char)
2483                .collect();
2484            let b: String = (0..lb)
2485                .map(|_| (b'a' + (rng() % 4) as u8) as char)
2486                .collect();
2487            let ac: Vec<char> = a.chars().collect();
2488            let bc: Vec<char> = b.chars().collect();
2489            assert_eq!(levenshtein(&a, &b), naive(&ac, &bc), "({a:?}, {b:?})");
2490        }
2491    }
2492
2493    #[test]
2494    fn tiny_patterns_across_all_kinds() {
2495        let mut state = 0x85A3_08D3_243F_6A88_u64;
2496        let mut rng = move || {
2497            state = state
2498                .wrapping_mul(6364136223846793005)
2499                .wrapping_add(1442695040888963407);
2500            (state >> 33) as usize
2501        };
2502        // Latin-1 (u8 non-ascii)
2503        for _ in 0..1000 {
2504            let la = 1 + (rng() % 16);
2505            let lb = 1 + (rng() % 20);
2506            let a: Vec<u8> = (0..la).map(|_| 0xc0 + (rng() % 10) as u8).collect();
2507            let b: Vec<u8> = (0..lb).map(|_| 0xc0 + (rng() % 10) as u8).collect();
2508            let ac: Vec<char> = a.iter().map(|&c| c as char).collect();
2509            let bc: Vec<char> = b.iter().map(|&c| c as char).collect();
2510            assert_eq!(compute_u8::<Lev, false>(&a, &b), naive(&ac, &bc));
2511        }
2512        // UCS-2 (u16)
2513        for _ in 0..1000 {
2514            let la = 1 + (rng() % 16);
2515            let lb = 1 + (rng() % 20);
2516            let a: Vec<u16> = (0..la).map(|_| 0x0500 + (rng() % 10) as u16).collect();
2517            let b: Vec<u16> = (0..lb).map(|_| 0x0500 + (rng() % 10) as u16).collect();
2518            let ac: Vec<char> = a
2519                .iter()
2520                .map(|&c| char::from_u32(c as u32).unwrap())
2521                .collect();
2522            let bc: Vec<char> = b
2523                .iter()
2524                .map(|&c| char::from_u32(c as u32).unwrap())
2525                .collect();
2526            assert_eq!(compute_sorted::<Lev, _>(&a, &b), naive(&ac, &bc));
2527        }
2528        // UCS-4 (u32)
2529        for _ in 0..1000 {
2530            let la = 1 + (rng() % 16);
2531            let lb = 1 + (rng() % 20);
2532            let a: Vec<u32> = (0..la).map(|_| 0x1f600 + (rng() % 10) as u32).collect();
2533            let b: Vec<u32> = (0..lb).map(|_| 0x1f600 + (rng() % 10) as u32).collect();
2534            let ac: Vec<char> = a.iter().map(|&c| char::from_u32(c).unwrap()).collect();
2535            let bc: Vec<char> = b.iter().map(|&c| char::from_u32(c).unwrap()).collect();
2536            assert_eq!(compute_sorted::<Lev, _>(&a, &b), naive(&ac, &bc));
2537        }
2538        // Mixed: u8 pattern, u16 text
2539        for _ in 0..1000 {
2540            let la = 1 + (rng() % 16);
2541            let lb = 1 + (rng() % 20);
2542            let a: Vec<u8> = (0..la).map(|_| b'a' + (rng() % 6) as u8).collect();
2543            let b: Vec<u16> = (0..lb)
2544                .map(|_| {
2545                    if rng() % 3 == 0 {
2546                        (b'a' + (rng() % 6) as u8) as u16
2547                    } else {
2548                        0x0300 + (rng() % 5) as u16
2549                    }
2550                })
2551                .collect();
2552            let ac: Vec<char> = a.iter().map(|&c| c as char).collect();
2553            let bc: Vec<char> = b
2554                .iter()
2555                .map(|&c| char::from_u32(c as u32).unwrap())
2556                .collect();
2557            assert_eq!(compute_sorted_mixed::<Lev, _, _>(&a, &b), naive(&ac, &bc));
2558        }
2559    }
2560
2561    /// The mixed-kind block scans must find a difference at every offset
2562    /// inside, at and past a block boundary, from either end, and count
2563    /// mismatches like the element-wise definition.
2564    #[test]
2565    fn mixed_scans_find_every_mismatch() {
2566        for n in [15usize, 16, 17, 31, 32, 33, 50] {
2567            let a: Vec<u8> = (0..n).map(|i| b'a' + (i % 26) as u8).collect();
2568            let wide: Vec<u16> = a.iter().map(|&c| c as u16).collect();
2569            assert_eq!(strip_affix_mixed(&a, &wide).0.len(), 0);
2570            for pos in 0..n {
2571                let mut b = wide.clone();
2572                b[pos] |= 0x100;
2573                let (sa, sb) = strip_affix_mixed(&a, &b);
2574                assert_eq!((sa, sb), (&a[pos..=pos], &b[pos..=pos]), "n={n}");
2575                assert_eq!(small_ub_mixed(&a, &b), Some(1), "n={n} pos={pos}");
2576            }
2577            let all: Vec<u16> = wide.iter().map(|&c| c | 0x100).collect();
2578            assert_eq!(small_ub_mixed(&a, &all), None, "n={n}");
2579        }
2580    }
2581
2582    #[test]
2583    fn tiny_patterns_mixed_kinds() {
2584        // Tiny path through the mixed-kind pipeline (u8 pattern, u16 text).
2585        let a: Vec<u8> = b"abc".to_vec();
2586        let b: Vec<u16> = "axc".encode_utf16().collect();
2587        assert_eq!(compute_sorted_mixed::<Lev, _, _>(&a, &b), 1);
2588        let c: Vec<u16> = "xyz".encode_utf16().collect();
2589        assert_eq!(compute_sorted_mixed::<Lev, _, _>(&a, &c), 3);
2590    }
2591
2592    /// Deterministic LCG used by the banded-kernel tests.
2593    pub(crate) fn lcg(state: &mut u64) -> usize {
2594        *state = state
2595            .wrapping_mul(6364136223846793005)
2596            .wrapping_add(1442695040888963407);
2597        (*state >> 33) as usize
2598    }
2599
2600    pub(crate) fn rand_string(state: &mut u64, len: usize, alpha: &[u8]) -> String {
2601        (0..len)
2602            .map(|_| alpha[lcg(state) % alpha.len()] as char)
2603            .collect()
2604    }
2605
2606    pub(crate) fn mutate(state: &mut u64, s: &str, edits: usize, alpha: &[u8]) -> String {
2607        let mut out: Vec<char> = s.chars().collect();
2608        for _ in 0..edits {
2609            let i = lcg(state) % out.len();
2610            match lcg(state) % 3 {
2611                0 => out[i] = alpha[lcg(state) % alpha.len()] as char,
2612                1 => out.insert(i, alpha[lcg(state) % alpha.len()] as char),
2613                _ => {
2614                    out.remove(i);
2615                }
2616            }
2617        }
2618        out.into_iter().collect()
2619    }
2620
2621    #[test]
2622    fn banded_long_strings_match_oracle() {
2623        // m > 512 (w > 8) exercises the banded kernel: small edit counts hit
2624        // the optimistic pass, larger ones the retry/fallback logic.
2625        let mut state = 0xDEAD_BEEF_CAFE_F00D_u64;
2626        let alpha = b"abcdefghij";
2627        for &(len, edits) in &[(520, 1), (700, 3), (900, 30), (1500, 60), (900, 120)] {
2628            let a = rand_string(&mut state, len, alpha);
2629            let b = mutate(&mut state, &a, edits, alpha);
2630            let ac: Vec<char> = a.chars().collect();
2631            let bc: Vec<char> = b.chars().collect();
2632            assert_eq!(
2633                levenshtein(&a, &b),
2634                naive(&ac, &bc),
2635                "len {len} edits {edits}"
2636            );
2637        }
2638        // Substitution-only edits keep the Hamming upper bound tight while
2639        // d > 64, forcing the second (guaranteed) banded pass.
2640        let a = rand_string(&mut state, 700, alpha);
2641        let mut bc: Vec<char> = a.chars().collect();
2642        for _ in 0..100 {
2643            let i = lcg(&mut state) % bc.len();
2644            bc[i] = alpha[lcg(&mut state) % alpha.len()] as char;
2645        }
2646        let b: String = bc.iter().collect();
2647        let ac: Vec<char> = a.chars().collect();
2648        assert_eq!(levenshtein(&a, &b), naive(&ac, &bc));
2649    }
2650
2651    #[test]
2652    fn banded_shifted_and_dissimilar() {
2653        let mut state = 0x1234_5678_9ABC_DEF0_u64;
2654        let alpha = b"abcdefghij";
2655        // A 1-char front shift defeats both affix stripping and the Hamming
2656        // bound, but the optimistic band still catches d = 2 directly.
2657        let a = rand_string(&mut state, 800, alpha);
2658        let b = format!("x{}", &a[..a.len() - 1]);
2659        let (ac, bc): (Vec<char>, Vec<char>) = (a.chars().collect(), b.chars().collect());
2660        assert_eq!(levenshtein(&a, &b), naive(&ac, &bc));
2661        // Unrelated strings: early-abort then full-matrix path.
2662        let c = rand_string(&mut state, 750, alpha);
2663        let d = rand_string(&mut state, 640, alpha);
2664        let (cc, dc): (Vec<char>, Vec<char>) = (c.chars().collect(), d.chars().collect());
2665        assert_eq!(levenshtein(&c, &d), naive(&cc, &dc));
2666        // Large length difference: banding is skipped outright (t0 too wide).
2667        let e = rand_string(&mut state, 600, alpha);
2668        let f = mutate(&mut state, &e, 5, alpha) + &rand_string(&mut state, 700, alpha);
2669        let (ec, fc): (Vec<char>, Vec<char>) = (e.chars().collect(), f.chars().collect());
2670        assert_eq!(levenshtein(&e, &f), naive(&ec, &fc));
2671    }
2672
2673    #[test]
2674    fn banded_non_ascii_long_strings() {
2675        // Banded kernel through the hash-based (non-UCS-1) peq pipeline.
2676        let mut state = 0x0F0F_0F0F_1111_2222_u64;
2677        let alpha: Vec<char> = "ぁあぃいぅうぇえおかがきぎく".chars().collect();
2678        let a: String = (0..640)
2679            .map(|_| alpha[lcg(&mut state) % alpha.len()])
2680            .collect();
2681        let mut bc: Vec<char> = a.chars().collect();
2682        for _ in 0..7 {
2683            let i = lcg(&mut state) % bc.len();
2684            match lcg(&mut state) % 3 {
2685                0 => bc[i] = alpha[lcg(&mut state) % alpha.len()],
2686                1 => bc.insert(i, alpha[lcg(&mut state) % alpha.len()]),
2687                _ => {
2688                    bc.remove(i);
2689                }
2690            }
2691        }
2692        let b: String = bc.iter().collect();
2693        let ac: Vec<char> = a.chars().collect();
2694        assert_eq!(levenshtein(&a, &b), naive(&ac, &bc));
2695    }
2696
2697    #[test]
2698    fn ucs2_pipeline_matches_oracle() {
2699        let cases: &[(&str, &str, usize)] = &[
2700            ("日本語", "日本", 1),
2701            ("日本語のテスト", "日本語のテスツ", 1),
2702            ("한국어", "한국", 1),
2703            ("中文", "中英文", 1),
2704        ];
2705        for &(a, b, expected) in cases {
2706            let au: Vec<u16> = a.encode_utf16().collect();
2707            let bu: Vec<u16> = b.encode_utf16().collect();
2708            assert_eq!(
2709                compute_sorted::<Lev, _>(&au, &bu),
2710                expected,
2711                "u16 ({a:?}, {b:?})"
2712            );
2713            assert_eq!(levenshtein(a, b), expected, "u32 ({a:?}, {b:?})");
2714        }
2715    }
2716
2717    #[test]
2718    fn ucs4_pipeline_matches_oracle() {
2719        let cases: &[(&str, &str, usize)] = &[
2720            ("🦀🐍", "🐍🦀", 2),
2721            ("🎉🎊🎈", "🎊🎈", 1),
2722            ("😀😁😂", "😀😂", 1),
2723        ];
2724        for &(a, b, expected) in cases {
2725            let au: Vec<u32> = a.chars().map(|c| c as u32).collect();
2726            let bu: Vec<u32> = b.chars().map(|c| c as u32).collect();
2727            assert_eq!(
2728                compute_sorted::<Lev, _>(&au, &bu),
2729                expected,
2730                "u32 ({a:?}, {b:?})"
2731            );
2732            assert_eq!(levenshtein(a, b), expected, "u32 ({a:?}, {b:?})");
2733        }
2734    }
2735
2736    #[test]
2737    fn oracle_random_cases() {
2738        let long_a = "a".repeat(70);
2739        let long_b = "a".repeat(65);
2740        let long_c = "a".repeat(80);
2741        let long_d = "b".repeat(80);
2742        let cases: &[(&str, &str)] = &[
2743            ("hello", "world"),
2744            ("kitten", "sitting"),
2745            ("intention", "execution"),
2746            ("abcdefghijklmnopqrstuvwxyz", "zyxwvutsrqponmlkjihgfedcba"),
2747            ("", ""),
2748            ("a", ""),
2749            ("", "a"),
2750            ("aaaa", "aaaa"),
2751            ("abababab", "babababa"),
2752            ("the quick brown fox", "the quik brwn fx"),
2753            ("résumé", "résume"),
2754            ("日本語のテスト", "日本語のテスツ"),
2755            (&long_a, &long_b),
2756            (&long_c, &long_d),
2757        ];
2758        for &(a, b) in cases {
2759            let ac: Vec<char> = a.chars().collect();
2760            let bc: Vec<char> = b.chars().collect();
2761            assert_eq!(levenshtein(a, b), naive(&ac, &bc), "({a:?}, {b:?})");
2762        }
2763    }
2764
2765    #[test]
2766    fn bench_multiword_raw() {
2767        use std::time::Instant;
2768        // Simulate the benchmark strings after affix stripping.
2769        let s1_full = b"Lets pretend things woulda been no different".repeat(8);
2770        let s2_full = b"Lets pretend Marshall Mathers never picked up a pen".repeat(8);
2771        // strip 13-char common prefix "Lets pretend "
2772        let short = &s1_full[13..];
2773        let long = &s2_full[13..];
2774        assert!(short.len() <= long.len());
2775        let n = 10_000u32;
2776        let mut sink = 0usize;
2777        let t0 = Instant::now();
2778        for _ in 0..n {
2779            sink += hyrro_multiword_direct::<Lev, 256, _, _>(short, long);
2780        }
2781        let us = t0.elapsed().as_secs_f64() * 1e6 / n as f64;
2782        #[cfg(debug_assertions)]
2783        eprintln!("NOTE: debug build — release will be ~10-20x faster");
2784        eprintln!(
2785            "raw Rust hyrro_multiword_direct: {:.3} μs/call  (sink={})",
2786            us, sink
2787        );
2788        assert!(sink > 0); // prevent DCE
2789    }
2790
2791    // -----------------------------------------------------------------------
2792    // Indel distance / ratio
2793    // -----------------------------------------------------------------------
2794
2795    /// Reference O(m·n) indel DP (insert and delete only), used as oracle.
2796    pub(crate) fn naive_indel(a: &[char], b: &[char]) -> usize {
2797        let (m, n) = (a.len(), b.len());
2798        let mut dp = vec![vec![0usize; n + 1]; m + 1];
2799        for (i, row) in dp.iter_mut().enumerate() {
2800            row[0] = i;
2801        }
2802        for (j, val) in dp[0].iter_mut().enumerate() {
2803            *val = j;
2804        }
2805        for i in 1..=m {
2806            for j in 1..=n {
2807                dp[i][j] = if a[i - 1] == b[j - 1] {
2808                    dp[i - 1][j - 1]
2809                } else {
2810                    (dp[i - 1][j] + 1).min(dp[i][j - 1] + 1)
2811                };
2812            }
2813        }
2814        dp[m][n]
2815    }
2816
2817    fn check_indel(a: &str, b: &str, expected: usize) {
2818        assert_eq!(indel(a, b), expected, "({a:?}, {b:?})");
2819        assert_eq!(indel(b, a), expected, "({b:?}, {a:?}) symmetry");
2820    }
2821
2822    #[test]
2823    fn indel_basic() {
2824        check_indel("", "", 0);
2825        check_indel("", "abc", 3);
2826        check_indel("abc", "abc", 0);
2827        check_indel("abc", "xyz", 6);
2828        check_indel("a", "b", 2);
2829        // kitten/sitting: LCS "ittn" (4), so 6 + 7 - 8.
2830        check_indel("kitten", "sitting", 5);
2831        check_indel("résumé", "resume", 4);
2832        // Affix stripping must not disturb the result.
2833        check_indel("abcXdef", "abcYdef", 2);
2834        check_indel("prefix_middle_suffix", "prefix_MIDDLE_suffix", 12);
2835    }
2836
2837    #[test]
2838    fn ratio_matches_indel_convention() {
2839        let r = |a: &str, b: &str| -> f64 {
2840            let total = a.chars().count() + b.chars().count();
2841            if total == 0 {
2842                1.0
2843            } else {
2844                1.0 - indel(a, b) as f64 / total as f64
2845            }
2846        };
2847        assert!((r("", "") - 1.0).abs() < 1e-12);
2848        assert!((r("abc", "abc") - 1.0).abs() < 1e-12);
2849        // Values cross-checked against rapidfuzz.fuzz.ratio / Levenshtein.ratio.
2850        assert!((r("kitten", "sitting") - 0.615_384_615_384_615_4).abs() < 1e-12);
2851        assert!((r("résumé", "resume") - 2.0 / 3.0).abs() < 1e-12);
2852        // Sharing no character is exactly 0.0, unlike a max-length or
2853        // Levenshtein-over-sum normalization.
2854        assert!(r("abc", "xyz").abs() < 1e-12);
2855        assert!(r("a", "b").abs() < 1e-12);
2856    }
2857
2858    /// The LCS kernels must agree with the DP across every length regime:
2859    /// tiny (<= 8), single-word (<= 64), multi-word (<= 512), and the
2860    /// heap-backed `lcs_full` path (> 512), over ASCII and astral alphabets.
2861    #[test]
2862    fn indel_matches_oracle_across_length_regimes() {
2863        let mut state = 0x5EED_1234_ABCD_0001_u64;
2864        for alpha in [&b"ab"[..], &b"abcdefghij"[..]] {
2865            for &len in &[1usize, 5, 8, 9, 33, 64, 65, 100, 200, 513, 700] {
2866                for &edits in &[1usize, 3, 12] {
2867                    let a = rand_string(&mut state, len, alpha);
2868                    // `mutate` indexes modulo the current length, so keep at
2869                    // least one character alive even if every edit deletes.
2870                    let b = mutate(&mut state, &a, edits.min(len - 1).max(1), alpha);
2871                    let (ac, bc): (Vec<char>, Vec<char>) =
2872                        (a.chars().collect(), b.chars().collect());
2873                    assert_eq!(indel(&a, &b), naive_indel(&ac, &bc), "({a:?}, {b:?})");
2874
2875                    // Same content promoted to UCS-4 must agree.
2876                    let map = |s: &str| -> String {
2877                        s.chars()
2878                            .map(|c| char::from_u32(0x1_0000 + c as u32).unwrap())
2879                            .collect()
2880                    };
2881                    let (wa, wb) = (map(&a), map(&b));
2882                    assert_eq!(indel(&wa, &wb), naive_indel(&ac, &bc), "astral ({a:?})");
2883                }
2884            }
2885        }
2886    }
2887
2888    /// Unrelated strings of every length regime: no common affix to strip and
2889    /// a disjoint alphabet, so the kernels run over their full width.
2890    #[test]
2891    fn indel_disjoint_alphabets() {
2892        for &len in &[1usize, 8, 64, 65, 300, 600] {
2893            let a: String = std::iter::repeat_n('a', len).collect();
2894            let b: String = std::iter::repeat_n('b', len).collect();
2895            check_indel(&a, &b, 2 * len);
2896        }
2897    }
2898
2899    /// Pairs whose lengths differ several-fold, in every length regime of the
2900    /// shorter string: embedded unchanged (exactly `n - m` apart), unrelated,
2901    /// and mutated inside random padding; ASCII and astral alphabets.
2902    #[test]
2903    fn very_different_lengths_match_oracle() {
2904        let mut state = 0xD1FF_E4E7_0000_0001_u64;
2905        let astral = |s: &str| -> String {
2906            s.chars()
2907                .map(|c| char::from_u32(0x1_0000 + c as u32).unwrap())
2908                .collect()
2909        };
2910        for alpha in [&b"ab"[..], &b"abcdefghij"[..]] {
2911            for &(m, n) in &[
2912                (0usize, 40usize),
2913                (1, 200),
2914                (7, 300),
2915                (9, 150),
2916                (40, 900),
2917                (64, 700),
2918                (65, 1000),
2919                (200, 1300),
2920                (530, 1700),
2921            ] {
2922                let short = rand_string(&mut state, m, alpha);
2923                let pad = rand_string(&mut state, n - m, alpha);
2924                for cut in [0, pad.len() / 3, pad.len()] {
2925                    let long = format!("{}{short}{}", &pad[..cut], &pad[cut..]);
2926                    assert_eq!(levenshtein(&short, &long), n - m, "embedded m={m} n={n}");
2927                }
2928                let edited = if m > 0 {
2929                    mutate(&mut state, &short, m / 10 + 1, alpha)
2930                } else {
2931                    String::new()
2932                };
2933                let unrelated = rand_string(&mut state, n, alpha);
2934                let padded = format!("{}{edited}{}", &pad[..pad.len() / 2], &pad[pad.len() / 2..]);
2935                for long in [unrelated, padded] {
2936                    let (sc, lc): (Vec<char>, Vec<char>) =
2937                        (short.chars().collect(), long.chars().collect());
2938                    let d = naive(&sc, &lc);
2939                    assert_eq!(levenshtein(&short, &long), d, "m={m} n={n}");
2940                    assert_eq!(levenshtein(&long, &short), d, "m={m} n={n} swapped");
2941                    let (sa, la) = (astral(&short), astral(&long));
2942                    assert_eq!(levenshtein(&sa, &la), d, "astral m={m} n={n}");
2943                }
2944            }
2945        }
2946    }
2947
2948    /// The public API, which transcodes each `&str` into the narrowest
2949    /// encoding that holds it, matches the oracles for every pair of
2950    /// encodings and every length regime.
2951    #[test]
2952    fn public_api_matches_oracle() {
2953        let alphabets: [&[char]; 4] = [
2954            &['a', 'b', 'c', 'd'],
2955            &['a', 'é', 'ß', 'ü'],
2956            &['a', 'é', 'Ω', '日'],
2957            &['a', 'é', '日', '😀'],
2958        ];
2959        let mut state = 0x5eed_u64;
2960        for a1 in alphabets {
2961            for a2 in alphabets {
2962                for &(m, n) in &[(0, 0), (0, 5), (3, 7), (40, 45), (70, 90), (600, 640)] {
2963                    let s1: String = (0..m).map(|_| a1[lcg(&mut state) % a1.len()]).collect();
2964                    let s2: String = (0..n).map(|_| a2[lcg(&mut state) % a2.len()]).collect();
2965                    let (c1, c2): (Vec<char>, Vec<char>) =
2966                        (s1.chars().collect(), s2.chars().collect());
2967                    assert_eq!(distance(&s1, &s2), naive(&c1, &c2), "{s1:?} {s2:?}");
2968                    let total = m + n;
2969                    let expected = if total == 0 {
2970                        1.0
2971                    } else {
2972                        1.0 - naive_indel(&c1, &c2) as f64 / total as f64
2973                    };
2974                    assert_eq!(ratio(&s1, &s2), expected, "{s1:?} {s2:?}");
2975                }
2976            }
2977        }
2978    }
2979
2980    #[test]
2981    fn hyrro_64_hash_matches_oracle() {
2982        let mut state = 0x1234_5678_9abc_def0u64;
2983        let alphabet: Vec<u32> = (0x3041..0x3090).collect();
2984        for &m in &[1, 2, 7, 16, 25, 31, 32, 33, 48, 63, 64] {
2985            for &n in &[m, m + 1, m + 5, 2 * m + 10] {
2986                let pat: Vec<u32> = (0..m)
2987                    .map(|_| {
2988                        state = state.wrapping_mul(6364136223846793005).wrapping_add(1);
2989                        alphabet[(state as usize) % alphabet.len()]
2990                    })
2991                    .collect();
2992                let text: Vec<u32> = (0..n)
2993                    .map(|_| {
2994                        state = state.wrapping_mul(6364136223846793005).wrapping_add(1);
2995                        alphabet[(state as usize) % alphabet.len()]
2996                    })
2997                    .collect();
2998                let pc: Vec<char> = pat.iter().map(|&c| char::from_u32(c).unwrap()).collect();
2999                let tc: Vec<char> = text.iter().map(|&c| char::from_u32(c).unwrap()).collect();
3000                let expected = naive(&pc, &tc);
3001                let got = hyrro_64_hash::<Lev, _, _>(&pat, &text);
3002                assert_eq!(got, expected, "m={m}, n={n}");
3003            }
3004        }
3005    }
3006}