Rust, score ≈ 22
This uses a dynamic programming approach whose running time seems to scale as approximately \$\tilde O(2^{1.5n})\$, rather than the \$\tilde O(2^{2n})\$ of a brute-force search. On my Ryzen 7 1800X (8 cores/16 threads), it gets through \$1 \le n \le 21\$ in 1.7 minutes, \$1 \le n \le 22\$ in 5.1 minutes.
Now using SIMD for the inner loop.
src/main.rs
use fxhash::FxBuildHasher;
use itertools::izip;
use rayon::prelude::*;
use std::arch::x86_64::*;
use std::collections::HashMap;
use std::hash::{Hash, Hasher};
use std::mem;
use typed_arena::Arena;
#[global_allocator]
static ALLOC: mimallocator::Mimalloc = mimallocator::Mimalloc;
type Distance = i8;
type Count = u32;
type Total = u64;
#[derive(Debug)]
struct Distances(__m128i);
impl PartialEq for Distances {
fn eq(&self, other: &Distances) -> bool {
unsafe {
let x = _mm_xor_si128(self.0, other.0);
_mm_testz_si128(x, x) != 0
}
}
}
impl Eq for Distances {}
impl Hash for Distances {
fn hash<H: Hasher>(&self, state: &mut H) {
unsafe {
_mm_extract_epi64(self.0, 0).hash(state);
_mm_extract_epi64(self.0, 1).hash(state);
}
}
}
fn main() {
let splat0 = unsafe { _mm_set1_epi8(0) };
let splat1 = unsafe { _mm_set1_epi8(1) };
let splatff = unsafe { _mm_set1_epi8(!0) };
let splat7f = unsafe { _mm_set1_epi8(0x7f) };
let seq = unsafe { _mm_set_epi8(15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0) };
let grow0 = unsafe {
_mm_set_epi8(
-0x80, 0, 0x20, 0, 0x8, 0, 0x2, 0, 0, 0x40, 0, 0x10, 0, 0x4, 0, 0x1,
)
};
let grow1 = unsafe {
_mm_set_epi8(
0, 0x40, 0, 0x10, 0, 0x4, 0, 0x1, -0x80, 0, 0x20, 0, 0x8, 0, 0x2, 0,
)
};
for n in 1 as Distance.. {
if n > 31
|| (n as Count).leading_zeros() < n as u32
|| (n as Total).leading_zeros() < 2 * n as u32
{
break;
}
let total: Total = (0u32..1 << (n - 1))
.into_par_iter()
.map(|a| {
let mut a_sym = a.reverse_bits();
a_sym ^= (a_sym >> 31).wrapping_neg();
a_sym >>= 32 - n as usize;
if a_sym < a {
return 0;
}
let arena = Arena::<Distances>::new();
let stride = (n as usize + 16) / 16 * 16;
let idx = |i: Distance, j: Distance| i as usize + stride * j as usize;
let both = |[x, y]: [Distance; 2]| x.max(y);
let mut worst = vec![[stride as Distance; 2]; idx(0, n + 1)];
for j in 0..=n {
worst[idx(n, j)] = [n - j; 2];
}
for i in (0..n).rev() {
worst[idx(i, n)] = [n - i; 2];
let good = (a >> i & 1) as usize;
let bad = good ^ 1;
for j in (0..n).rev() {
worst[idx(i, j)][good] = both(worst[idx(i + 1, j + 1)]);
worst[idx(i, j)][bad] = 1 + worst[idx(i + 1, j)][bad]
.min(both(worst[idx(i, j + 1)]))
.min(both(worst[idx(i + 1, j + 1)]));
}
}
let worst: &[Distances] = arena.alloc_extend(
worst
.into_iter()
.map(both)
.collect::<Box<[Distance]>>()
.chunks(16)
.map(|chunk| {
Distances(unsafe {
_mm_loadu_si128(chunk as *const [i8] as *const __m128i)
})
}),
);
let mut states: HashMap<&[Distances], Count, FxBuildHasher> = HashMap::default();
let mut new_states = HashMap::default();
states.insert(
arena.alloc_extend(
(0..n + 1)
.step_by(16)
.map(|i| Distances(unsafe { _mm_add_epi8(_mm_set1_epi8(i), seq) })),
),
1,
);
let bvs: Vec<_> = [a, !a]
.iter()
.map(|b| {
arena.alloc_extend((0..n + 1).step_by(16).map(|i| unsafe {
let x = _mm_set1_epi16(((b << 1) >> i) as i16);
Distances(_mm_xor_si128(
_mm_cmpeq_epi8(
_mm_or_si128(
_mm_and_si128(x, grow0),
_mm_and_si128(_mm_alignr_epi8(x, x, 1), grow1),
),
splat0,
),
splatff,
))
}))
})
.collect();
for j in 1..=n {
new_states.reserve(2 * states.len());
let worst_slice = &worst[idx(0, j) / 16..idx(0, j + 1) / 16];
for (state, count) in states.drain() {
for bv in &bvs {
let mut x = j;
let mut y = n.into();
let mut bound = n;
let new_state: &mut [Distances] =
arena.alloc_extend(izip!(&**bv, state, worst_slice).map(
|(&Distances(bc), &Distances(yc), &Distances(wc))| unsafe {
let o = _mm_min_epi8(
_mm_add_epi8(yc, splat1),
_mm_sub_epi8(
_mm_insert_epi8(_mm_slli_si128(yc, 1), y, 0),
bc,
),
);
y = _mm_extract_epi8(yc, 15);
let o = _mm_sub_epi8(o, seq);
let o = _mm_min_epi8(o, _mm_set1_epi8(x));
let o = _mm_sub_epi8(splat7f, o);
let o = _mm_max_epu8(o, _mm_slli_si128(o, 1));
let o = _mm_max_epu8(o, _mm_slli_si128(o, 2));
let o = _mm_max_epu8(o, _mm_slli_si128(o, 4));
let o = _mm_max_epu8(o, _mm_slli_si128(o, 8));
let o = _mm_sub_epi8(splat7f, o);
x = _mm_extract_epi8(o, 15) as i8 + 16;
let o = _mm_add_epi8(o, seq);
let z = _mm_add_epi8(o, wc);
let z = _mm_min_epi8(z, _mm_srli_si128(z, 1));
let z = _mm_min_epi8(z, _mm_srli_si128(z, 2));
let z = _mm_min_epi8(z, _mm_srli_si128(z, 4));
let z = _mm_min_epi8(z, _mm_srli_si128(z, 8));
bound = bound.min(_mm_extract_epi8(z, 0) as i8);
Distances(o)
},
));
let bound = unsafe { _mm_set1_epi8(bound) };
for (i, Distances(x)) in (0..).step_by(16).zip(&mut *new_state) {
*x = unsafe {
_mm_min_epi8(
*x,
_mm_sub_epi8(
bound,
_mm_abs_epi8(_mm_add_epi8(_mm_set1_epi8(i - j), seq)),
),
)
};
}
*new_states.entry(&*new_state).or_insert(0) += count;
}
}
mem::swap(&mut states, &mut new_states);
}
let control = unsafe { _mm_insert_epi8(splatff, (n % 16).into(), 0) };
Total::from(
states
.into_iter()
.map(|(state, count)| unsafe {
count
* _mm_extract_epi8(
_mm_shuffle_epi8(state[n as usize / 16].0, control),
0,
) as Count
})
.sum::<Count>(),
) * if a_sym == a { 1 } else { 2 }
})
.sum();
let shift = total.trailing_zeros();
println!(
"{} {}/{}",
n,
total >> shift,
(1 as Total) << (2 * n as u32 - 1 - shift),
);
}
}
Cargo.toml
[package]
name = "levenshtein"
version = "0.1.0"
authors = ["Anders Kaseorg <andersk@mit.edu>"]
edition = "2018"
[profile.release]
lto = true
codegen-units = 1
[dependencies]
fxhash = "0.2.1"
itertools = "0.8.2"
mimallocator = "0.1.3"
rayon = "1.3.0"
typed-arena = "2.0.0"
Running
RUSTFLAGS='-C target-cpu=native' cargo build --release
target/release/levenshtein
Output
(With cumulative timing data prepended by ts -s %.s.)
0.000008 1 1/2
0.000150 2 1/1
0.000219 3 47/32
0.000282 4 243/128
0.000344 5 1179/512
0.000413 6 2755/1024
0.000476 7 12561/4096
0.000538 8 56261/16384
0.000598 9 124329/32768
0.000660 10 2175407/524288
0.000721 11 589839/131072
0.000782 12 40664257/8388608
0.000843 13 174219279/33554432
0.006964 14 742795299/134217728
0.068070 15 1576845897/268435456
0.310136 16 13340661075/2147483648
1.062122 17 14062798725/2147483648
3.586745 18 59125997473/8589934592
11.265840 19 123976260203/17179869184
33.691822 20 259354089603/34359738368
101.514674 21 8662782598909/1099511627776
307.427106 22 72199426617073/8796093022208
956.299101 23 150173613383989/17592186044416
3077.477731 24 1247439983177201/140737488355328
10276.205241 25 5173410986415247/562949953421312
34550.754308 26 5356540527479769/562949953421312
Static core2 build for Anush
