Patterns in Permutations

Question

This fastest-code challenge is based partly on this MSE question and exists to extend some OEIS sequences, and create others. If I extend or create sequences based on this challenge, I'll link to this challenge so that folks can see where the values come from. If you'd like direct credit, include your first and last name so I can cite you on the entry.

These are now tables in the OEIS as A342840 (1342), A342860 (2413), A342861 (1324), A342862 (2143), A342863 (1243), A342864 (1432), and A342865 (1234).

Definition

Let \$\sigma \in S_n\$ be a permutation on \$n\$ letters, and let \$\tau \in S_m\$ be a permution on \$m \leq n\$ letters, both read as words: $$ \sigma = \sigma_1\sigma_2 \dots \sigma_n \text{ and } \tau = \tau_1\tau_2 \dots \tau_m $$ The permutation statistic \$\operatorname{pat}_\tau(\sigma)\$ counts the number of length-\$m\$ subsequences of \$\sigma\$ that are in the same relative order as \$\tau\$.

Challenge

This fastest-code challenge will have you take a permutation \$\tau \in S_m\$ and an integer \$n \geq m\$ and return a (0-indexed) list where position \$i\$ has the number of permutations \$\sigma \in S_n\$ such that \$\operatorname{pat}_{\tau}(\sigma) = i\$. If done correctly, the list will sum to \$n!\$.

Scoring

I will run your code on each of the following for one minute each for increasing values of \$n\$:

the only pattern in \$S_1\$: \$\tau = 1\$;
essentially the only pattern in \$S_2\$: \$\tau = 12\$;
the two essentially different patterns in \$S_3\$: \$\tau \in \{123, 132\}\$; and
the seven essentially distinct patterns of \$S_4\$: \$ \tau \in \{1234, 1243, 1324, 1342, 1432, 2143, 2413\}\$.

Your score will be the largest value of \$n\$ that can be computed within a minute for every pattern. For example, if you can compute up to

\$n = 1000\$ for \$\tau = 123\$ within a minute,
\$n = 10\$ for \$\tau = 1342\$ within a minute, and
\$n = 11\$ for \$\tau = 2413\$ within a minute,

then your score would be \$10\$.

Examples

Suppose that \$\sigma = 7532146\$ and \$\tau = 4213\$. Then $$ \operatorname{pat}_{4213}(7532146) = 13 $$ because \$\sigma = 7532146\$ has thirteen subsequences that are in the same relative order as \$\tau = 4213\$, \$7536\$, \$7526\$, \$7516\$, \$7546\$, \$7324\$, \$7326\$, \$7314\$, \$7316\$, \$7214\$, \$7216\$, \$5324\$, \$5314\$, and \$5214\$.

Test Data

 n | tau   | target output
---+-------+--------------------------------------------------------------
 6 |  1234 | [513, 102, 63, 10,  6, 12, 8, 0, 0,      5, 0, 0, 0, 0, 0, 1]
 5 |   123 | [ 42,  27, 24,  7,  9,  6, 0, 4, 0,      0, 1]
 9 |     1 | [  0,   0,  0,  0,  0,  0, 0, 0, 0, 362880]
 5 |   132 | [ 42,  21, 23, 14, 12,  5, 3]
 6 |  4213 | [512,  77, 69, 30, 21,  5, 6]
 4 |    12 | [  1,   3,  5,  6,  5,  3, 1]
 6 | 53214 | [694,  18,  7,  0,  1]

Note: For \$\tau = 12\$, the target output should be the \$n\$-th row of OEIS sequence A008302. For \$\tau = 1\$, the target output should be \$n\$ \$0\$s followed by \$n!\$.

I didn't quite understand what "subsequences that are in the same relative order as tau" means here. Does it mean the subsequence contains tau as its subsequence? Or any two letters in the subsequence will either 1. one or two of them not in tau, or 2. both them are in tau, and the order of these two letters is as same as in both subsequence or tau? — tsh, Mar 25, 2021 at 8:05
What are \$S_n\$ and \$S_m\$? Should I understand that they are two sets, and whenever \$m\le n\$, we have \$S_m \subseteq S_n\$? — tsh, Mar 25, 2021 at 8:07
\$S_n\$ is just the set of a all permutations of \$\{1,2,3,\dots,n\}\$. Strictly speaking, \$S_m\$ is not a subset of this unless \$m=n\$. But I suspect I may be misunderstanding what you’re getting at. — Peter Kagey, Mar 25, 2021 at 17:04
@tsh, "Subsequences that are in the same relative order as tau" means that if \$\tau\$ is $4123$, then we're looking for length-4 subsequences where the first number is the biggest, the second is the smallest, the third is the second-smallest, and the last number is the second-largest. — Peter Kagey, Mar 25, 2021 at 20:05
\$1342\$ and \$1423\$ are not essentially distinct, since they’re inverses. — Anders Kaseorg, Mar 28, 2021 at 4:28

Anders Kaseorg · Accepted Answer · 2021-03-28 16:59:26Z

Rust, score ≈ 13

This enumerates permutations σ recursively, sharing as much computation as possible at branches of the search tree, so that only constant work is needed at each leaf.

Times in seconds measured on my Ryzen 7 1800X (8 cores/16 threads):

τ   n = 10    11    12     13      14       15
──────────────────────────────────────────────
1     0.00  0.02  0.25   2.93   38.72   660.86
12    0.00  0.04  0.47   5.89   79.08  1268.18
123   0.00  0.06  0.67   8.47  115.61  1836.36
132   0.01  0.10  1.18  15.50  218.39  3443.56
1234  0.00  0.05  0.59   7.65  104.88  1643.29
1243  0.00  0.06  0.81  10.72  152.60  2417.75
1324  0.01  0.09  1.10  14.52  208.59  3280.18
1342  0.01  0.10  1.21  16.59  242.01  3871.64
1432  0.01  0.09  1.19  16.21  237.93  3808.25
2143  0.00  0.06  0.81  10.72  151.09  2409.76
2413  0.01  0.09  1.20  16.07  232.74  3717.70
──────────────────────────────────────────────
max   0.01  0.10  1.21  16.59  242.01  3871.64

Build and run with:

cargo build --release
target/release/permutations 11 '1 3 4 2'

Or try it online! (But note, the TIO version lacks parallelism and many other optimizations, so it’s much slower, due to being unable to use external crates.)

All outputs for n = 0…15.

Cargo.toml

[package]
name = "permutations"
version = "0.1.0"
authors = ["Anders Kaseorg <[email protected]>"]
edition = "2018"

[dependencies]
rayon = "1.5.0"
thread_local = "1.1.3"
typed-arena = "2.0.1"

[profile.release]
lto = "y"
panic = "abort"

src/main.rs

use rayon::prelude::*;
use std::cell::Cell;
use std::env;
use std::iter::Zip;
use std::slice;
use thread_local::ThreadLocal;
use typed_arena::Arena;

enum UpdateIterator<'a> {
    Leave(u16, slice::Iter<'a, (u16, u16)>),
    Take(
        u16,
        u16,
        Zip<slice::Iter<'a, (u16, u16)>, slice::Iter<'a, u16>>,
    ),
}

impl Iterator for UpdateIterator<'_> {
    type Item = (u16, u16);

    fn next(&mut self) -> Option<(u16, u16)> {
        match *self {
            UpdateIterator::Leave(i, ref mut bounds) => bounds
                .next()
                .map(|&(low, high)| (low + (low > i) as u16, high + (high >= i) as u16)),
            UpdateIterator::Take(i, tau0, ref mut zip) => {
                zip.next().map(|(&(low, high), &tau1)| {
                    if tau0 < tau1 {
                        (low.max(i) + 1, high + (high >= i) as u16)
                    } else {
                        (low + (low > i) as u16, high.min(i))
                    }
                })
            }
        }
    }

    fn size_hint(&self) -> (usize, Option<usize>) {
        match self {
            UpdateIterator::Leave(_, bounds) => bounds.size_hint(),
            UpdateIterator::Take(_, _, zip) => zip.size_hint(),
        }
    }
}

fn update(
    tau: &[u16],
    remaining: usize,
    states: &[(&[(u16, u16)], isize)],
    i: u16,
    mut callback: impl FnMut(UpdateIterator, isize),
) {
    for &(bounds, count) in states {
        if bounds.len() < remaining {
            callback(UpdateIterator::Leave(i, bounds.iter()), count);
        }

        if let Some((&(low0, high0), bounds1)) = bounds.split_first() {
            if low0 <= i && i <= high0 {
                callback(
                    UpdateIterator::Take(
                        i,
                        tau[tau.len() - bounds.len()],
                        bounds1.iter().zip(&tau[tau.len() - bounds.len() + 1..]),
                    ),
                    count,
                );
            }
        }
    }
}

fn search_step(
    tau: &[u16],
    remaining: usize,
    states: &[(&[(u16, u16)], isize)],
    arena_len: usize,
    i: u16,
    search: impl FnOnce(&[(&[(u16, u16)], isize)], usize),
) {
    let arena = Arena::with_capacity(2 * arena_len + 1 - states.len());
    let mut new_states = Vec::with_capacity(states.len() * 2);
    update(tau, remaining, states, i, |bounds, count| {
        new_states.push((&*arena.alloc_extend(bounds), count))
    });
    new_states.sort_by_key(|&(bounds, _)| bounds);
    new_states.dedup_by(|(bounds0, count0), (bounds1, count1)| {
        bounds0 == bounds1 && {
            *count1 += *count0;
            true
        }
    });
    search(&new_states, arena.len());
}

fn search(
    n: usize,
    tau: &[u16],
    k: usize,
    states: &[(&[(u16, u16)], isize)],
    arena_len: usize,
    output: &[Cell<u64>],
) {
    if n - k == 1 {
        let mut deltas = vec![0; n + 1];
        let mut total = 0;
        for &(bounds, count) in states {
            match bounds {
                [] => total += count,
                &[(low, high)] => {
                    deltas[low as usize] += count;
                    deltas[high as usize + 1] -= count;
                }
                _ => {}
            }
        }
        for delta in &deltas[..n] {
            total += delta;
            output[total as usize].set(output[total as usize].get() + 1);
        }
    } else if n - k == 2 {
        for i in 0..=k as u16 {
            let mut deltas = vec![0; n + 1];
            let mut total = 0;
            update(tau, n - k, states, i, |mut bounds, count| {
                if let Some((low, high)) = bounds.next() {
                    deltas[low as usize] += count;
                    deltas[high as usize + 1] -= count;
                    debug_assert!(bounds.next().is_none());
                } else {
                    total += count;
                }
            });
            for delta in &deltas[..n] {
                total += delta;
                output[total as usize].set(output[total as usize].get() + 1);
            }
        }
    } else {
        for i in 0..=k as u16 {
            search_step(tau, n - k, states, arena_len, i, |states, arena_len| {
                search(n, tau, k + 1, states, arena_len, output)
            })
        }
    }
}

fn par_search(
    n: usize,
    tau: &[u16],
    k: usize,
    states: &[(&[(u16, u16)], isize)],
    arena_len: usize,
    comb: usize,
    outputs: &ThreadLocal<Vec<Cell<u64>>>,
) {
    if n - k <= 4 {
        search(
            n,
            tau,
            k,
            states,
            arena_len,
            outputs.get_or(|| vec![Cell::new(0); comb + 1]),
        );
    } else {
        (0..=k as u16).into_par_iter().for_each(|i| {
            search_step(tau, n - k, states, arena_len, i, |states, arena_len| {
                par_search(n, tau, k + 1, states, arena_len, comb, outputs)
            })
        });
    }
}

fn counts(n: usize, tau: &[u16]) -> Vec<u64> {
    if n > tau.len() {
        let mut comb = 1;
        for i in 0..tau.len() {
            comb = comb * (n - i) / (i + 1);
        }
        let outputs = ThreadLocal::new();
        par_search(
            n,
            &tau,
            0,
            &[(&vec![(0, 0); tau.len()], 1)],
            tau.len(),
            comb,
            &outputs,
        );
        let mut output = vec![0; comb + 1];
        for thread_output in outputs {
            for (x, y) in output.iter_mut().zip(thread_output) {
                *x += y.get();
            }
        }
        output.truncate(
            output
                .iter()
                .rposition(|&count| count != 0)
                .map_or(0, |i| i + 1),
        );
        output
    } else if n == tau.len() {
        vec![(1..=n as u64).product::<u64>() - 1, 1]
    } else {
        vec![(1..=n as u64).product()]
    }
}

fn main() {
    let args: Vec<String> = env::args().collect();
    if let [_, n, tau] = &*args {
        let n: usize = n.parse().unwrap();
        let tau: Vec<u16> = tau.split_whitespace().map(|n| n.parse().unwrap()).collect();
        for count in counts(n, &tau) {
            println!("{}", count);
        }
    } else {
        panic!("usage: permutations n 'τ₁ τ₂ … τₘ'")
    }
}

Arnauld · Accepted Answer · 2021-03-25 01:19:27Z

JavaScript (Node.js), Score: 10 (estimated)

For now, this is a very simple implementation to make sure I understood the challenge correctly. The score is probably just 10, which is quite bad.

function count(a, b, c = [], p = 0, i = 0) {
  return (
    p == b.length ? 1 :
    i == a.length ? 0 :
    (
      b.every((v, k) => k >= p || Math.sign(b[p] - v) == Math.sign(a[i] - c[k])) &&
      count(a, b, [...c, a[i]], p + 1, i + 1)
    ) +
    count(a, b, c, p, i + 1)
  );
}

function perm(a) {
  let r = [];

  (function P(a, p = []) {
    if(a.length) {
      a.forEach((v, i) => P(a.filter(_ => i--), [...p, v]));
    }
    else {
      r.push(p);
    }
  })(a);

  return r;
}

function f(n, t) {
  let list = [];

  perm([...Array(n)].map((_, i) => i + 1)).forEach(p => {
    let x = count(p, t);

    if(!list[x]) {
      for(let i = 0; i <= x; i++) list[i] |= 0;
    }
    list[x]++;
  });
  return list;
}

Try it online!

Noodle9 · Accepted Answer · 2021-03-26 20:54:33Z

C++ (clang), score 11

#include "patterns_in_permutations.hpp"

int main(int argc, char *argv[])
{
    unsigned int running_time = 60;
    unsigned int n = 4;
    std::unordered_map<size_t, std::vector<std::string>> all_gens;
    const auto start = std::chrono::steady_clock::now();
    for (const auto finish = start + std::chrono::seconds{running_time}; n < 12 && std::chrono::steady_clock::now() < finish; ++n)
    {
        Perms p(n);
        const std::vector<std::string> gens = p.Generate();
        all_gens[n] = gens;
        Runner1(n);
    }
    const std::chrono::duration<double> duration = std::chrono::steady_clock::now() - start;
    std::cout << "It took " << duration.count() << " seconds.\n";
    Runner(Runner2, running_time, 4, "12", all_gens);
    Runner(Runner3, running_time, 4, "123", all_gens);
    Runner(Runner3, running_time, 4, "132", all_gens);
    Runner(Runner4, running_time, 4, "1234", all_gens);
    Runner(Runner4, running_time, 4, "1243", all_gens);
    Runner(Runner4, running_time, 4, "1324", all_gens);
    Runner(Runner4, running_time, 4, "1342", all_gens);
    Runner(Runner4, running_time, 4, "1423", all_gens);
    Runner(Runner4, running_time, 4, "1432", all_gens);
    Runner(Runner4, running_time, 4, "2143", all_gens);
    Runner(Runner4, running_time, 4, "2413", all_gens);

    return 0;
}