Contar o número de matrizes Hankelable

12

fundo

Uma matriz binária de Hankel é uma matriz com diagonais de inclinação constantes (diagonais inclinadas positivas) contendo apenas 0s e 1s. Por exemplo, uma matriz Hankel binária 5x5 se parece com

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

Onde a, b, c, d, e, f, g, h, iestão 0ou 1.

Vamos definir uma matriz M como Hankelable se houver uma permutação da ordem das linhas e colunas de M para que M seja uma matriz de Hankel. Isso significa que é possível aplicar uma permutação na ordem das linhas e uma possivelmente diferente nas colunas.

O desafio

O desafio é contar quantos Hankelable n por nmatrizes existem para todosn até como um valor maior possível.

Resultado

Para cada número inteiro n de 1 para cima, imprima o número de Hankelablen por nmatrizes com entradas que são 0ou1 .

Para n = 1,2,3,4,5as respostas devem ser2,12,230,12076,1446672 . (Obrigado ao orlp pelo código para produzi-los.)

Prazo

Vou executar o seu código na minha máquina e pará-lo após 1 minuto. O código que gera as respostas corretas até o maior valor de n ganha. O limite de tempo é para tudo, desde n = 1o maior valor den para o qual você responde.

O vencedor será a melhor resposta até o final do sábado, 18 de abril.

Desempate

No caso de um empate por um tempo máximo n, cronometrarei quanto tempo leva para dar as saídas n+1e a mais rápida vencer. No caso em que eles sejam executados ao mesmo tempo em até um segundo n+1, a primeira submissão vence.

Línguas e bibliotecas

Você pode usar qualquer idioma que tenha um compilador / intérprete disponível gratuitamente / etc. para Linux e quaisquer bibliotecas que também estão disponíveis gratuitamente para Linux.

Minha maquina

Os horários serão executados na minha máquina. Esta é uma instalação padrão do ubuntu em um processador AMD FX-8350 de oito núcleos em uma placa-mãe Asus M5A78L-M / USB3 (soquete AM3 +, 8GB DDR3). Isso também significa que eu preciso poder executar seu código. Como conseqüência, use apenas software livre facilmente disponível e inclua instruções completas sobre como compilar e executar seu código.

Notas

Eu recomendo não repetir todas as matrizes n por n e tentar detectar se cada uma tem a propriedade que eu descrevo. Primeiro, há muitos e, segundo, parece não haver uma maneira rápida de fazer essa detecção .

Entradas principais até agora

  • n = 8 por Peter Taylor. Java
  • n = 5 por orlp. Pitão
Comunidade
fonte
4
Eu realmente gosto da palavra "Hankelable".
Alex A.
3
Para n=6o total é 260357434. Eu acho que a pressão da memória é um problema maior que o tempo da CPU.
31515 Peter
Esta é uma pergunta incrível. Eu fui completamente nerd.
Alexander-brett

Respostas:

7

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;
        }
    }
}

Salvar como HankelCombinatorics.java, compilar como javac HankelCombinatorics.java, executar como java -Xmx2G HankelCombinatorics.

Com NUM_THREADS = 4a minha máquina quad-core fica 20420819767436por n=8em 50 a 55 segundos decorridos, com uma boa quantidade de variabilidade entre as execuções; Espero que ele gerencie facilmente o mesmo em sua máquina octa-core, mas levará uma hora ou mais para obtê-lo n=9.

Como funciona

Dado n, existem matrizes 2^(2n-1)binárias nx nHankel. As linhas podem ser permutadas de n!maneiras e as colunas de n!maneiras. Tudo o que precisamos fazer é evitar a contagem dupla ...

Se você calcular a soma de cada linha, nem permutar as linhas nem permutar as colunas alterará o conjunto múltiplo de somas. Por exemplo

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

possui soma de linha multiset {3, 3, 2, 2, 2} e, assim como todas as matrizes Hankelable derivadas dela. Isso significa que podemos agrupar as matrizes de Hankel por esses multisets de soma de linhas e manipular cada grupo independentemente, explorando vários núcleos de processador.

Há também uma simetria explorável: as matrizes com mais zeros do que aquelas estão em bijeção com as matrizes com mais zeros.

Dupla contagem ocorre quando matriz de Hankel M_1com permutação de linha r_1e coluna de permutação c_1corresponde matriz de Hankel M_2com permutação de linha r_2e coluna de permutação c_2(com um máximo de dois, mas não todos os três M_1 = M_2, r_1 = r_2, c_1 = c_2). As linhas e colunas permutações são independentes, por isso, se nós aplicamos fila permutação r_1para M_1e linha de permutação r_2para M_2, as colunas como multisets devem ser iguais. Portanto, para cada grupo, calculo todos os conjuntos múltiplos de colunas obtidos aplicando uma permutação de linha a uma matriz no grupo. A maneira mais fácil de obter uma representação canônica dos multisets é classificar as colunas, o que também é útil na próxima etapa.

Tendo obtido os diversos conjuntos de colunas distintos, precisamos descobrir quantas n!permutações de cada um são únicas. Nesse ponto, a contagem dupla só pode ocorrer se um determinado conjunto múltiplo de colunas tiver colunas duplicadas: o que precisamos fazer é contar o número de ocorrências de cada coluna distinta no conjunto múltiplo e calcular o coeficiente multinomial correspondente. Como as colunas são classificadas, é fácil fazer a contagem.

Finalmente, adicionamos todos eles.

A complexidade assintótica não é trivial para calcular com precisão total, porque precisamos fazer algumas suposições sobre os conjuntos. Avaliamos a ordem dos 2^(2n-2) n!vários conjuntos de colunas, levando n^2 ln ntempo para cada um (incluindo a classificação); se o agrupamento não leva mais que um ln nfator, temos complexidade de tempo Theta(4^n n! n^2 ln n). Mas desde que os fatores exponenciais dominam completamente os polinomiais, é isso Theta(4^n n!) = Theta((4n/e)^n).

Peter Taylor
fonte
Isso é muito impressionante. Você poderia dizer algo sobre o algoritmo que você usou?
3

Python2 / 3

Abordagem bastante ingênua, em uma linguagem lenta:

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

Execute digitando python script.py.

orlp
fonte
Você tem o idioma listado como Python 2/3, mas para que ele funcione no Python 2, você não precisa from __future__ import print_function(ou algo assim)?
Alex A.
2
@AlexA. Normalmente, sim, mas não neste caso. Considere o comportamento do Python2 ao digitar return(1). Agora substitua returnpor print:)
orlp 13/04/2015
Legal! Eu aprendo algo novo todos os dias. :)
Alex A.
2

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

Nem tão rápido quanto o de Peter - é uma configuração bastante impressionante que ele conseguiu! Agora, com muito mais código copiado da internet. Uso:

$ ghc -threaded hankell.hs
$ ./hankell
alexander-brett
fonte
Uma resposta Haskell é sempre bem-vinda. Obrigado.
@Lembik - como está o meu na sua máquina?
22415 alexander-brett