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

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

        let set_index_key =
            |&amp;set_index: &amp;usize| (set_descs[set_index].0.abs(), &amp;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 (&amp;set_index0, &amp;set_index1) in sorted_set_indexes.iter().tuple_windows() {
                if set_index_key(&amp;set_index0) &lt; set_index_key(&amp;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 &amp; 1) as Coeff;
            let set_ranks = &amp;set_ranks[..];
            let set_descs = &amp;set_descs[..];
            move |&amp;set_index: &amp;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 = |&amp;scol0: &amp;u32, &amp;scol1: &amp;u32| {
            scol_ranks[scol0 as usize]
                .cmp(&amp;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 (&amp;scol0, &amp;scol1) in sorted_scols.iter().tuple_windows() {
            if scol_ranks[scol0 as usize] &lt; scol_ranks[scol1 as usize] {
                next_rank += 1;
            } else if scol_cmp(&amp;scol0, &amp;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 &amp; 1 &lt;&lt; (scol &gt;&gt; 1) == 0 &amp;&amp; {
            seen |= 1 &lt;&lt; (scol &gt;&gt; 1);
            true
        }
    });
    assert_eq!(seen, !(!0 &lt;&lt; n));

    sorted_scols
}

fn count(n: u32) -&gt; Count {
    static ERR: &amp;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 = |&amp;set: &amp;u32| {
                if set.count_ones() == k {
                    !set &amp; !1 &lt;&lt; i &amp; !(!0 &lt;&lt; n) != 0
                } else {
                    set &amp; 1 &lt;&lt; i == 0
                }
            };
            let enlarge_set = |&amp;set: &amp;u32| set.count_ones() == k &amp;&amp; set &amp; 1 &lt;&lt; i == 0;
            let new_sets: Vec&lt;u32&gt; = cur_sets
                .iter()
                .copied()
                .filter(keep_set)
                .chain(
                    cur_sets
                        .iter()
                        .copied()
                        .filter(enlarge_set)
                        .map(|set| set | 1 &lt;&lt; i),
                )
                .collect();
            let mut new_set_indexes = vec![!0; 1 &lt;&lt; n];
            for (set_index, &amp;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&lt;Coeff&gt; = cur_sets
                        .iter()
                        .copied()
                        .zip(&amp;cur_state)
                        .filter(|(set, _)| keep_set(set))
                        .map(|(_, &amp;coeff)| coeff)
                        .chain(
                            cur_sets
                                .iter()
                                .copied()
                                .zip(&amp;cur_state)
                                .filter(|(set, _)| enlarge_set(set))
                                .map(|(set, &amp;coeff)| {
                                    let mut new_coeff = coeff * elt;
                                    if set &amp; !(!0 &lt;&lt; i) != 0 {
                                        new_coeff +=
                                            cur_state[cur_set_indexes[set as usize | 1 &lt;&lt; i]];
                                    }
                                    new_coeff
                                }),
                        )
                        .collect();

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

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

                    let new_state: Vec&lt;Coeff&gt; = new_sets
                        .iter()
                        .map(|&amp;set| {
                            let mut permuted_set = 0;
                            let mut sign = 0;
                            for (col, &amp;permuted_scol) in sorted_scols.iter().enumerate() {
                                sign ^= set &gt;&gt; col &amp; permuted_scol &amp; 1;
                                permuted_set |= (set &gt;&gt; col &amp; 1) &lt;&lt; (permuted_scol &gt;&gt; 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(&amp;vec![0]).unwrap_or(&amp;0)
}

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


Cargo.toml

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