text:
Find the number of n-by-n (-1, 0, 1) matrices with zero permanent as quickly as possible
The permanent of an \$n\$-by-\$n\$ matrix \$A = (a_{i,j})\$ is defined as:
$$\operatorname{perm}(A)=\sum_{\sigma\in S_n}\prod_{i=1}^n a_{i,\sigma(i)}$$
For a fixed \$n\$, consider the \$n\$-by-\$n\$ matrices whose entries are taken from \$\{-1, 0, +1\}\$ .
Task
For each \$n\$ from 1 upwards, output the number of \$n\$-by-\$n\$ \$(-1, 0, 1)\$ matrices with zero permanent.
The output format looks like this:
n = 1, 1
n = 2, 33
n = 3, 7555
n = 4, 13482049
n = 5, 186481694371
n = 6, ...

Score
Your score is the largest \$n\$ your code gets to in 20 minutes in total on my computer, this is not 20 minutes per \$n\$.
If two entries get the same \$n\$ score, then the winning entry will be the one that gets to the highest \$n\$ in the shortest time on my machine. If the two best entries are equal on this criterion too then the winner will be the answer submitted first.
You are in no way allowed to hard-code solutions.
Languages and libraries
You can use any freely available language and libraries you like. I must be able to run your code so please include a full explanation for how to run/compile your code in windows if at all possible.
My Machine
The timings will be run on my machine, Intel Core i7-9750H 6 cores, 12 threads, The machine runs Windows 10.
Updated
Brute forcing every matrix (of which there are \$3^{n^2}\$) is not necessarily a particularly good answer, especially since there are many symmetries that preserve the value of the permanent. Thanks for @isaacg.

We have MathJax, just use \$ for the delimiter.

Possible duplicate of this post and this post.

It's only a duplicate if you can't think of a better way than brute force. This is fastest-algorithm, and surely it can be done better than by checking all matrices!

This is definitely not a duplicate - brute forcing every matrix (of which there are 3^(n^2)) is not necessarily a particularly good answer, especially since there are many symmetries that preserve the value of the permanent. Exploiting features like that will be necessary to do well, so this is not a duplicate.

I'm cautiously reopening this based on what isaacg points out about counting vs computing. I ask that people not post solutions basically of the form "generate all matrices and count their permanents" because those are unlikely to be competitive and should go on the compute-permanent challenge.

If there will be a non brute force answer then indeed this will be both a very interesting question & answer.

Some brute force results by Mathematica: 2 -> 33, 3 -> 7555, 4 -> no longer brute-forceable.

FWIW, I have $13482049$ for $n=4$.

According to this OEIS draft of A330134, the answer is $186481694371$ for $n=5$. (The sequence will be here once approved.)

Rust, \$A(5) = 186481694371\$ in 0.2 s, \$A(6) = 19733690332538577\$ in 580 s
(Unofficial times on a Core i7-10710U with 6 cores/12 threads.)
src/main.rs
use itertools::Itertools;
use num_integer::gcd;
use rayon::prelude::*;
use std::cmp::Ordering;
use std::collections::HashMap;
use std::convert::TryFrom;
use std::sync::Mutex;

type Count = u64;
type Coeff = i32;

fn sort_scols(n: u32, i: u32, sets: &[u32], state: &[Coeff]) -> Vec<u32> {
    let mut scol_ranks: Vec<u32> = (0..2 * n).map(|scol| (scol >> 1 > i) as u32).collect();
    let mut sorted_scols: Vec<u32> = (0..2 * n).collect();
    let mut scol_descs = vec![vec![]; 2 * n as usize];
    for (set_index, (&set, &coeff)) in sets.iter().zip(state).enumerate() {
        if coeff != 0 {
            for col in 0..n {
                if set & 1 << col != 0 {
                    scol_descs[(col as usize) << 1].push(set_index);
                    scol_descs[(col as usize) << 1 | 1].push(set_index);
                }
            }
        }
    }
    let mut sorted_set_indexes: Vec<_> = (0..sets.len()).collect();
    let mut set_ranks: Vec<u32> = vec![!0; sets.len()];
    let mut set_descs = vec![(0, vec![]); sets.len()];

    while {
        for (&coeff, set_desc) in state.iter().zip(&mut set_descs) {
            set_desc.0 = coeff;
            set_desc.1.clear();
        }
        let mut seen = 0;
        for &scol in &sorted_scols {
            if seen & 1 << (scol >> 1) == 0 {
                seen |= 1 << (scol >> 1);
                if scol & 1 == 0 {
                    for (&set, set_desc) in sets.iter().zip(&mut set_descs) {
                        if set_desc.0 != 0 && set & 1 << (scol >> 1) != 0 {
                            set_desc.1.push(scol_ranks[scol as usize]);
                        }
                    }
                } else {
                    for (&set, set_desc) in sets.iter().zip(&mut set_descs) {
                        if set_desc.0 != 0 && set & 1 << (scol >> 1) != 0 {
                            set_desc.0 = -set_desc.0;
                            set_desc.1.push(scol_ranks[scol as usize]);
                        }
                    }
                }
            }
        }

        let set_index_key =
            |&set_index: &usize| (set_descs[set_index].0.abs(), &set_descs[set_index].1);
        sorted_set_indexes.sort_by_key(set_index_key);
        if sets.len() != 0 {
            let mut next_set_rank = 0;
            set_ranks[sorted_set_indexes[0]] = next_set_rank;
            for (&set_index0, &set_index1) in sorted_set_indexes.iter().tuple_windows() {
                if set_index_key(&set_index0) < set_index_key(&set_index1) {
                    next_set_rank += 1;
                }
                set_ranks[set_index1] = next_set_rank;
            }
        }

        let scol_set_index_key = |scol: u32| {
            let sign = 1 - 2 * (scol & 1) as Coeff;
            let set_ranks = &set_ranks[..];
            let set_descs = &set_descs[..];
            move |&set_index: &usize| (set_ranks[set_index], set_descs[set_index].0 * sign)
        };
        for (scol, scol_desc) in scol_descs.iter_mut().enumerate() {
            scol_desc.sort_by_key(scol_set_index_key(scol as u32));
        }
        let scol_cmp = |&scol0: &u32, &scol1: &u32| {
            scol_ranks[scol0 as usize]
                .cmp(&scol_ranks[scol1 as usize])
                .then_with(|| {
                    scol_descs[scol0 as usize]
                        .iter()
                        .map(scol_set_index_key(scol0))
                        .cmp(
                            scol_descs[scol1 as usize]
                                .iter()
                                .map(scol_set_index_key(scol1)),
                        )
                })
        };
        sorted_scols.sort_by(scol_cmp);

        let mut new_scol_ranks = vec![!0; 2 * n as usize];
        let mut next_rank = 0;
        let mut changed = false;
        new_scol_ranks[sorted_scols[0] as usize] = next_rank;
        for (&scol0, &scol1) in sorted_scols.iter().tuple_windows() {
            if scol_ranks[scol0 as usize] < scol_ranks[scol1 as usize] {
                next_rank += 1;
            } else if scol_cmp(&scol0, &scol1) == Ordering::Less {
                changed = true;
                next_rank += 1;
            }
            new_scol_ranks[scol1 as usize] = next_rank;
        }
        scol_ranks = new_scol_ranks;

        changed
    } {}

    let mut seen = 0;
    sorted_scols.retain(|scol| {
        seen & 1 << (scol >> 1) == 0 && {
            seen |= 1 << (scol >> 1);
            true
        }
    });
    assert_eq!(seen, !(!0 << n));

    sorted_scols
}

fn count(n: u32) -> Count {
    static ERR: &str = "n is too large";
    Count::checked_pow(3, TryFrom::try_from(n * n).expect(ERR)).expect(ERR);
    (1..=Coeff::try_from(n).expect(ERR)).fold(1, |x: Coeff, y| x.checked_mul(y).expect(ERR));

    let mut cur_sets = vec![0];
    let mut cur_set_indexes = vec![0];

    let mut cur_states = HashMap::new();
    cur_states.insert(vec![1], 1);

    for k in 0..n {
        for i in 0..n {
            let keep_set = |&set: &u32| {
                if set.count_ones() == k {
                    !set & !1 << i & !(!0 << n) != 0
                } else {
                    set & 1 << i == 0
                }
            };
            let enlarge_set = |&set: &u32| set.count_ones() == k && set & 1 << i == 0;
            let new_sets: Vec<u32> = cur_sets
                .iter()
                .copied()
                .filter(keep_set)
                .chain(
                    cur_sets
                        .iter()
                        .copied()
                        .filter(enlarge_set)
                        .map(|set| set | 1 << i),
                )
                .collect();
            let mut new_set_indexes = vec![!0; 1 << n];
            for (set_index, &set) in new_sets.iter().enumerate() {
                new_set_indexes[set as usize] = set_index;
            }

            let new_states = Mutex::new(HashMap::new());
            cur_states.into_par_iter().for_each(|(cur_state, count)| {
                for elt in -1..=1 {
                    let mut new_state: Vec<Coeff> = cur_sets
                        .iter()
                        .copied()
                        .zip(&cur_state)
                        .filter(|(set, _)| keep_set(set))
                        .map(|(_, &coeff)| coeff)
                        .chain(
                            cur_sets
                                .iter()
                                .copied()
                                .zip(&cur_state)
                                .filter(|(set, _)| enlarge_set(set))
                                .map(|(set, &coeff)| {
                                    let mut new_coeff = coeff * elt;
                                    if set & !(!0 << i) != 0 {
                                        new_coeff +=
                                            cur_state[cur_set_indexes[set as usize | 1 << i]];
                                    }
                                    new_coeff
                                }),
                        )
                        .collect();

                    let factor = new_state.iter().copied().fold(0, gcd);
                    if factor > 1 {
                        for coeff in &mut new_state {
                            *coeff /= factor;
                        }
                    }

                    let sorted_scols = sort_scols(n, i, &new_sets, &new_state);

                    let new_state: Vec<Coeff> = new_sets
                        .iter()
                        .map(|&set| {
                            let mut permuted_set = 0;
                            let mut sign = 0;
                            for (col, &permuted_scol) in sorted_scols.iter().enumerate() {
                                sign ^= set >> col & permuted_scol & 1;
                                permuted_set |= (set >> col & 1) << (permuted_scol >> 1);
                            }
                            (1 - 2 * sign as Coeff)
                                * new_state[new_set_indexes[permuted_set as usize] as usize]
                        })
                        .collect();

                    *new_states.lock().unwrap().entry(new_state).or_insert(0) += count;
                }
            });

            cur_sets = new_sets;
            cur_set_indexes = new_set_indexes;
            cur_states = new_states.into_inner().unwrap();
        }
    }

    *cur_states.get(&vec![0]).unwrap_or(&0)
}

fn main() {
    for n in 1..=6 {
        println!("{} {}", n, count(n));
    }
}

Cargo.toml
[package]
name = "permanent"
version = "0.1.0"
authors = ["Anders Kaseorg <andersk@mit.edu>"]
edition = "2018"

[dependencies]
num-integer = "0.1.41"
rayon = "1.2.1"
itertools = "0.8.2"

Building and running
$ cargo build --release
$ time target/release/permanent
1 1
2 33
3 7555
4 13482049
5 186481694371
6 19733690332538577
5995.30user 282.67system 9:42.97elapsed 1076%CPU (0avgtext+0avgdata 6479028maxresident)k
0inputs+0outputs (0major+3259222minor)pagefaults 0swaps

How it works
This counts matrices using dynamic programming so we don’t need to enumerate all of them directly.  Basically, the state of a partially filled matrix is represented by the coefficients of an expression for the permanent as a polynomial in the remaining unfilled elements.  Specifically, it corresponds to a partial execution of the following algorithm for calculating the permanent:
\$\operatorname{perm} A = c_{\{0, 1, \dotsc, n - 1\}}\$, where \$c_\varnothing = 1\$, \$c_S = \sum_{i \in S} c_{S \smallsetminus \{i\}}a_{|S| - 1, i}\$ for \$S \ne \varnothing\$.
For example, here’s the \$n = 4\$ calculation written out explicitly:
c = 1

c0 = c * a00
c1 = c * a01
c2 = c * a02
c3 = c * a03

c01  = c1 * a10; c02  = c2 * a10; c03  = c3 * a10
c01 += c0 * a11; c12  = c2 * a11; c13  = c3 * a11
c02 += c0 * a12; c12 += c1 * a12; c23  = c3 * a12
c03 += c0 * a13; c13 += c1 * a13; c23 += c2 * a13

c012  = c12 * a20; c013  = c13 * a20; c023  = c23 * a20
c012 += c02 * a21; c013 += c03 * a21; c123  = c23 * a21
c012 += c01 * a22; c023 += c03 * a22; c123 += c13 * a22
c013 += c01 * a23; c023 += c02 * a23; c123 += c12 * a23

c0123  = c123 * a30
c0123 += c023 * a31
c0123 += c013 * a32
c0123 += c012 * a33

return c0123

After executing each line, we can forget the filled elements a… and represent the current state with just the c… variables that are still live.  This representation collapses many different partially filled matrices into the same coefficient vector.  For example, swapping any two filled rows results in the same coefficient vector.
The sort_scols function further collapses the state space by using the symmetries of column permutations and column sign flipping to try to find a canonical representation for each state.  It doesn’t do a perfect job, but this imperfection only impacts performance, not correctness.

C# (\$\textrm{A330134}(5)=186481694371\$ in under 10 seconds)
using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;

namespace Sandbox
{
    // See https://codegolf.stackexchange.com/q/196712
    static class Permanent
    {
        const int numThreads = 12;

        static void Main(string[] args)
        {
            var sw = new System.Diagnostics.Stopwatch();
            sw.Start();
            for (int n = 1; n < 8; n++)
            {
                Console.WriteLine($"{n}\t{_PPCG196712_LaplaceMT(n)}\t{sw.ElapsedMilliseconds}ms");
            }
        }

        internal static BigInteger _PPCG196712_LaplaceMT(int n)
        {
            // If we have a zero row then the permanent is zero.
            // If the first zero row is number i then the previous i rows have (3^n)-1 possibilities each, and the remaining ones have 3^n.
            BigInteger result = 0;
            int threePowN = (int)BigInteger.Pow(3, n);
            for (int i = 0; i < n; i++)
            {
                result += BigInteger.Pow(threePowN - 1, i) * BigInteger.Pow(threePowN, n - i - 1);
            }

            // Let's call a row containing a permutation of 1^1 0^{n-1} a "singleton row".
            // Under Laplace expansion it simplifies to a single {n-1}x{n-1} permanent.
            // If we have two singleton rows with the same non-zero index, that {n-1}x{n-1} permanent has a zero row, so the nxn permanent must be zero.
            // If we have multiple distinct singleton rows, we still need to worry about double-counting.
            // But all this effort would give roughly a 1.5% improvement at n=6, so is probably not worth it.

            // Now the interesting rows. We cluster similar rows as a technique to optimise column permutations.
            var clusters = new List<Cluster>();
            for (int ones = 1; ones <= n; ones++)
                for (int negs = 0; negs <= ones && ones + negs <= n; negs++)
                    clusters.Add(new Cluster(n, ones, negs));

            // Consider partitions of n-1 by cluster.
            foreach (var partition in PartitionAssignments(clusters, n-1))
            {
                var subclusters = partition.OrderBy(tuple =>
                {
                    var (cluster, k) = tuple;
                    // What's the speedup factor of pivoting on this cluster?
                    return k < cluster.DistinctCounts.Count ? cluster.DistinctCounts[k] / (double)cluster.NaiveCounts[k] : 1;
                }).ToList();

                // We pivot on the first cluster, and take a Cartesian product.
                result += Enumerable.Range(0, numThreads).AsParallel().
                    Select(threadId => CountByPartition(n, subclusters, 0, new sbyte[n - 1][], 0, 1, threadId)).
                    Aggregate(BigInteger.Zero, (a, b) => a + b);
            }

            return result;
        }

        private static BigInteger CountByPartition(int n, IReadOnlyList<(Cluster cluster, int repeat)> clusters, int clusterIdx, sbyte[][] matrix, int rowIdx, BigInteger weight, int threadId)
        {
            if (clusterIdx == clusters.Count) return weight * Count_Laplace(matrix) << (n - 1);

            int k = clusters[clusterIdx].repeat;
            if (clusterIdx == 0 && clusters[0].repeat < clusters[0].cluster.DistinctCounts.Count)
            {
                // Pivot
                weight *= Binom(n - 1, k) * Factorial(n - 1 - k); // Choose the pivoted rows and then permute the others freely

                BigInteger sum = 0;
                var reps = clusters[0].cluster.Representatives[k];

                int delta = reps.Count >= numThreads || clusterIdx == clusters.Count - 1 ? numThreads : 1;
                int skip = delta == 1 ? 0 : threadId;

                for (int repIdx = skip; repIdx < reps.Count; repIdx += delta)
                {
                    var (prefixWeight, prefix) = reps[repIdx];
                    for (int i = 0; i < k; i++) matrix[rowIdx + i] = prefix[i];
                    sum += CountByPartition(n, clusters, clusterIdx + 1, matrix, rowIdx + k, weight * prefixWeight, delta == 1 ? threadId : -1);
                }
                return sum;
            }

            {
                if (clusterIdx == 0) weight = Factorial(n - 1); // Permute all rows freely

                var rows = clusters[clusterIdx].cluster.Rows;
                var indices = new int[k];
                BigInteger sum = 0;

                int skip = threadId > 0 ? threadId : 0;
                int delta = threadId >= 0 ? numThreads : 1;

                while (true)
                {
                    if (skip == 0)
                    {
                        for (int i = 0; i < k; i++)
                        {
                            var j = indices[i];
                            var inRow = rows[j];
                            matrix[rowIdx + i] = inRow;
                        }

                        // Adjust the weight to account for duplicates in this cluster.
                        var recurWeight = weight;
                        int currRun = 1;
                        for (int i = 1; i < k; i++)
                        {
                            if (indices[i] == indices[i - 1]) currRun++;
                            else { recurWeight /= Factorial(currRun); currRun = 1; }
                        }
                        recurWeight /= Factorial(currRun);

                        sum += CountByPartition(n, clusters, clusterIdx + 1, matrix, rowIdx + k, recurWeight, -1);
                    }

                    skip--;
                    if (skip < 0) skip += delta;

                    if (indices[0] == rows.Count - 1) break;
                    NextOrderedMultiset(indices, rows.Count);
                }

                return sum;
            }
        }

        private static IEnumerable<IReadOnlyList<(T, int)>> PartitionAssignments<T>(IReadOnlyList<T> elts, int n)
        {
            IEnumerable<IReadOnlyList<(T, int)>> Inner(int m, int maxPart, int off, List<(T, int)> prefix, ISet<int> indicesUsed)
            {
                if (m == 0)
                {
                    yield return new List<(T, int)>(prefix);
                    yield break;
                }

                for (int part = Math.Min(m, maxPart); part > 0; part--)
                {
                    // When we reduce the part, we re-enable skipped elts.
                    for (int i = part == maxPart ? off : 0; i < elts.Count; i++)
                    {
                        if (indicesUsed.Contains(i)) continue;

                        prefix.Add((elts[i], part));
                        indicesUsed.Add(i);

                        foreach (var soln in Inner(m - part, part, i + 1, prefix, indicesUsed)) yield return soln;

                        indicesUsed.Remove(i);
                        prefix.RemoveAt(prefix.Count - 1);
                    }
                }
            }

            return Inner(n, n, 0, new List<(T, int)>(), new HashSet<int>());
        }

        private static void NextOrderedMultiset(int[] indices, int @base)
        {
            int j = indices.Length - 1;
            while (indices[j] == @base - 1) j--;
            indices[j]++;
            j++;
            while (j < indices.Length) { indices[j] = indices[j - 1]; j++; }
        }

        private static BigInteger Count_Laplace(sbyte[][] matrix)
        {
            int m = matrix.Length;
            int n = m + 1;
            BigInteger[] sums = new BigInteger[m];
            uint gray = 0;
            int sign = 1;
            BigInteger[] weights = new BigInteger[n];
            for (uint i = 0; i < 1u << n; i++)
            {
                BigInteger term = 1;
                foreach (var sum in sums) term *= sum;

                for (int j = 0; j < n; j++)
                {
                    if (((gray >> j) & 1) == 1) weights[j] += sign * term;
                }

                // Gray code update
                int flipPos = (i + 1).CountTrailingZeros();
                if (flipPos == n) break;

                uint bit = 1u << flipPos;
                if ((gray & bit) == bit)
                {
                    for (int j = 0; j < m; j++) sums[j] -= matrix[j][flipPos];
                }
                else
                {
                    for (int j = 0; j < m; j++) sums[j] += matrix[j][flipPos];
                }
                gray ^= bit;
                sign = -sign;
            }

            // Count weighted subsets which total zero, optimising with meet-in-the-middle
            var distribL = SubsetSums(weights, 0, n >> 1);
            var distribR = SubsetSums(weights, n >> 1, n - (n >> 1));
            BigInteger total = -1; // Discount the zero column solution
            foreach (var kvp in distribL)
            {
                if (distribR.TryGetValue(-kvp.Key, out var v2)) total += kvp.Value * v2;
            }

            return total;
        }

        private static Dictionary<BigInteger, BigInteger> SubsetSums(BigInteger[] weights, int off, int len)
        {
            var distrib = new Dictionary<BigInteger, BigInteger> { [0] = 1 };
            for (int i = 0; i < len; i++)
            {
                var weight = weights[off + i];
                var nextDistrib = new Dictionary<BigInteger, BigInteger>(distrib);
                foreach (var kvp in distrib)
                {
                    nextDistrib.TryGetValue(kvp.Key - weight, out var tmp);
                    nextDistrib[kvp.Key - weight] = tmp + kvp.Value;

                    nextDistrib.TryGetValue(kvp.Key + weight, out tmp);
                    nextDistrib[kvp.Key + weight] = tmp + kvp.Value;
                }

                distrib = nextDistrib;
            }

            return distrib;
        }

        class Cluster
        {
            internal readonly IReadOnlyList<sbyte[]> Rows;

            internal readonly IReadOnlyList<int> NaiveCounts;
            internal readonly IReadOnlyList<int> DistinctCounts;
            internal readonly IReadOnlyList<IReadOnlyList<(int weight, IReadOnlyList<sbyte[]> prefix)>> Representatives;

            internal Cluster(int n, int ones, int negs)
            {
                if (ones == 0) throw new ArgumentOutOfRangeException(nameof(ones));
                if (negs > ones) throw new ArgumentOutOfRangeException(nameof(negs));

                Rows = BuildRows(n, ones, negs);

                // Build pairs and triples, and calculate their advantage factors.
                var naiveCounts = new List<int> { 1, Rows.Count };
                var distinctCounts = new List<int> { 1, 1 };
                var representatives = new List<IReadOnlyList<(int, IReadOnlyList<sbyte[]>)>>
                {
                    new List<(int, IReadOnlyList<sbyte[]>)>{ (1, new List<sbyte[]>()) },
                    new List<(int, IReadOnlyList<sbyte[]>)>{ (Rows.Count, new List<sbyte[]> { Rows.First() }) }
                };

                for (int k = 2; k < 4; k++)
                {
                    var counts = new Dictionary<BigInteger, int>();
                    var reps = new Dictionary<BigInteger, IReadOnlyList<sbyte[]>>();
                    int uncollapsed = 0;

                    var indices = new int[k];
                    while (true)
                    {
                        var subset = new sbyte[k][];
                        for (int i = 0; i < k; i++) subset[i] = Rows[indices[i]];

                        BigInteger count = Factorial(k);
                        int currRun = 1;
                        for (int i = 1; i < k; i++)
                        {
                            if (indices[i] == indices[i - 1]) currRun++;
                            else { count /= Factorial(currRun); currRun = 1; }
                        }
                        count /= Factorial(currRun);

                        uncollapsed += (int)count;
                        var encoding = CanonicalEncoding(subset);
                        counts.TryGetValue(encoding, out var oldCount);
                        counts[encoding] = oldCount + (int)count;
                        if (oldCount == 0) reps[encoding] = (sbyte[][])subset.Clone();

                        if (indices[0] == Rows.Count - 1) break;
                        NextOrderedMultiset(indices, Rows.Count);
                    }

                    naiveCounts.Add(uncollapsed);
                    distinctCounts.Add(reps.Count);
                    representatives.Add(counts.Select(kvp => (kvp.Value, reps[kvp.Key])).ToList());
                }

                NaiveCounts = naiveCounts;
                DistinctCounts = distinctCounts;
                Representatives = representatives;
            }

            private static List<sbyte[]> BuildRows(int n, int ones, int negs)
            {
                var rows = new List<sbyte[]>();
                void Append(sbyte[] row)
                {
                    // Ensure that we don't have both a row and its negative by testing canonicity.
                    if (ones > negs || Array.IndexOf(row, (sbyte)1) < Array.IndexOf(row, (sbyte)-1)) rows.Add((sbyte[])row.Clone());
                }

                var perm = new sbyte[n];
                for (int i = 0; i < negs; i++) perm[i] = -1;
                for (int i = n - ones; i < n; i++) perm[i] = 1;

                while (true)
                {
                    Append(perm);

                    // Next permutation
                    int k;
                    for (k = n - 2; k >= 0; k--)
                    {
                        if (perm[k] < perm[k + 1]) break;
                    }
                    if (k == -1) return rows;

                    int l;
                    for (l = n - 1; l > k; l--)
                    {
                        if (perm[k] < perm[l]) break;
                    }

                    var tmp = perm[k]; perm[k] = perm[l]; perm[l] = tmp;

                    for (int i = k + 1, j = n - 1; i < j; i++, j--)
                    {
                        tmp = perm[i]; perm[i] = perm[j]; perm[j] = tmp;
                    }
                }
            }

            private static BigInteger CanonicalEncoding(sbyte[][] rows)
            {
                return rows.Permutations().Select(_CanonicalEncodingInner).Min();
            }

            private static BigInteger _CanonicalEncodingInner(sbyte[][] rows)
            {
                BigInteger @base = BigInteger.Pow(3, rows.Length);
                BigInteger[] digits = new BigInteger[rows[0].Length];
                for (int i = 0; i < digits.Length; i++)
                {
                    for (int j = 0; j < rows.Length; j++)
                    {
                        digits[i] = digits[i] * 3 + rows[j][i] + 1;
                    }
                }
                return digits.OrderBy(x => x).Aggregate(BigInteger.Zero, (accum, digit) => accum * @base + digit);
            }
        }

        private static readonly int[] _DeBruijn32 = new int[]
        {
            0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8,
            31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9
        };

        // Returns 32 for input 0.
        public static int CountTrailingZeros(this uint i)
        {
            if (i == 0) return 32;

            int j = (int)i;
            i = (uint)(j & -j);
            var idx = (i * 0x077CB531U) >> 27;
            return _DeBruijn32[idx];
        }

        public static IEnumerable<T[]> Permutations<T>(this IEnumerable<T> elts)
        {
            T[] arr = elts.ToArray();
            if (arr.Length == 0)
            {
                yield return arr;
                yield break;
            }

            int[] indices = Enumerable.Range(0, arr.Length).ToArray();
            while (true)
            {
                yield return arr.ToArray();
                // Next permutation
                int i = indices.Length - 1;
                while (i > 0 && indices[i - 1] > indices[i]) i--;

                if (i == 0) yield break;

                int j = indices.Length - 1;
                while (indices[j] < indices[i - 1]) j--;

                _DoubleSwap(arr, indices, i - 1, j);
                for (j = indices.Length - 1; i < j; i++, j--) _DoubleSwap(arr, indices, i, j);
            }
        }

        private static void _DoubleSwap<T>(T[] elts, int[] indices, int i, int j)
        {
            var t1 = elts[i];
            elts[i] = elts[j];
            elts[j] = t1;

            var t2 = indices[i];
            indices[i] = indices[j];
            indices[j] = t2;
        }

        private static BigInteger Factorial(int n) => n < 2 ? BigInteger.One : n * Factorial(n - 1);

        private static BigInteger Binom(int n, int k)
        {
            BigInteger result = 1;
            for (int i = 0; i < k; i++) result = result * (n - i) / (i + 1);
            return result;
        }
    }
}

Online demo
This employs the following tricks:

If there's a zero row, the permanent is zero. These can be counted very fast, and then we can assume there's no zero row.
Multiplying a row by \$-1\$ inverts the sign of the permanent. So we can assume a canonical sign for each row and then multiply by \$2^n\$.
Permuting the rows leaves the permanent unchanged, so we can assume a canonical permutation of the rows and then multiply by a suitable factor.
The permanent can be computed in better-than-brute-force time using Ryser's formula with Gray code enumeration of subsets of columns.
Laplace expansion of the last row turns a factor of \$\frac{3^n - 1}{2}\$ into an instance of the subset sum problem, which can be tackled using dynamic programming with meet-in-the-middle.
We can also permute columns, which is what Cluster is for. I think I've pushed the Cluster approach about as far as it can go, but I would like to find an efficient way of replacing it, because if we take the first row [1, 1, 0, -1, -1, -1] (with weight 60) then it would be nice to exploit the symmetry of the first two columns and the last three columns in the next cluster. The problem is that doing this naïvely adds a factor of \$n!\$, which is more expensive than the saving.
The permutations are independent, so after a bit of pre-processing we have an embarrassingly parallel problem.

I compiled with VS2017 Community targetting .Net Framework 4.7.1 and compiled for "Release" with default settings. You may prefer to use .Net Core.
\$n=6\$ is still technically out of range - I calculated it in just under 3 hours on a 4-core AMD with numThreads = 8.

Just as a permanent is zero if there is a zero row, it is zero if there are two rows with zeros at the same n-1 positions, or three rows with zeros at the same n-2 positions, and so on. These cases are not so easy to count, but we can count the matrices with non-zero permanent instead.

@ChristianSievers, interesting idea, although I don't immediately see how to apply it. The idea that I'm currently pursuing is finding a way to effectively permute columns as well as rows.

Wolfram Language(Mathematica), \$\textrm{A330134}(3)=7555\$ in under 3 seconds
Naive brute force. So slow.
Try it online!
ClearAll[zeroPermanentMatricesCount]
zeroPermanentMatricesCount[n_] :=
  Module[{allMatrices},
    (* Generate all n x n matrices with entries -1, 0, 1 *)
    allMatrices = Tuples[Tuples[{-1, 0, 1}, n], n];

    (* Count matrices with zero permanent *)
    Count[allMatrices, m_ /; Permanent[m] == 0]
  ]
Do[Print["n = ", n, ", ", zeroPermanentMatricesCount[n]], {n, 1, 3}]


text:
Roll for Initiative!
Roll for Initiative!
Introduction
In tabletop games like Dungeons and Dragons, when you begin a battle, all involved parties roll for initiative. In DnD 5e, this is 1d20 + DEX + Other bonuses, where DEX is the bonus given by your Dexterity stat. The characters that roll higher numbers go first. We'll use a similar, deterministic system in this challenge.
The Challenge
Write a program or function that, when given a list of characters, will output a list of characters in order of initiative.
A character is defined as this:
character = {
    name: "name" // a string
    statblock: [SPD, DEX, WHT] // a list of numbers
                               // DEX = dexterity, SPD = speed, WHT = weight
}

The formula for initiative is the following:
$$\text{Initiative} = \left\lfloor{ \frac{\text{SPD}^2}{\sqrt{\lvert\text{DEX}\rvert}} }\right\rfloor - \text{WHT}$$
Input
A list of characters, unsorted. This can be a JSON object, a list of lists, a list of dictionaries, a series of strings etc.
It is guaranteed that all names will be unique.
Output
A list of characters, or character names, sorted by initiative order from highest to lowest, based on the above formula.
Rules

IO can be used with any convenient method.
Standard loopholes are forbidden.
This is code-golf, so smallest program wins!

Sample IO
Input --> Output
[[Name, SPD, DEX, WHT], ...]
    --> [[Name, SPD, DEX, WHT], ...] (or [Name, Name, ...])
---------
[[Alice,1,2,3],[Bob,10,5,0],[Charlie,3,2,1]]
    --> [Bob, Charlie, Alice]
// Alice = -3, Bob = 44, Charlie = 5

[[Z,1,1,1],[B,1,1,1],[XY,5,1,1]]
    --> [XY, Z, B]
// Retain the order of characters from the input if they have the same initiative.
// Z = 0, B = 0, XY = 24

[[Neg,-3,-3,-1],[SomeNeg,5,-2,-4],[NoNeg,4,6,8]]
    --> [SomeNeg, Neg, NoNeg]
// Negative values are valid.
// Neg = 6, SomeNeg = 21, NoNeg = -2

[[Flo,1.5,2.5,3.5],[MoreFlo,2,2.5,3.5]]
    --> [[MoreFlo,2,2.5,3.5], [Flo,1.5,2.5,3.5]]
// Floats are also valid.
// Flo = -2.5, MoreFlo = -1.5

[[Lonely,1,2,3]]
    --> [[Lonely,1,2,3]]
// Input with 1 item.

[]
    --> []
// Empty input leads to empty output.

[[Foo,2,1,4], [Baz,5,5,4], [Bar,5,1,4]]
    --> [Bar, Baz, Foo]
// Foo = 0, Bar = 21, Baz = 7

[['Adam', 2, 4, 1], ['Eve', 2, 3, 1]]
--> [Adam, Eve]
// Adam = 2, Eve = 2
// If you do not floor, you'll end up with [Eve, Adam] (Adam = 2, Eve ~= 2.3)

Sandbox link

Curious - where did this initiative formula come from?

Something I came up with on the fly while writing this, it's not really based on anything.

So it is. Thanks for catching it!

Can we sort the input in place (i.e whatever is passed into the function will be sorted when the function finishes)?

So small dexterity is good here?

For some reason ignoring DEX parameter does not affect the final result.

I just added a new test case. You should not ignore DEX. If you look at Baz in the last test case, you'll get 7 if you use DEX, and 21 if you don't.

@bigyihsuan But you get Bar, Baz, Foo even if ignore it: check here, or do I miss something here?

@game0ver Your formula fails for [['A', 1, 1, 5], ['B', 2, 50, 5]]. With the correct formula, you get [A, B], your formula returns [B, A]

@BlackOwlKai thanks, your are right, didn't check that... Also I ended up doing some math and it turned out that the formula I used was completely arbitrary. There is an alternative way to write the formula but it doesn't save bytes unfortunately...

Suggested testcase: [['Adam', 2, 4, 1], ['Eve', 2, 3, 1]]. Both have an initiative of 2, so input order should be retained and the output should be ['Adam', 'Eve']. But if a solution does not round down the first part, Adam stays at 2, but Eve gets to ~2.3, therefore the solution outputs ['Eve', 'Adam']

JavaScript (ES7),  73  60 bytes
-13 bytes by just sorting the input without isolating the names, as suggested by @asgallant
Expects a list of [SPD, DEX, WHT, Name].
a=>a.sort((a,b)=>(g=([s,d,w])=>~(s*s/(d*d)**.25)+w)(a)-g(b))

Try it online!
How?
To save a few bytes:

we compute \$-\text{Initiative}-1\$ instead and sort the list the other way around
we compute \$\sqrt{|\text{DEX}|}\$ as \$(\text{DEX}^2)^{\frac{1}{4}}\$

Outputting a list of characters (not just their names) is valid according to the challenge, so you can save a few bytes by dropping the .map call.

05AB1E, 10 8 bytes
Σ`nsÄt÷-

Input in the format and order ["name", WHT, DEX, SPD], and output are the same inner lists sorted from highest to lowest initiative.
-2 bytes thanks to @Grimmy.
Try it online or verify all test cases.
Explanation:
The sort builtin Σ sorts from lowest to highest by default, so I'm using the following modified formula to not only save bytes, but sort from highest to lowest at the same time with this builtin (thanks @Grimmy):
$$Initiative = WHT - \left\lfloor\frac{SPD^2}{\sqrt{\lvert DEX\rvert}}\right\rfloor$$
Σ         # Sort the (implicit) input-list by:
 `        #  Pop and push all values of the inner list separated to the stack
  n       #  Square the SPD at the top of the stack
   s      #  Swap so the DEX is now at the top of the stack
    Ä     #  Take its absolute value
     t    #  Then its square-root
      ÷   #  And integer-divide the squared SPD by this
       -  #  Subtract this from the WHT (to sort in descending instead of ascending order)
          #  (after which the result is output implicitly)

s-( can be just -.

Japt, 20 18 bytes
ñÈÌ-(Xg1 ²/Xg2a ¬f

Try it
Simply calculates the formula for each character and sorts. Negates the formula to sort in the correct order.
Takes input as a list of lists [Name, SPD, DEX, WHT]. (I could save a byte taking input as [SPD, Name, DEX, WHT] but that feels cheaty.)

Python 3, 63 59 bytes
lambda c:sorted(c,key=lambda x:x[3]-x[1]**2//abs(x[2]**.5))

Try it online!
Pretty straigtforward. Sorts by - initiative to get the right output order.
Python 2, 60 56 bytes
lambda c:sorted(c,key=lambda(a,b,c,d):d-b*b//abs(c)**.5)

Try it online!
thanks to @Chas Brown

In python2, you can shave off 2 bytes via unpacking. (This is one of the few things I don't like in Python 3).

I think you are supposed to return a list of just the names,not the entire records.

@GalenIvanov The Output section of the challenge specifies you can return the whole Character, or the characters name. This answer has chosen the first option.

@Potato44 Yes, the Sample IO shows different possibilities.

@ChasBrown Thank you, I tried this because I thought I had seen it somewhere, but it resulted in an error and I didn't think that it could be 2-only. Got it to 60 with b**2 -> b*b.

@BlackOwlKai You don't need int(...), you can save 5 bytes... The same goes also for python2

K (oK), 28 bytes
f:{(!x)@>{-z-(x*x)%%y|-y}.'.x}

Try it online!
-10 thanks to ngn :-)

(y;-1*y)0>y -> y|-y, %: -> %, (-z)+ -> -z-

@ngn aghhh, the y sign switch was obviously doing too much work in my implementation. your solution seems to obvious in hindsight. i won't forget that!

Ruby, 44 bytes
->d{d.sort_by{|_,x,y,z|-x*x/y.abs()**1/2+z}}

Try it online!
-3 bytes thanks to @79037662

You can get rid of the reverse by sorting "backwards": https://tio.run/##VczNC4IwHMbx@/6K4UnXb@qmRheD6t6l4xix@VIDS/HlIOrfbiYSBM/pw5en7nQ/5/FMj@mQuk1Zt3fdD6MZqRFMEsI9I7h0lW5shxDm8Z0RgZymWeEYC4SxsE6FSTILGHAIJKx0LvUCPkTgb3J5qrowSxYsGZNIIv07uGYPC2iwjm39rXxlq0dAOdBw42u5Ygh7OHxPqtq8W5y7iSoKWzmo6trmH7UzfwA

@79037662 thanks! Good point! I also got rid of some extra bytes by replacing indexing with variables!

See the last test case, which fails if you ignore DEX.

No it doesn't, it gives Bar, Baz, Foo which is the test case...

Add++, 24 bytes
D,g,@#~,p2^$|.5^/i_
L,§g

Try it online!
Inputs and outputs a list of lists representing the characters, in the same order given in the challenge ([Name SPD DEX WHT])
How it works
On the second line, we define an anonymous function that takes in an argument and sorts it on the result of g, defined on the first list.
g starts by reversing its argument and splatting the elements to its stack: #~. It then removes the name and squares the first stat, Speed, (p2^) before square rooting the absolute value of the second stat, Dexterity, ($|.5^). Finally, we integer divide the two (/i) and subtract from Weight (_). This is the opposite to the formula (where Weight is subtracted from the division), as it will cause the sort function in the second line to sort in the order we want. As Add++ uses a stable sort, order is retained if two characters have the same initiative.

APL+WIN, 41 bytes
Prompts for a n x 4 matrix where n is the number of rows 1 per character and the columns are Name, SPD, DEX, WHT. The output is a 1 column matrix of character names
m←⎕⋄⊃m[⍒⌊((m[;2]*2)÷(|m[;3])*.5)-m[;4];1]

Try it online!Courtesy of Dyalog Classic

Red, 108 107 bytes
-1 byte thanks to raznagul
func[b][extract next sort/skip collect[foreach e
b[keep 0 -(e/2 ** 2 /(e/3 ** 2 ** .25)- e/4)keep e/1]]2 2]

Try it online!

I don't know red, but shouldn't you be able to replace e/2 * e/2 with e/2 ** 2?

@raznagul Yes,of course. You can see that I did it for e/3 but apparently I forgot to apply it to e/2. Thanks!

Google Sheets, 49 37 bytes
-12 bytes thanks to Black Owl Kai
=Sort(A:A,D:D-B:B^2/Sqrt(Abs(C:C)),1)

Input is in the range A1:D with the four columns being player, speed, dexterity, and weight.
Explanation:
D:D-B:B^2/Sqrt(Abs(C:C)) do the math bit. The values are actually Initiaive * -1 so we can sort them ascending and get the order we want. Otherwise, all the blank rows would be sorted before the player names.
=Sort(A:A,~,1) sorts the player list based on the math bit output.
Test cases:

=Sort(A:A;D:D-B:B^2/Sqrt(Abs(C:C));1) gets rid of the if

Jelly, 20 bytes
ḊA2¦U*Ø½j.¤UŒH:/€IµÞ

Try it online!
Well, this is horribly hacked together but it works. The Footer in the TIO link simply runs the code over each list then returns just the names. The code itself returns the names and all stats in the same format as the input, namely [Name, SPD, DEX, WHT]
Will post an explanation after any more golfing

Python 3.8, 81 bytes
This solution uses the formula as described in the question.
As input a dictionary is expected in the following form: {"name":[SPD, DEX, WHT],...}
lambda d:sorted(n:={_:(i:=-a*a/abs(b)**.5+c)for _,(a,b,c)in d.items()},key=n.get)

Try it online!

See the last test case, which fails if you ignore DEX.

No it doesn't, it gives Bar, Baz, Foo which is the test case...

Charcoal, 36 bytes
ＦＡ⊞υ⟦⁻⌊∕Ｘ⊟ι²₂↔⊟ι⊟ι⊟ι⟧Ｗ⌈υ«≔Φυ¬⁼ικυ⟦⊟ι

Try it online! Takes quadruples in the order name, weight, speed, dexterity. Explanation:
ＦＡ

Loop over the quadruples.
⊞υ⟦⁻⌊∕Ｘ⊟ι²₂↔⊟ι⊟ι⊟ι⟧

Square the dexterity, floor divide by the square rooted absolute speed and subtract the weight, then create a tuple with the result and the name and push that to the predefined list.
Ｗ⌈υ«

While there are still tuples in the list, take the maximum.
≔Φυ¬⁼ικυ

Remove the tuple from the list.
⟦⊟ι

Output the name on its own line.

C++ (clang), 324 \$\cdots\$ 232 204 bytes
Saved 18 46 bytes thanks to cielingcat!!!
#import<regex>
#import<cmath>
#define i(s)int(*s**s/sqrt(abs(s[1])))-s[2]
struct C{std::string n;float s[3];};void f(std::vector<C>&v){std::stable_sort(&v[0],&*end(v),[](C a,C b){return i(a.s)>i(b.s);});}

Try it online!
Ungolfed:
#include <string>
#include <algorithm>
#include <cmath>
#include <vector>
using namespace std;

struct C {
    string n;
    float s[3];
};
float i(C c) {
    return floor(c.s[0]*c.s[0] / sqrt(abs(c.s[1]))) - c.s[2];
}
void f(vector<C>&v) {
    stable_sort(begin(v), end(v),
        [](C a,C b) {
            return i(a) > i(b);
        });
}

Jelly, 15 14 bytes
ṛ²:A½$}ɗ_@3ƭ/Þ

Try it online!
A monadic link that takes a list of characters as its argument. Each list is structured as [Name, SPD, DEX, WHT]. Returns a list of characters in the desired order.

APL(NARS), chars 40, bytes 80
{⍵≡⍬:⍬⋄↑¨⍵[⍒{(4⊃⍵)-⍨⌊(2*⍨2⊃⍵)÷√∣3⊃⍵}¨⍵]}

test:
  f←{⍵≡⍬:⍬⋄↑¨⍵[⍒{(4⊃⍵)-⍨⌊(2*⍨2⊃⍵)÷√∣3⊃⍵}¨⍵]}
  ⎕fmt ,⊂(('Lonely') 1 2 3)
┌1────────────────┐
│┌4──────────────┐│
││┌6──────┐      ││
│││ Lonely│ 1 2 3││
││└───────┘ ~ ~ ~2│
│└∊──────────────┘3
└∊────────────────┘
  ⎕fmt f,⊂(('Lonely') 1 2 3)
┌1────────┐
│┌6──────┐│
││ Lonely││
│└───────┘2
└∊────────┘
  ⎕fmt (('Alice') 1 2 3)(('Bob') 10 5 0)(('Charlie') 3 2 1)
┌3──────────────────────────────────────────────────┐
│┌4─────────────┐ ┌4────────────┐ ┌4───────────────┐│
││┌5─────┐      │ │┌3───┐       │ │┌7───────┐      ││
│││ Alice│ 1 2 3│ ││ Bob│ 10 5 0│ ││ Charlie│ 3 2 1││
││└──────┘ ~ ~ ~2 │└────┘ ~~ ~ ~2 │└────────┘ ~ ~ ~2│
│└∊─────────────┘ └∊────────────┘ └∊───────────────┘3
└∊──────────────────────────────────────────────────┘
  ⎕fmt f (('Alice') 1 2 3)(('Bob') 10 5 0)(('Charlie') 3 2 1)
┌3─────────────────────────┐
│┌3───┐ ┌7───────┐ ┌5─────┐│
││ Bob│ │ Charlie│ │ Alice││
│└────┘ └────────┘ └──────┘2
└∊─────────────────────────┘
  ⎕fmt f ('Z' 1 1 1)('B' 1 1 1)(('XY') 5 1 1)
┌3────────┐
│┌2──┐    │
││ XY│ Z B│
│└───┘ ¯ ¯2
└∊────────┘
  ⎕fmt f (('Neg') ¯3 ¯3 ¯1)(('SomeNeg') 5 ¯2 ¯4)(('NoNeg') 4 6 8)
┌3─────────────────────────┐
│┌7───────┐ ┌3───┐ ┌5─────┐│
││ SomeNeg│ │ Neg│ │ NoNeg││
│└────────┘ └────┘ └──────┘2
└∊─────────────────────────┘
  ⎕fmt f ⍬
┌0─┐
│ 0│
└~─┘
  ⎕fmt f (('Flo') 1.5 2.5 3.5)(('MoreFlo') 2 2.5 3.5)
┌2────────────────┐
│┌7───────┐ ┌3───┐│
││ MoreFlo│ │ Flo││
│└────────┘ └────┘2
└∊────────────────┘
  ⎕fmt f (('Foo') 2 1 4)(('Baz') 5 5 4)(('Bar') 5 1 4)
┌3───────────────────┐
│┌3───┐ ┌3───┐ ┌3───┐│
││ Bar│ │ Baz│ │ Foo││
│└────┘ └────┘ └────┘2
└∊───────────────────┘
  ⎕fmt f (('Adam') 2 4 1)(('Eve') 2 3 1)
┌2─────────────┐
│┌4────┐ ┌3───┐│
││ Adam│ │ Eve││
│└─────┘ └────┘2
└∊─────────────┘

Perl 5 -n, 72 bytes
sub i{int($_[1]**2/sqrt abs$_[2])-pop}say$$_[0]for sort{i(@$b)-i@$a}eval

Try it online!

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Who's next to me in the queue?
Problem 4 in the 2019 BMO, Round 1 describes the following setup:

There are \$2019\$ penguins waddling towards their favourite restaurant. As
  the penguins arrive, they are handed tickets numbered in ascending order
  from \$1\$ to \$2019\$, and told to join the queue. The first penguin starts the queue.
  For each \$n > 1\$ the penguin holding ticket number \$n\$ finds the greatest \$m < n\$
  which divides \$n\$ and enters the queue directly behind the penguin holding
  ticket number \$m\$. This continues until all \$2019\$ penguins are in the queue.

The second part of the question asked candidates to determine the penguins standing directly in front of, and directly behind, penguin \$33\$. This could be done by examining the patterns in the queue, considering prime factors: see the online video solutions for more information.

The Challenge
Your task is to design a program or function which, given a positive integer \$k\$ representing the penguin with ticket number \$k\$, outputs the ticket numbers of the penguins directly before and after this penguin.
For example, penguin \$33\$, stands directly behind \$1760\$ and directly in front of \$99\$, so the program should output, in some reasonable format, \$[1760, 99]\$.

Rules

The input will be an integer in the range \$1 < k \le 2019\$.
Your program should output two integers, in any reasonable format, representing the ticket numbers of the penguins before and after.
These can be output in any order, (front first or behind first) but this order must be consistent.
The penguin will not be at the front or back of the queue: so you don't have to handle the edge cases of \$k = 1\$ or \$k = 1024\$.
As penguins find it difficult to read human glyphs, your program should be as short as possible. This is a code-golf - so the shortest program (in bytes) wins!

Test Cases
These outputs are given in the format [front, behind].

33   -> [1760, 99]
512  -> [256, 1024]
7    -> [1408, 49]
56   -> [28, 112]
1387 -> [1679, 1241]
2019 -> [673, 1346]
2017 -> [1, 2011]
2    -> [1536, 4]

@game0ver It doesn't, it joined the queue behind n = 11, but then subsequent penguins whose numbers are multiples of 11 but not 33 pushed in front of it.

Ah, thanks, my mistake, I though that $m$ was supposed to be 1760.

Perl 6, 79 76 71 bytes
{($!=sort {[R,] -$_,{first $_%%*,$_^..1}...1},^2020)[+(@$!...$_)X-2,0]}

Try it online!
Sorts numbers 1 to 2019 lexicographically by the reversed, negated sequence of largest divisors, then finds the neighbors. Example mappings:
1760 => (-1 -11 -55 -110 -220 -440 -880 -1760)
  33 => (-1 -11 -33)
  99 => (-1 -11 -33 -99)

In other words, the negated cumulative product of prime factors in descending order.

JavaScript (ES6),  105 ... 101  100 bytes
A straightforward implementation that builds the full queue and looks for the penguin with the input ticket within it.
Returns [front, behind].
k=>[eval("for(n=q=[];m=++n%2020;q.splice(q.indexOf(m),0,n))while(n%--m)q")[i=q.indexOf(k)+1],q[i-2]]

Try it online!
Commented
The code in eval() builds and returns the queue in reverse order:
for(                      // outer loop:
  n = q = [];             //   q[] = queue, n = counter (initially zero'ish)
  m = ++n % 2020;         //   increment n; set m to n mod 2020, so that we
                          //   stop when n = 2020
  q.splice(               //   after each iteration:
    q.indexOf(m), 0, n    //     insert n before m in q[] (should be *after* m to build
  )                       //     the queue from first to last, but it's shorter this way)
) while(n % --m)          //   inner loop: decrement m until it divides n
    q                     //     dummy loop statement so that q[] is returned by eval()

which gives \$q[\:]=[1024,512,...,2017,1]\$.
Wrapper code:
k => [                    // k = input
  eval("...")             // build q[]
  [i = q.indexOf(k) + 1], // return the element after k ('in front of')
  q[i - 2]                // and the element before k ('behind')
]                         //

"straighforward" - Arnauld, 2019 :D I find this evil ^^

@Kaddath Did you mean I find this eval? :p

Yeah but not only, the way loops are used to build and return q are really clever!

Jelly, 15 bytes
⁽¥ØÆfṚNƊÞṣị"@Ø.

A monadic Link accepting an integer between 1 and 2019 * which yields a list of two integers, [in-front, behind].
* If 1 or 1024 is input the missing ticket number will be shown as 0.
Try it online!
How?
⁽¥ØÆfṚNƊÞṣị"@Ø. - Link: n
⁽¥Ø             - literal 2019
        Þ       - sort (range [1..2019]) by:
       Ɗ        -   last three links as a monad:
   Æf           -     prime factorisation (e.g. 1100 -> [2,2,5,5,11])
     Ṛ          -     reversed
      N         -     negated
         ṣ      - split at (n)
             Ø. - literal [0,1]
            @   - with swapped arguments:
           "    -   zip with:
          ị     -     index into (1-based & modular)
                      - i.e. [last of left list, first of right list]

Python 2, 137 133 121 bytes
q=[1]
for n in range(2,2020):j=q.index(max(m for m in q if n%m<1))+1;q=q[:j]+[n]+q[j:]
print q[q.index(input())-1::2][:2]

Try it online!
4 bytes thx to Arnauld; 11 bytes thx to FlipTack.
Another naive implementation.

05AB1E, 16 15 bytes
Ž7ëLΣÒR(}ûI¡€θ¨

Port of @JonathanAllan's Jelly answer, so make sure to upvote him!
-1 byte thanks to @Grimmy.
Outputs in the order [front, behind].
Try it online or verify all test cases.
Explanation:
Ž7ë              # Push compressed integer 2019
   L             # Pop and push a list in the range [1,2019]
    Σ            # Sort this list by:
     Ò           #  Get the prime factors (with duplicates)
      R          #  Reverse it
       (         #  Negate each inner value
        }û       # After the sort: palindromize this list ([a,b,c] → [a,b,c,b,a])
          I¡     # Split this pallindromized list by the input-integer
            €θ   # For each inner list: only leave the last value
              ¨  # And remove the last value, so the pair remains
                 # (after which this pair is output implicitly as result)

See this 05AB1E tip of mine (section How to compress large integers?) to understand why Ž7ë is 2019.

-1 using ûI¡€θ¨ (ûDIQÀÏ also works, for the same byte-count).

There's also ûʒQYsV, for a slightly weirder 15.

@Grimmy Thanks! Nice use of the palindromize builtin! And that third one is indeed weird. Took me a few looks to understand what was going on. But again pretty smart. :)

C++ (clang), 288 \$\cdots\$ 210 203 bytes
Saved a whopping 22 29 bytes thanks to ceilingcat!!!!
Saved a byte thanks to @AZTECCO!!!!
#import<bits/stdc++.h>
#define l find(begin(q),end(q)
using v=std::vector<int>;v f(int p){v q{1};for(int n=1,d,i;d=i=++n<2020;q.insert(l,d),n))for(;++i<n&d<2;)d=n%i?d:n/i;auto t=l,p);return{*++t,t[-2]};}

Try it online!
Ungolfed:
#include <vector>
#include <algorithm>

auto f(int p)
{
    std::vector<int> q{1};
    for (int n = 2; n <= 2019; ++n)
    {
        int d = 1;
        for (int i = 2; i < n; ++i)
        {
            if (n % i == 0)
            {
                d = n / i;
                break;
            }
        }
        auto it = find(begin(q), end(q), d);
        q.insert(it, n);
    }
    auto it = find(begin(q), end(q), p);
    std::vector<int> r = {*(it + 1), *(it - 1)};
    return r;
}

Perl 5, 88 bytes
sub f{$_=1;for$n(2..2019){map{$n%$_ or$m=$_}1..$n/2;s/\b$m\b/$&,$n/}/(\d+),$_[0],(\d+)/}

Try it online!
Same ungolfed:
sub f {
  $_=1;                             # $_ is the comma separated queue "array string"
                                    # ...with 1 as the first penguin
  for $n (2..2019){                 # process the rest of the penguins 2-2019
    map { $n%$_ or $m=$_ } 1..$n/2; # find max $m divisible by current $n
    s/\b$m\b/$&,$n/                 # search-replace to place current $n behind $m
  }
  /(\d+),$_[0],(\d+)/               # find and return the two numbers:
                                    # the one before and the one after the input
                                    # parameter $_[0] in the $_ array string
}

R, 90 89 bytes
for(i in 2:2019)T=append(T,i,match(max(which(!i%%1:(i-1))),T));T[match(scan(),T)+c(-1,1)]

Try it online!
Builds up the list into T and extracts the relevant values.
As a bonus, it works for 1 and almost works for 1024, as R's 1-based indexing returns nothing for index 0 and NA for index 2020.

Jelly, 23 21 bytes
ḍƇṀ=⁸k⁸j
⁽¥ØRç/ṣ⁸ṪḢƭ€

Try it online!
Naïve answer that simply generates the list of penguins and finds the relevant place. A pair of links which, when called as a monad, takes an integer as its argument and returns a pair of integers indicating the penguins before and after in the list of penguins.

Ruby, 104 101 95 bytes
->n{*w=1;(2..2019).map{|z|w.insert(1+w.index(w.select{|x|z%x<1}.max),z)};w[w.index(n)-1,3]-[n]}

Try it online!

Retina 0.8.2, 109 bytes
^
2019$*_
_
$`_¶
(_+)\1+¶
$1 $&
+m`^(_+)¶((_+¶)*)\1 (_+¶)
$1¶$4$2
_+
$.&
^(.+¶)*(.+)¶(.+)(¶.+)(¶.+)*¶\3$
$2$4

Try it online! Note: Link only goes to 201 as 2019 is too slow for TIO. Explanation:
^
2019$*_

Prepend 2019 in unary.
_
$`_¶

Count from 1 to 2019.
(_+)\1+¶
$1 $&

Precede each number with its highest proper factor (except 1, which remains unchanged).
+m`^(_+)¶((_+¶)*)\1 (_+¶)
$1¶$4$2

Insert in turn each number after its highest proper factor.
_+
$.&

Convert to decimal.
^(.+¶)*(.+)¶(.+)(¶.+)(¶.+)*¶\3$
$2$4

Extract the predecessor and successor of the original input.
118 byte version that's fast enough to handle 2019 on TIO:
^
2019$*_
_
$`_¶
(_+)\1+\b
$.1;$.&
_
1
m{`^(\d+)¶\1;
$1¶
}`^(\d+¶)(\d+;\d+¶)
$2$1
^(.+¶)*(.+)¶(.+)(¶.+)(¶.+)*¶\3$
$2$4

Try it online! Explanation:
^
2019$*_
_
$`_¶

Prepend all of the unary numbers from 1 to 2019.
(_+)\1+\b
$.1;$.&

Convert to decimal and prepend the largest factor.
_
1

Except 1, which just becomes the first entry in the list.
m{`^(\d+)¶\1;
$1¶
}`^(\d+¶)(\d+;\d+¶)
$2$1

Bubble up each line until it finds its insertion point.
^(.+¶)*(.+)¶(.+)(¶.+)(¶.+)*¶\3$
$2$4

Extract the predecessor and successor of the original input.

Japt, 36 bytes
1k
#É9ò2@iUaXâ n- gJÑ)ÄX
g[J1]c+UaNg

Try it
1k               assign [1] to U

#É9ò2            pass range 2...2019 to function:
     @i ...  X   insert each value in U at index:
 Ua               find max divisor of X in U:
   Xâ              divisors
      n-           sorted
         gJÑ)Ä     2nd to last

g                get elements at
 [J1]             [-1,1]
     c+UaNg       add index of 1st input in U

Scala, 163 bytes
n=>{val a=(2 to 2019).map(x=>(x,((1 to x/2).filter(x%_==0)).max)).foldLeft(Seq(1))((l,e)=>l.patch(l.indexOf(e._2)+1,Seq(e._1),0));val i=a indexOf n;(a(i-1),a(i+1))

Try it online!

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Generate an ASCII art representation of a ribbon plot
Sandbox post
In this code-golf challenge, you are provided with a series of y values at even x-values as well as a console width and height. The desired output is an ASCII-art representation of a ribbon plot of the y-values scaled to the desired size. If all values are positive or negative, the x-axis should be displayed at the bottom or top respectively. Otherwise, the x-axis should be displayed at the appropriate position with the graph extending to the top and bottom.
The data should be scaled horizontally by splitting the data into sublists
 of equal (+/- 1) size and taking the arithmetic mean of each sublist. There will never be fewer values than columns, but can be more.
Inspired by this question which was posed as a tips in C question. Most similar question is here but there are some important differences in terms of input (floats versus equation) and output style.
Default loopholes apply.
Input/output

Default I/O rules apply.
Default loopholes forbidden.
The y-values will be a list, array or equivalent of floating point numbers at least as long as the width of the console. These may be positive or negative. Choice of float type is at the discretion of the answerer. Y values should be rounded to the nearest integer after scaling to the relevant size; the default rounding method for 0.5s supported by the language is fine (rounding to even, rounding away from zero, etc.) (See Wikipedia for a discussion of different rounding methods if this is an unfamiliar topic.) However, truncation/floor/ceiling are not permitted approaches.
The range of y values will always be non-zero (i.e. there will be more than 1 y value, and there will be at least 2 different values).
The console width and height will be two integers, with the height always at least 2.
The output can be a return value or be output to STDOUT. If opting for a return value, this can be a single newline-separated string, a list of strings, a list of lists of characters or a matrix of characters.
Any distinct character can be chosen for the filled area, but the blank area should be a space character.
A separate non-space character should be used to indicate the x-axis.
Trailing space on each line is optional.

Worked example of how the data is scaled to width
Floats: 1.0,2.0,3.0,4.0,5.0
Width: 3

Split into three pieces: [1,2],[3,4],[5]
Take arithmetic mean of each: [1.5,3.5,5]

Examples of full input/output
Input
Floats: [1.4330127018922192, 1.4740546219716173, 1.46726631732171, 1.4118917879852302, 1.3095749615197134, 1.1642846441770978, 0.98209964435211, 0.770867389064186, 0.5397574401179837, 0.29873801741462563, 0.05800851721072442, -0.17257624372014382, -0.38405389828888165, -0.5688402657187754, -0.7211228325469043, -0.8371379250265838, -0.915315134687127, -0.9562812209887335, -0.9627248583907679, -0.939132622452629, -0.8914149315594007, -0.8264477348143321, -0.7515611018951998, -0.6740091786833334, -0.6004570095761375, -0.5365184287338843, -0.48637567624364625, -0.45250583664959154, -0.4355319935803172, -0.43420862005149424, -0.445541730745404, -0.46503530191616943, -0.4870470099638535, -0.5052290073886498, -0.5130237273726901, -0.5041809629881258, -0.47326095136808916, -0.41608900827923856, -0.330130352731954, -0.21475893339455088, -0.07140098112303239, 0.09645778693148394, 0.28340343297327375, 0.4823565362517651, 0.6849507176587676, 0.8820035800170399, 1.0640505126834876, 1.2219075413197402, 1.3472273300168627, 1.433012701892219]
Width: 24
Height: 24

Output:
#
#                      #
##                     #
##                    ##
##                    ##
##                    ##
###                   ##
###                  ###
###                  ###
###                  ###
####                ####
####                ####
####                ####
####                ####
------------------------
    ###############
     ##############
     ##############
     #############
     #######   ##
      #####
      #####
       ###
       ##

Input
Floats: [0.5, 3, 7.5, 14, 22.5, 33, 45.5, 60, 76.5, 95, 115.5, 138, 162.5, 189, 217.5, 248, 280.5, 315, 351.5, 390, 430.5, 473, 517.5, 564, 612.5, 663, 715.5, 770, 826.5, 885, 945.5, 1008, 1072.5, 1139, 1207.5, 1278, 1350.5, 1425, 1501.5, 1580, 1660.5, 1743, 1827.5, 1914, 2002.5, 2093, 2185.5, 2280, 2376.5, 2475]
Width: 24
Height: 24

Output
                       #
                       #
                      ##
                      ##
                     ###
                    ####
                    ####
                   #####
                  ######
                  ######
                 #######
                ########
               #########
               #########
              ##########
             ###########
            ############
           #############
          ##############
        ################
       #################
     ###################
  ######################
------------------------

Input
Floats: [-4.0,0.0,0.0,4.0,2.0]
Width: 5
Height: 5

Output
   #
   ##
-----
#
#

Input
Floats: [0.0,1.0,2.0,3.0]
Width: 4
Height: 4

Output
    #
   ##
  ###
 ----

@JonathanAllan The input isn't sorted in the first example, nor in the third example. I've provided a worked example of how the splitting of data should work, but in Jelly it would be equivalent to œsÆm€ where the left argument is the floats and the right argument the width.

Ahha the example makes sense. Thanks.

I may be missing something but judging by the example with only positive values it looks like maybe you're flooring rather than rounding (I see un-rounded, normalised values like [0.036, 0.229, 0.522, 0.848, 1.254, 1.739, 2.303, 2.946, 3.668, 4.47, 5.35, 6.31, 7.349, 8.467, 9.664, 10.94, 12.296, 13.73, 15.244, 16.837, 18.509, 20.26, 22.09, 24.0])

@JonathanAllan thanks, fixed, I think. Could you check please?

It would be great if anyone who downvotes could post a comment to explain why; I’ve answered lots of questions here but this my first attempt at asking one!

I assume we're allowed to not have rounding built into the programs, which gives slightly different results in comparison to your test cases? Take the first test case for example: in your output it has 2 bars reaching all the way to the bottom, whereas both mine and @JonathanAllan's answers have 3 bars reaching all the way to the bottom. This is - I assume - because you've rounded when scaling, where me and Jonathan didn't.

If you split 7 values into 5 pieces, what's the desired result?

@Neil Something like [[1,2],[3,4],[5],[6],[7]]

Jelly,  40  33 bytes
-7 thanks to NickKennedy!
œsÆm÷_ṂṀƊ$×⁵’¤ær0r0FṀ‘,_¥ƲṬSz0o⁶Y

A full program printing the result to STDOUT with 1s for bars and 2s for the x-axis.
Try it online!
How?
œsÆm÷_ṂṀƊ$×⁵’¤ær0r0FṀ‘,_¥ƲṬSz0o⁶Y - Main Link: values, width
œs                                - split (values) into (width) chunks
  Æm                              - arithmetic mean (vectorises)
         $                        - last two links as a monad:
        Ɗ                         -   last three links as a monad:
      Ṃ                           -     minimum
     _                            -     (left) subtract (minimum) (vectorises)
       Ṁ                          -     maximum
    ÷                             -   divide (vectorises)
             ¤                    - nilad followed by link(s) as a nilad:
           ⁵                      -   third program argument (height)
            ’                     -   decrement
          ×                       - mutiply (vectorises)
              ær0                 - round to 0 decimal places (vectorises)
                 r0               - inclusive range to zero (vectorises)
                         Ʋ        - last four links as a monad:
                   F              -   flatten
                    Ṁ             -   maximum
                     ‘            -   increment
                        ¥         -   last two links as a monad:
                       _          -     subtract (vectorises)
                      ,           -     pair
                          Ṭ       - untruth (vectorises)  (e.g. 3 -> [0,0,1])
                           S      - sum (vectorises)
                            z0    - transpose with filler zero
                              o⁶  - logical OR with a space character (i.e. replace zeros)
                                Y - join with newline characters
                                  - implicit, smashing print

Nice. I’ll allow the leading whitespace line.

Is it ok to help improve an answer to my own question? If so, you can save 7 bytes using an alternative second half: https://tio.run/##AVwAo/9qZWxsef//xZNzw4Ztw7df4bmC4bmAxookw5figbXigJnCpMOmcjByMEbhuYDigJgsX8KlxrLhuaxTejBv4oG2Wf///y01LDEsMiwzLDQsNSw2LDcsOCw5/zX/NQ

I had planned to come back to this as per the "I feel like this is probably a little golf-able (what with all the Zzs)...", that's why I had not written up a code break-down either. Looks like you've done it for me; It's a really nice way to go about it - thanks!

05AB1E, 41 bytes
äÅAWUZX-I/I<LR*X+©δ@øD®dÏs®d_Ï_¹Å8s)J˜0ð:

Inputs in the order width, [floats], height.
Outputs as a list of lines with 1 for the bars and 8 for the x-axis.
Also doesn't round at all for more accurate scaling (similar output as the Jelly answer, which differs from the test cases for some bars).
Try it online or verify all test cases. (Uses » in the footer to join the resulting list of lines by newlines to pretty-print it. The test suite will have the input-order as height, width, [floats] instead.)
Explanation:
ä                # Split the second (implicit) input-list `floats` into the first (implicit)
                 # input-integer `width` amount of equal-sized parts
 ÅA              # Take the arithmetic mean of every inner sublist
   W             # Push the minimum of this list (without popping the list itself)
    U            # Pop and store it in variable `X`
   Z             # Push the maximum of this list (again without popping the list itself)
    X-           # Subtract the minimum `X` from it
      I/         # Divide it by the third input `height`
        I<L      # Push a list in the range [1, `height`-1]
           R     # Reverse it to [`height`-1, 1]
            *    # Multiply each value by the (max-min)/height we calculated
             X+  # And add the minimum `X`
               © # Store this list of y-axis steps in variable `®` (without popping)
   δ             # Apply double-vectorized on the arithmetic mean and y-step lists:
    @            #  Check >= among the two
     ø           # Zip/transpose this matrix (swapping rows/columns)
      D          # Duplicate it
       ®         # Push the y-steps again from variable `®`
        d        # Check for each whether its non-negative (>= 0)
         Ï       # And only leave those inner sublists of the matrix (top part)
      s          # Swap to get the duplicated list again
       ®d_Ï      # Do the same, but this time for all negative values (bottom part)
           _     # Invert each (0 becomes 1; everything else becomes 0)
      ¹Å8        # Push a list with the first input `width` amount of 8s
         s       # Swap so the stack order is: top-part, x-axis, bottom-part
          )      # Wrap these three into a list
           J     # Join each inner(-most) list of digits to a single string
            ˜    # Flatten it to a single list of string lines
             0ð: # And replace all 0s with spaces
                 # (after which the resulting list of lines is output implicitly as result)

@Grimmy Although it works for most test cases, I'm afraid it fails for the second: current answer vs your suggestion. Your split-by method will only result in a single inner list instead of two when the bars are only positive or only negative.

Actually, your answer has incorrect rounding as well. You round away from 0, while the spec requires rounding to nearest integer (only the x.5 case is allowed to round either way). I posted my approach as a separate answer, since it's different enough.

@Grimmy I'm indeed aware the test cases are slightly different than the one in the challenge description. I actually left a comment about that on the challenge as well. The Jelly answer has (or maybe had - haven't checked his latest version) the same issues as my answer by not rounding at all and simply using the floats for the checks. Regardless, nice 26 byter! :)

Haskell, 324 312 298 bytes
import Data.List
l=genericLength
(!)=replicate
0#_=[]
w#x=take n x:(w-1)#(drop n x)where n=ceiling(l x/w)
d x=transpose$(++).(!'#')<*>(!' ').(maximum x-).max 0<$>x
w?h=((++).reverse.d<*>((round w)!'-':).d.map(*(-1))).(\x->round.(*(((h-1)/).((-).maximum<*>minimum).(0:))x)<$>x).map((/).sum<*>l).(w#)

Try it online!
A lot longer than the other answers, but I would already be happy to get it under 300 bytes. It also gives a slightly different result for test-case 2 (the 16th column is 11 high instead of 10), but I'm not sure how to fix that. I guess it's the difference in rounding methods? Haskell rounds 0.5 to even as far as I know.

05AB1E, 30 26 bytes
äÅAÐδ-àI<//òÝMαεðIиāΘyǝ}øJ

Try it online! or validate all test cases.
# scale and round the values
ä                # split the list of values in width sublists
 ÅA              # arithmetic mean of each sublist
   Ð             # triplicate
    δ-           # double-vectorized subtraction
      à          # maximum
       I</       # divide by (height - 1)
          /      # divide the list by the result
           ò     # round to nearest integer

# draw the plot
M                # maximum of the stack
 α               # absolute difference of each value with this max
  ε       }      # for each difference y:
   ðIи           #  a space character, repeated width times
      āΘ         #  list [1, 0, 0, 0, ...] (width elements)
        yǝ       #  replace characters at the indices in y
           ø     # transpose
            J    # join

Charcoal, 49 bytes
⊞θ⁰≧×∕⊖ζ⁻⌈θ⌊θθ⊟θＦη«≔÷Ｌθ⁻ηιι≦∕ΣＥι⊟θιＰ↑⁺⊘³ιＰ↓⁻⊘³ι←¹

Try it online! Link is to verbose version of code. Uses | as the drawing character. I'm not sure I've got all of the desired calculations, but I am now using the latest array splitting, which allowed me to save 9 bytes. Explanation:
⊞θ⁰≧×∕⊖ζ⁻⌈θ⌊θθ⊟θ

Temporarily add a 0 to the array so that the values can be scaled to the desired height.
Ｆη«

Loop over the width. Although this loop counts up, the array is actually processed from right-to-left.
≔÷Ｌθ⁻ηιι

Calculate how many array elements to pop.
≦∕ΣＥι⊟θι

Pop that many elements and calculate the mean.
Ｐ↑⁺⊘³ιＰ↓⁻⊘³ι←¹

Print the desired amount upwards or downwards as applicable, then print the X-axis in a leftward direction.