Can you calculate the average Levenshtein distance exactly? The Levenshtein distance between two strings is the minimum number of single character insertions, deletions, or substitutions to convert one string into the other one.

The challenge is to compute the average Levenshtein distance between two independent and uniformly random chosen binary strings of length n each.  Your output must be exact but can be given in any easy to understand human readable form. 

Examples:

These are the answer for n up to 24.

1 1/2 
2 1 
3 47/32 
4 243/128 
5 1179/512
6 2755/1024 
7 12561/4096 
8 56261/16384
9 124329/32768 
10 2175407/524288 
11 589839/131072
12 40664257/8388608 
13 174219279/33554432 
14 742795299/134217728
15 1576845897/268435456
16 13340661075/2147483648
17 14062798725/2147483648
18 59125997473/8589934592
19 123976260203/17179869184
20 259354089603/34359738368
21 8662782598909/1099511627776
22 72199426617073/8796093022208
23 150173613383989/17592186044416
24 1247439983177201/140737488355328


Score

Your score is the highest value of  you can reach. Where humanly possible, I will run your code on my Linux machine for 10 minutes and then kill the job to get the score.

Notes

As always this should be a competition per language. I will maintain a leaderboard that shows the best score for each language used in an answer.  I will also give a bounty of 50 points for the first answer to get to n = 20. 

My CPU is an Intel(R) Xeon(R) CPU X5460.

Leaderboard


n = 18 in Python+numba by Shamis (timing pending...).
n = 19 in Java by Bob Genom (278 seconds).
n = 19 in C by ngn (257 seconds).
n = 21 in Rust by Anders Kaseorg (297 seconds). 150 point bonus awarded.
Rust, score ≈ 22

This uses a dynamic programming approach (I’ve added an explanation here) 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(&amp;self, other: &amp;Distances) -&gt; 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&lt;H: Hasher&gt;(&amp;self, state: &amp;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 &gt; 31
            || (n as Count).leading_zeros() &lt; n as u32
            || (n as Total).leading_zeros() &lt; 2 * n as u32
        {
            break;
        }

        let total: Total = (0u32..1 &lt;&lt; (n - 1))
            .into_par_iter()
            .map(|a| {
                let mut a_sym = a.reverse_bits();
                a_sym ^= (a_sym &gt;&gt; 31).wrapping_neg();
                a_sym &gt;&gt;= 32 - n as usize;
                if a_sym &lt; a {
                    return 0;
                }

                let arena = Arena::&lt;Distances&gt;::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 &gt;&gt; i &amp; 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: &amp;[Distances] = arena.alloc_extend(
                    worst
                        .into_iter()
                        .map(both)
                        .collect::&lt;Box&lt;[Distance]&gt;&gt;()
                        .chunks(16)
                        .map(|chunk| {
                            Distances(unsafe {
                                _mm_loadu_si128(chunk as *const [i8] as *const __m128i)
                            })
                        }),
                );

                let mut states: HashMap&lt;&amp;[Distances], Count, FxBuildHasher&gt; = 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&lt;_&gt; = [a, !a]
                    .iter()
                    .map(|b| {
                        arena.alloc_extend((0..n + 1).step_by(16).map(|i| unsafe {
                            let x = _mm_set1_epi16(((b &lt;&lt; 1) &gt;&gt; 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 = &amp;worst[idx(0, j) / 16..idx(0, j + 1) / 16];
                    for (state, count) in states.drain() {
                        for bv in &amp;bvs {
                            let mut x = j;
                            let mut y = n.into();
                            let mut bound = n;

                            let new_state: &amp;mut [Distances] =
                                arena.alloc_extend(izip!(&amp;**bv, state, worst_slice).map(
                                    |(&amp;Distances(bc), &amp;Distances(yc), &amp;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(&amp;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(&amp;*new_state).or_insert(0) += count;
                        }
                    }
                    mem::swap(&amp;mut states, &amp;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::&lt;Count&gt;(),
                ) * if a_sym == a { 1 } else { 2 }
            })
            .sum();

        let shift = total.trailing_zeros();
        println!(
            "{} {}/{}",
            n,
            total &gt;&gt; shift,
            (1 as Total) &lt;&lt; (2 * n as u32 - 1 - shift),
        );
    }
}


Cargo.toml

[package]
name = "levenshtein"
version = "0.1.0"
authors = ["Anders Kaseorg &lt;andersk@mit.edu&gt;"]
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