Count the number of Hankelable matrices

Question

Background

A binary Hankel matrix is a matrix with constant skew-diagonals (positive sloping diagonals) containing only 0s and 1s. E.g. a 5x5 binary Hankel matrix looks like

a b c d e
b c d e f
c d e f g
d e f g h
e f g h i

where a, b, c, d, e, f, g, h, i are either 0 or 1.

Let's define a matrix M as Hankelable if there is a permutation of the order of rows and columns of M so that M is a Hankel matrix. This means one can apply one permutation to the order of the rows and a possibly different one to the columns.

The Challenge

The challenge is to count how many Hankelable n by n matrices there are for all n up to as large a value as possible.

Output

For each integer n from 1 upwards, output the number of Hankelable n by n matrices with entries which are 0 or 1.

For n = 1,2,3,4,5 the answers should be 2,12,230,12076,1446672. (Thanks to orlp for the code to produce these.)

Time limit

I will run your code on my machine and stop it after 1 minute. The code that outputs the correct answers up to the largest value of n wins. The time limit is for everything from n = 1 to the biggest value of n for which you give an answer.

The winner will be the best answer by the end of Saturday 18 April.

Tie Breaker

In the case of a tie for some maximum n I will time how long it takes to give the outputs up to n+1 and the quickest one wins. In the case that they run in the same time to within a second up to n+1, the first submission wins.

Languages and libraries

You can use any language which has a freely available compiler/interpreter/etc. for Linux and any libraries which are also freely available for Linux.

My machine

The timings will be run on my machine. This is a standard ubuntu install on an AMD FX-8350 Eight-Core Processor on an Asus M5A78L-M/USB3 Motherboard (Socket AM3+, 8GB DDR3). This also means I need to be able to run your code. As a consequence, only use easily available free software and please include full instructions how to compile and run your code.

Notes

I recommend against iterating over all n by n matrices and trying to detect if each one has the property I describe. First, there are too many and second, there seems to be no quick way to do this detection.

Leading entries so far

• n = 8 by Peter Taylor. Java
• n = 5 by orlp. Python

Answer 1

Java (n=8)

import java.util.*;
import java.util.concurrent.*;

public class HankelCombinatorics {
    public static final int NUM_THREADS = 8;

    private static final int[] FACT = new int[13];
    static {
        FACT[0] = 1;
        for (int i = 1; i < FACT.length; i++) FACT[i] = i * FACT[i-1];
    }

    public static void main(String[] args) {
        long prevElapsed = 0, start = System.nanoTime();
        for (int i = 1; i < 12; i++) {
            long count = count(i), elapsed = System.nanoTime() - start;
            System.out.format("%d in %dms, total elapsed %dms\n", count, (elapsed - prevElapsed) / 1000000, elapsed / 1000000);
            prevElapsed = elapsed;
        }
    }

    @SuppressWarnings("unchecked")
    private static long count(int n) {
        int[][] perms = new int[FACT[n]][];
        genPermsInner(0, 0, new int[n], perms, 0);

        // We partition by canonical representation of the row sum multiset, discarding any with a density > 50%.
        Map<CanonicalMatrix, Map<CanonicalMatrix, Integer>> part = new HashMap<CanonicalMatrix, Map<CanonicalMatrix, Integer>>();
        for (int m = 0; m < 1 << (2*n-1); m++) {
            int density = 0;
            int[] key = new int[n];
            for (int i = 0; i < n; i++) {
                key[i] = Integer.bitCount((m >> i) & ((1 << n) - 1));
                density += key[i];
            }
            if (2 * density <= n * n) {
                CanonicalMatrix _key = new CanonicalMatrix(key);
                Map<CanonicalMatrix, Integer> map = part.get(_key);
                if (map == null) part.put(_key, map = new HashMap<CanonicalMatrix, Integer>());
                map.put(new CanonicalMatrix(m, perms[0]), m);
            }
        }

        List<Job> jobs = new ArrayList<Job>();
        ExecutorService pool = Executors.newFixedThreadPool(NUM_THREADS);

        for (Map.Entry<CanonicalMatrix, Map<CanonicalMatrix, Integer>> e : part.entrySet()) {
            Job job = new Job(n, perms, e.getKey().sum() << 1 == n * n ? 0 : 1, e.getValue());
            jobs.add(job);
            pool.execute(job);
        }

        pool.shutdown();
        try {
            pool.awaitTermination(1, TimeUnit.DAYS); // i.e. until it's finished - inaccurate results are useless
        }
        catch (InterruptedException ie) {
            throw new IllegalStateException(ie);
        }

        long total = 0;
        for (Job job : jobs) total += job.subtotal;
        return total;
    }

    private static int genPermsInner(int idx, int usedMask, int[] a, int[][] perms, int off) {
        if (idx == a.length) perms[off++] = a.clone();
        else for (int i = 0; i < a.length; i++) {
            int m = 1 << (a[idx] = i);
            if ((usedMask & m) == 0) off = genPermsInner(idx+1, usedMask | m, a, perms, off);
        }
        return off;
    }

    static class Job implements Runnable {
        private volatile long subtotal = 0;
        private final int n;
        private final int[][] perms;
        private final int shift;
        private final Map<CanonicalMatrix, Integer> unseen;

        public Job(int n, int[][] perms, int shift, Map<CanonicalMatrix, Integer> unseen) {
            this.n = n;
            this.perms = perms;
            this.shift = shift;
            this.unseen = unseen;
        }

        public void run() {
            long result = 0;
            int[][] perms = this.perms;
            Map<CanonicalMatrix, Integer> unseen = this.unseen;
            while (!unseen.isEmpty()) {
                int m = unseen.values().iterator().next();
                Set<CanonicalMatrix> equiv = new HashSet<CanonicalMatrix>();
                for (int[] perm : perms) {
                    CanonicalMatrix canonical = new CanonicalMatrix(m, perm);
                    if (equiv.add(canonical)) {
                        result += canonical.weight() << shift;
                        unseen.remove(canonical);
                    }
                }
            }

            subtotal = result;
        }
    }

    static class CanonicalMatrix {
        private final int[] a;
        private final int hash;

        public CanonicalMatrix(int m, int[] r) {
            this(permuteRows(m, r));
        }

        public CanonicalMatrix(int[] a) {
            this.a = a;
            Arrays.sort(a);

            int h = 0;
            for (int i : a) h = h * 37 + i;
            hash = h;
        }

        private static int[] permuteRows(int m, int[] perm) {
            int[] cols = new int[perm.length];
            for (int i = 0; i < perm.length; i++) {
                for (int j = 0; j < cols.length; j++) cols[j] |= ((m >> (perm[i] + j)) & 1L) << i;
            }
            return cols;
        }

        public int sum() {
            int sum = 0;
            for (int i : a) sum += i;
            return sum;
        }

        public int weight() {
            int prev = -1, count = 0, weight = FACT[a.length];
            for (int col : a) {
                if (col == prev) weight /= ++count;
                else {
                    prev = col;
                    count = 1;
                }
            }
            return weight;
        }

        @Override public boolean equals(Object obj) {
            // Deliberately unsuitable for general-purpose use, but helps catch bugs faster.
            CanonicalMatrix that = (CanonicalMatrix)obj;
            for (int i = 0; i < a.length; i++) {
                if (a[i] != that.a[i]) return false;
            }
            return true;
        }

        @Override public int hashCode() {
            return hash;
        }
    }
}

Save as HankelCombinatorics.java, compile as javac HankelCombinatorics.java, run as java -Xmx2G HankelCombinatorics.

With NUM_THREADS = 4 on my quad-core machine it gets 20420819767436 for n=8 in 50 to 55 seconds elapsed, with a fair amount of variability between runs; I expect that it should easily manage the same on your octa-core machine but will take an hour or more to get n=9.

How it works

Given n, there are 2^(2n-1) binary nxn Hankel matrices. The rows can be permuted in n! ways, and the columns in n! ways. All we need to do is to avoid double-counting...

If you calculate the sum of each row, then neither permuting the rows nor permuting the columns changes the multiset of sums. E.g.

0 1 1 0 1
1 1 0 1 0
1 0 1 0 0
0 1 0 0 1
1 0 0 1 0

has row sum multiset {3, 3, 2, 2, 2}, and so do all Hankelable matrices derived from it. This means that we can group the Hankel matrices by these row sum multisets and then handle each group independently, exploiting multiple processor cores.

There's also an exploitable symmetry: the matrices with more zeroes than ones are in bijection with the matrices with more ones than zeroes.

Double-counting occurs when Hankel matrix M_1 with row permutation r_1 and column permutation c_1 matches Hankel matrix M_2 with row permutation r_2 and column permutation c_2 (with up to two but not all three of M_1 = M_2, r_1 = r_2, c_1 = c_2). The row and column permutations are independent, so if we apply row permutation r_1 to M_1 and row permutation r_2 to M_2, the columns as multisets must be equal. So for each group, I calculate all of the column multisets obtained by applying a row permutation to a matrix in the group. The easy way to get a canonical representation of the multisets is to sort the columns, and that is also useful in the next step.

Having obtained the distinct column multisets, we need to find how many of the n! permutations of each are unique. At this point, double-counting can only occur if a given column multiset has duplicate columns: what we need to do is to count the number of occurrences of each distinct column in the multiset and then calculate the corresponding multinomial coefficient. Since the columns are sorted, it's easy to do the count.

Finally we add them all up.

The asymptotic complexity isn't trivial to calculate to full precision, because we need to make some assumptions about the sets. We evaluate on the order of 2^(2n-2) n! column multisets, taking n^2 ln n time for each (including the sorting); if grouping doesn't take more than a ln n factor, we have time complexity Theta(4^n n! n^2 ln n). But since the exponential factors completely dominate the polynomial ones, it's Theta(4^n n!) = Theta((4n/e)^n).

Answer 2

Python2/3

Pretty naive approach, in a slow language:

import itertools

def permute_rows(m):
    for perm in itertools.permutations(m):
        yield perm

def permute_columns(m):
    T = zip(*m)
    for perm in itertools.permutations(T):
        yield zip(*perm)

N = 1
while True:
    base_template = ["abcdefghijklmnopqrstuvwxyz"[i:i+N] for i in range(N)]

    templates = set()
    for c in permute_rows(base_template):
        for m in permute_columns(c):
            templates.add("".join("".join(row) for row in m))

    def possibs(free, templates):
        if free == 2*N - 1:
            return set(int(t, 2) for t in templates)

        s = set()
        for b in "01":
            new_templates = set(t.replace("abcdefghijklmnopqrstuvwxyz"[free], b) for t in templates)
            s |= possibs(free + 1, new_templates)

        return s

    print(len(possibs(0, templates)))
    N += 1

Run by typing python script.py.

Answer 3

Haskell

import Data.List
import Data.Hashable
import Control.Parallel.Strategies
import Control.Parallel
import qualified Data.HashSet as S

main = mapM putStrLn $ map (show.countHankellable) [1..]

a§b=[a!!i|i<-b]

hashNub :: (Hashable a, Eq a) => [a] -> [a]
hashNub l = go S.empty l
    where
      go _ []     = []
      go s (x:xs) = if x `S.member` s then go s xs
                                    else x : go (S.insert x s) xs

pmap = parMap rseq

makeMatrix :: Int->[Bool]->[[Bool]]
makeMatrix n vars = [vars§[i..i+n-1]|i<-[0..n-1]]

countHankellable :: Int -> Int
countHankellable n = let
    s = permutations [0..n-1]
    conjugates m = concat[permutations[r§q|r<-m]|q<-s]
    variableSets = sequence [[True,False]|x<-[0..2*(n-1)]]
 in
    length.hashNub.concat.pmap (conjugates.makeMatrix n ) $ variableSets

Nowhere near as fast as Peter's - that's a pretty impressive setup he's got there! Now with much more code copied from the internet. Usage:

$ ghc -threaded hankell.hs
$ ./hankell